Changeset 3597 for palm/trunk/SOURCE/land_surface_model_mod.f90

Timestamp:

Dec 4, 2018 8:40:18 AM (6 years ago)

Author:

maronga

Message:

revised calculation of near surface air potential temperature

File:

• palm/trunk/SOURCE/land_surface_model_mod.f90 (modified) (15 diffs)

palm/trunk/SOURCE/land_surface_model_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Added pt_2m / theta_2m. Removed unncessary _eb strings.
+! 
+! 3486 2018-11-05 06:20:18Z maronga
 ! Bugfix for liquid water treatment on pavements
 ! 
@@ -490 +493 @@
 
     USE basic_constants_and_equations_mod,                                     &
-        ONLY:  c_p, g, lv_d_cp, l_v, magnus, rho_l, r_d, r_v, rd_d_rv
+        ONLY:  c_p, g, lv_d_cp, l_v, kappa, magnus, rho_l, r_d, r_v, rd_d_rv
 
     USE calc_mean_profile_mod,                                                 &
@@ -530 +533 @@
 
     USE surface_mod,                                                           &
-        ONLY :  ind_pav_green, ind_veg_wall, ind_wat_win, surf_lsm_h,          &
-                surf_lsm_v, surf_type, surface_restore_elements
+        ONLY :  ind_pav_green, ind_veg_wall, ind_wat_win,                      &
+                surf_lsm_h, surf_lsm_v, surf_type, surface_restore_elements
 
     IMPLICIT NONE
@@ -750 +753 @@
               c_soil_av,        & !< average of c_soil
               c_veg_av,         & !< average of c_veg
-              lai_av,           & !< average of lai
+              lai_av,           & !< average of lai       
               qsws_liq_av,      & !< average of qsws_liq
               qsws_soil_av,     & !< average of qsws_soil
               qsws_veg_av,      & !< average of qsws_veg
               r_s_av              !< average of r_s
-                   
-
+  
 !
 !-- Predefined Land surface classes (vegetation_type)
@@ -1000 +1002 @@
 !-- Public functions
     PUBLIC lsm_boundary_condition, lsm_check_data_output,                      &
-           lsm_check_data_output_pr,                                           &
+           lsm_check_data_output_pr, lsm_calc_pt_near_surface,                 &
            lsm_check_parameters, lsm_define_netcdf_grid, lsm_3d_data_averaging,& 
            lsm_data_output_2d, lsm_data_output_3d, lsm_energy_balance,         &
@@ -1021 +1023 @@
     END INTERFACE lsm_boundary_condition
 
+    INTERFACE lsm_calc_pt_near_surface
+       MODULE PROCEDURE lsm_calc_pt_near_surface
+    END INTERFACE lsm_calc_pt_near_surface
+    
     INTERFACE lsm_check_data_output
        MODULE PROCEDURE lsm_check_data_output
@@ -2340 +2346 @@
     IF ( horizontal  .AND.  .NOT. constant_roughness )  CALL calc_z0_water_surface
     
+    
     CONTAINS
 !------------------------------------------------------------------------------!
@@ -2391 +2398 @@
                      
        ENDDO
-
+  
     END SUBROUTINE calc_q_surface
-
-
-
+       
  END SUBROUTINE lsm_energy_balance
-
+    
+   
 
 !------------------------------------------------------------------------------!
@@ -4923 +4929 @@
        ALLOCATE ( surf_lsm_h%r_s(1:surf_lsm_h%ns)                 )
        ALLOCATE ( surf_lsm_h%r_canopy_min(1:surf_lsm_h%ns)        )
+       ALLOCATE ( surf_lsm_h%pt_2m(1:surf_lsm_h%ns)               )
        ALLOCATE ( surf_lsm_h%vegetation_surface(1:surf_lsm_h%ns)  )
        ALLOCATE ( surf_lsm_h%water_surface(1:surf_lsm_h%ns)       )
@@ -5742 +5749 @@
                 ENDDO
              ENDIF
-
+             
           CASE DEFAULT
              CONTINUE
@@ -6090 +6097 @@
 
           IF ( mode == 'xy' ) grid = 'zs'
-
+         
        CASE ( 'qsws_liq*_xy' )        ! 2d-array
           IF ( av == 0 ) THEN
@@ -6208 +6215 @@
           IF ( mode == 'xy' )  grid = 'zs'
 
+
        CASE DEFAULT
           found = .FALSE.
@@ -6377 +6385 @@
        WRITE ( 14 )  qsws_veg_av
     ENDIF
-
+    
     IF ( ALLOCATED( t_soil_av ) )  THEN
        CALL wrd_write_string( 't_soil_av' )
@@ -6676 +6684 @@
               t_soil_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =           &
                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
-
+                 
            CASE ( 'lsm_start_index_h', 'lsm_start_index_v'  )   
                 IF ( k == 1 )  THEN
@@ -7196 +7204 @@
     END SUBROUTINE calc_z0_water_surface
 
-
-
+    
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculates 2m temperature for data output at coarse resolution
+!------------------------------------------------------------------------------!
+    SUBROUTINE lsm_calc_pt_near_surface
+
+       IMPLICIT NONE
+
+       INTEGER(iwp)                          :: i, j, k, l, m   !< running indices
+
+
+       DO  m = 1, surf_lsm_h%ns
+
+          i = surf_lsm_h%i(m)
+          j = surf_lsm_h%j(m)
+          k = surf_lsm_h%k(m)
+
+          surf_lsm_h%pt_2m(m) = surf_lsm_h%pt_surface(m) + surf_lsm_h%ts(m) / kappa &
+                             * ( log( 2.0_wp /  surf_lsm_h%z0h(m) )                 &
+                               - psi_h( 2.0_wp / surf_lsm_h%ol(m) )                 &
+                               + psi_h( surf_lsm_h%z0h(m) / surf_lsm_h%ol(m) ) )
+
+       ENDDO
+
+    END SUBROUTINE lsm_calc_pt_near_surface
+
+    
+    
+!
+!-- Integrated stability function for heat and moisture
+    FUNCTION psi_h( zeta )
+
+           USE kinds
+
+       IMPLICIT NONE
+
+       REAL(wp)            :: psi_h !< Integrated similarity function result
+       REAL(wp)            :: zeta  !< Stability parameter z/L
+       REAL(wp)            :: x     !< dummy variable
+
+       REAL(wp), PARAMETER :: a = 1.0_wp            !< constant
+       REAL(wp), PARAMETER :: b = 0.66666666666_wp  !< constant
+       REAL(wp), PARAMETER :: c = 5.0_wp            !< constant
+       REAL(wp), PARAMETER :: d = 0.35_wp           !< constant
+       REAL(wp), PARAMETER :: c_d_d = c / d         !< constant
+       REAL(wp), PARAMETER :: bc_d_d = b * c / d    !< constant
+
+
+      IF ( zeta < 0.0_wp )  THEN
+         x = SQRT( 1.0_wp  - 16.0_wp * zeta )
+         psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp )
+      ELSE
+         psi_h = - b * ( zeta - c_d_d ) * EXP( -d * zeta ) - (1.0_wp          &
+                 + 0.66666666666_wp * a * zeta )**1.5_wp - bc_d_d             &
+                 + 1.0_wp
+!
+!--       Old version for stable conditions (only valid for z/L < 0.5)
+!--       psi_h = - 5.0_wp * zeta
+       ENDIF
+
+   END FUNCTION psi_h
+    
  END MODULE land_surface_model_mod