Ignore:
Timestamp:
Mar 10, 2016 11:01:04 AM (6 years ago)
Author:
maronga
Message:

added support for water and paved surfaced in land surface model / minor changes

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/surface_layer_fluxes.f90

    r1758 r1788  
    1919! Current revisions:
    2020! ------------------
    21 !
     21! Added parameter z0q which replaces z0h in the similarity functions for
     22! humidity.
     23! Syntax layout improved.
    2224!
    2325! Former revisions:
     
    117119    USE arrays_3d,                                                             &
    118120        ONLY:  e, kh, nr, nrs, nrsws, ol, pt, q, ql, qr, qrs, qrsws, qs, qsws, &
    119                shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h
     121               shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h, z0q
    120122
    121123    USE cloud_parameters,                                                      &
     
    234236!--    Use either Newton iteration or a lookup table for the bulk Richardson
    235237!--    number to calculate the Obukhov length
    236        ELSEIF ( most_method == 'newton' .OR. most_method == 'lookup' )  THEN
     238       ELSEIF ( most_method == 'newton'  .OR. most_method == 'lookup' )  THEN
    237239
    238240          CALL calc_uv_total
     
    302304!--       Check for roughness heterogeneity. In that case terminate run and
    303305!--       inform user
    304           IF ( MINVAL( z0h ) /= MAXVAL( z0h ) .OR.                             &
     306          IF ( MINVAL( z0h ) /= MAXVAL( z0h )  .OR.                            &
    305307               MINVAL( z0  ) /= MAXVAL( z0  ) )  THEN
    306308             terminate_run_l = .TRUE.
     
    347349!
    348350!--       Calculate stretching factor for the free convection range
    349           DO WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
     351          DO  WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
    350352             regr_old = regr
    351353             regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr )   &
     
    533535!--             Ensure that the bulk Richardson number and the Obukhov
    534536!--             lengtH have the same sign
    535                 IF ( rib(j,i) * ol(j,i) < 0.0_wp .OR.                        &
     537                IF ( rib(j,i) * ol(j,i) < 0.0_wp  .OR.                         &
    536538                     ABS( ol(j,i) ) == ol_max )  THEN
    537539                   IF ( rib(j,i) > 0.0_wp ) ol(j,i) =  0.01_wp
     
    638640
    639641          !$OMP PARALLEL DO PRIVATE( k, z_mo )
    640           !# WARNING: does not work on GPU so far because of DO WHILE construct
     642          !# WARNING: does not work on GPU so far because of DO  WHILE construct
    641643          !!!!!!$acc kernels loop
    642644          DO  i = nxl, nxr
     
    656658                l = l_bnd
    657659                IF ( rib_tab(l) - rib(j,i) > 0.0_wp )  THEN
    658                    DO WHILE ( rib_tab(l-1) - rib(j,i) > 0.0_wp .AND. l > 0 )
     660                   DO  WHILE ( rib_tab(l-1) - rib(j,i) > 0.0_wp  .AND. l > 0 )
    659661                      l = l-1
    660662                   ENDDO
    661663                ELSE
    662                    DO WHILE ( rib_tab(l) - rib(j,i) < 0.0_wp                &
    663                               .AND. l < num_steps-1 )
     664                   DO  WHILE ( rib_tab(l) - rib(j,i) < 0.0_wp                  &
     665                              .AND.  l < num_steps-1 )
    664666                      l = l+1
    665667                   ENDDO
     
    669671!
    670672!--             Linear interpolation to find the correct value of z/L
    671                 ol(j,i) = ( ol_tab(l-1) + ( ol_tab(l) - ol_tab(l-1) )       &
    672                             / (  rib_tab(l) - rib_tab(l-1) )                &
     673                ol(j,i) = ( ol_tab(l-1) + ( ol_tab(l) - ol_tab(l-1) )          &
     674                            / (  rib_tab(l) - rib_tab(l-1) )                   &
    673675                            * ( rib(j,i) - rib_tab(l-1) ) )
    674676
     
    797799!
    798800!--       For a given surface temperature:
    799           IF ( large_scale_forcing .AND. lsf_surf )  THEN
     801          IF ( large_scale_forcing  .AND. lsf_surf )  THEN
    800802             !$OMP PARALLEL DO
    801803             !$acc kernels loop private( j, k )
     
    837839          IF ( constant_waterflux )  THEN
    838840!
    839 !--          For a given water flux in the Prandtl layer:
     841!--          For a given water flux in the surface layer
    840842             !$OMP PARALLEL DO
    841843             !$acc kernels loop private( j )
     
    847849
    848850          ELSE
    849              coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND.      &
     851             coupled_run = ( coupling_mode == 'atmosphere_to_ocean'  .AND.     &
    850852                             run_coupled )
    851853
    852              IF ( large_scale_forcing .AND. lsf_surf )  THEN
     854             IF ( large_scale_forcing  .AND. lsf_surf )  THEN
    853855                !$OMP PARALLEL DO
    854856                !$acc kernels loop private( j, k )
     
    862864
    863865             !$OMP PARALLEL DO PRIVATE( e_s, k, z_mo )
    864              !$acc kernels loop independent present( nzb_s_inner, ol, pt, q, qs, qv1, zu, zw, z0h ) private( e_s, j, k, z_mo )
     866             !$acc kernels loop independent present( nzb_s_inner, ol, pt, q, qs, qv1, zu, zw, z0q ) private( e_s, j, k, z_mo )
    865867             DO  i = nxlg, nxrg
    866868                !$acc loop independent
     
    882884                   IF ( cloud_physics )  THEN
    883885                      qs(j,i) = kappa * ( qv1(j,i) - q(k,j,i) )                &
    884                                         / ( LOG( z_mo / z0h(j,i) )             &
     886                                        / ( LOG( z_mo / z0q(j,i) )             &
    885887                                            - psi_h( z_mo / ol(j,i) )          &
    886                                             + psi_h( z0h(j,i) / ol(j,i) ) )
     888                                            + psi_h( z0q(j,i) / ol(j,i) ) )
    887889
    888890                   ELSE
    889891                      qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) )              &
    890                                         / ( LOG( z_mo / z0h(j,i) )             &
     892                                        / ( LOG( z_mo / z0q(j,i) )             &
    891893                                            - psi_h( z_mo / ol(j,i) )          &
    892                                             + psi_h( z0h(j,i) / ol(j,i) ) )
     894                                            + psi_h( z0q(j,i) / ol(j,i) ) )
    893895                   ENDIF
    894896
     
    901903!
    902904!--    If required compute qr* and nr*
    903        IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation )  THEN
     905       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )   &
     906       THEN
    904907
    905908          !$OMP PARALLEL DO PRIVATE( k, z_mo )
    906           !$acc kernels loop independent present( nr, nrs, nzb_s_inner, ol, qr, qrs, zu, zw, z0h ) private( j, k, z_mo )
     909          !$acc kernels loop independent present( nr, nrs, nzb_s_inner, ol, qr, qrs, zu, zw, z0q ) private( j, k, z_mo )
    907910          DO  i = nxlg, nxrg
    908911             !$acc loop independent
     
    913916
    914917                qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) )              &
    915                                  / ( LOG( z_mo / z0h(j,i) )                 &
     918                                 / ( LOG( z_mo / z0q(j,i) )                 &
    916919                                     - psi_h( z_mo / ol(j,i) )              &
    917                                      + psi_h( z0h(j,i) / ol(j,i) ) )
     920                                     + psi_h( z0q(j,i) / ol(j,i) ) )
    918921
    919922                nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) )              &
    920                                  / ( LOG( z_mo / z0h(j,i) )                 &
     923                                 / ( LOG( z_mo / z0q(j,i) )                 &
    921924                                     - psi_h( z_mo / ol(j,i) )              &
    922                                      + psi_h( z0h(j,i) / ol(j,i) ) )
     925                                     + psi_h( z0q(j,i) / ol(j,i) ) )
    923926             ENDDO
    924927          ENDDO
     
    10031006!
    10041007!--    Compute the vertical kinematic heat flux
    1005        IF ( .NOT. constant_heatflux .AND. ( simulated_time <=                &
    1006             skip_time_do_lsm .OR. .NOT. land_surface ) )  THEN
     1008       IF (  .NOT.  constant_heatflux .AND.  ( simulated_time <=               &
     1009            skip_time_do_lsm  .OR.  .NOT. land_surface ) )  THEN
    10071010          !$OMP PARALLEL DO
    10081011          !$acc kernels loop independent present( shf, ts, us )
     
    10181021!
    10191022!-- Compute the vertical water/scalar flux
    1020        IF ( .NOT. constant_waterflux .AND. ( humidity .OR. passive_scalar )    &
    1021             .AND. ( simulated_time <= skip_time_do_lsm .OR. .NOT.            &
    1022             land_surface ) )  THEN
     1023       IF (  .NOT.  constant_waterflux  .AND.  ( humidity  .OR.                &
     1024             passive_scalar )  .AND.  ( simulated_time <= skip_time_do_lsm     &
     1025            .OR.  .NOT.  land_surface ) )  THEN
    10231026          !$OMP PARALLEL DO
    10241027          !$acc kernels loop independent present( qs, qsws, us )
     
    10891092
    10901093       IF ( zeta < 0.0_wp )  THEN
    1091           x = SQRT( SQRT(1.0_wp  - 16.0_wp * zeta ) )
     1094          x = SQRT( SQRT( 1.0_wp  - 16.0_wp * zeta ) )
    10921095          psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2    &
    10931096                  * ( 1.0_wp + x**2 ) * 0.125_wp )
     
    11261129
    11271130       IF ( zeta < 0.0_wp )  THEN
    1128           x = SQRT(1.0_wp  - 16.0_wp * zeta )
     1131          x = SQRT( 1.0_wp  - 16.0_wp * zeta )
    11291132          psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp )
    11301133       ELSE
Note: See TracChangeset for help on using the changeset viewer.