Changeset 3165


Ignore:
Timestamp:
Dec 27, 2023, 9:30:11 AM (11 months ago)
Author:
llange
Message:

Mars PCM

  • The computation of Cdh; Cdv is now done for each sub-grid surfaces and not only the flat one.
  • No difference with former version results when nslope = 1

LL

Location:
trunk/LMDZ.MARS
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/changelog.txt

    r3164 r3165  
    44204420== 19/12/2023 == CS
    44214421Small fix following r3163
     4422
     4423== 27/12/2023 == LL
     4424The computation of Cdh; Cdv is now done for each sub-grid surfaces and not only the flat one.
     4425

  • trunk/LMDZ.MARS/libf/phymars/vdif_cd_mod.F90

    r3154 r3165  
    1616CONTAINS
    1717
    18    SUBROUTINE vdif_cd(ngrid,nlay,pz0,pg,pz,pp,pu,pv,wstar,pts,ph,virtual,mumean,pqvap,pqsurf,pcdv,pcdh)
     18   SUBROUTINE vdif_cd(ngrid,nlay,nslope,pz0,pg,pz,pp,pu,pv,wstar,pts,ph,virtual,mumean,pqvap,pqsurf,pcdv,pcdh)
    1919   
    2020   use comcstfi_h
    … …  
    7070!   Inputs:
    7171!   ----------
    72       INTEGER, INTENT(IN) :: ngrid,nlay                     ! Number of points in the physical grid and atmospheric grid
     72      INTEGER, INTENT(IN) :: ngrid,nlay,nslope              ! Number of points in the physical grid and atmospheric grid
    7373      REAL, INTENT(IN) :: pz0(ngrid)                        ! Surface Roughness (m)
    7474      REAL, INTENT(IN) :: pg                                ! Mars gravity (m/s^2)
    7575      REAL, INTENT(IN) :: pz(ngrid,nlay),pp(ngrid,nlay)     ! Layers pseudo altitudes (m) and pressure (Pa)               
    7676      REAL, INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay)     ! Zonal and Meriditionnal  winds (m/s)
    77       REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay)         ! Surface Temperature and atmospheric temperature (K)
     77      REAL, INTENT(IN) :: pts(ngrid,nslope),ph(ngrid,nlay)  ! Surface Temperature and atmospheric temperature (K)
    7878      REAL, INTENT(IN) :: wstar(ngrid)                      ! Vertical velocity due to turbulence (m/s)
    7979      LOGICAL, INTENT(IN) :: virtual                        ! Boolean to account for vapor flottability
    8080      REAL, INTENT(IN) :: mumean(ngrid)                     ! Molecular mass of the atmosphere (kg/mol)
    8181      REAL, INTENT(IN) :: pqvap(ngrid,nlay)                 ! Water vapor mixing ratio (kg/kg) to account for vapor flottability
    82       REAL, INTENT(IN) :: pqsurf(ngrid)                     ! Water ice frost on the surface (kg/m^2) to account for vapor flottability
     82      REAL, INTENT(IN) :: pqsurf(ngrid,nslope)              ! Water ice frost on the surface (kg/m^2) to account for vapor flottability
    8383     
    8484!   Outputs:
    8585!   ----------
    86       REAL, INTENT(OUT) :: pcdv(ngrid),pcdh(ngrid)          ! momentum and heat drag coefficient (1)
     86      REAL, INTENT(OUT) :: pcdv(ngrid,nslope)               ! momentum drag coefficient (1)
     87      REAL, INTENT(OUT) :: pcdh(ngrid,nslope)               ! heat drag coefficient (1)
    8788
    8889!   Local:
    8990!   ------
    9091
    91       INTEGER ig        ! Loop variable
    92 
    93       REAL karman,nu    ! Von Karman constant and fluid kinematic viscosity
     92      INTEGER ig,islope          ! Loop variable
     93
     94      REAL karman,nu             ! Von Karman constant and fluid kinematic viscosity
    9495     
    9596      LOGICAL firstcal
    … …  
    100101!$OMP THREADPRIVATE(karman,nu)
    101102
    102       REAL rib(ngrid)        ! Bulk Richardson number (1)
    103       REAL fm(ngrid)         ! stability function for momentum (1)
    104       REAL fh(ngrid)         ! stability function for heat (1)
     103      REAL rib(ngrid)            ! Bulk Richardson number (1)
     104      REAL fm(ngrid)             ! stability function for momentum (1)
     105      REAL fh(ngrid)             ! stability function for heat (1)
    105106
    106107! Formalism for stability functions  fm= 1/phim^2; fh = 1/(phim*phih) where
    … …  
    111112      REAL betam, betah, alphah, bm, bh, lambda
    112113
    113       REAL pz0t              ! initial thermal roughness length z0t for the iterative scheme (m)
    114       REAL ric               ! critical richardson number (1)
    115       REAL reynolds(ngrid)   ! Reynolds number (1)
    116       REAL prandtl(ngrid)    ! Prandtl number  (1)
    117       REAL pz0tcomp(ngrid)   ! computed z0t during the iterative scheme (m)
    118       REAL ite               ! Number of iteration (1)
    119       REAL itemax            ! Maximum number of iteration (1)
    120       REAL residual          ! Residual between pz0t and pz0tcomp (m)
    121       REAL tol_iter          ! Tolerance for the residual to ensure convergence (1=
    122 
    123       REAL zu2(ngrid)        ! Near surface winds (m/s)
    124      
    125       REAL cdn(ngrid)        ! neutral momentum  drag coefficient (1)
    126       REAL chn(ngrid)        ! neutral heat  drag coefficient (1)
    127       REAL z1z0,z1z0t        ! ratios z1/z0 and z1/z0T
    128       REAL z1,zcd0           ! Neutral roughness (m) and Cd/Ch coefficient when call richls is not called
    129       REAL tsurf_v(ngrid)    ! Virtual surface temperature, accounting for vapor flottability
    130       REAL temp_v(ngrid)     ! Potential virtual air temperature in the frist layer, accounting for vapor flottability
    131       REAL mu_h2o            ! Molecular mass of water vapor (kg/mol)
    132       REAL tol_frost         ! Tolerance to consider the effect of frost on the surface
    133       REAL qsat(ngrid)       ! saturated mixing ratio of water vapor
     114      REAL pz0t                  ! initial thermal roughness length z0t for the iterative scheme (m)
     115      REAL ric                   ! critical richardson number (1)
     116      REAL reynolds(ngrid)       ! Reynolds number (1)
     117      REAL prandtl(ngrid)        ! Prandtl number  (1)
     118      REAL pz0tcomp(ngrid)       ! computed z0t during the iterative scheme (m)
     119      REAL ite                   ! Number of iteration (1)
     120      REAL itemax                ! Maximum number of iteration (1)
     121      REAL residual              ! Residual between pz0t and pz0tcomp (m)
     122      REAL tol_iter              ! Tolerance for the residual to ensure convergence (1=
     123
     124      REAL zu2(ngrid)            ! Near surface winds (m/s)
     125     
     126      REAL cdn(ngrid)            ! neutral momentum  drag coefficient (1)
     127      REAL chn(ngrid)            ! neutral heat  drag coefficient (1)
     128      REAL z1z0,z1z0t            ! ratios z1/z0 and z1/z0T
     129      REAL z1,zcd0               ! Neutral roughness (m) and Cd/Ch coefficient when call richls is not called
     130      REAL tsurf_v(ngrid,nslope) ! Virtual surface temperature, accounting for vapor flottability
     131      REAL temp_v(ngrid)         ! Potential virtual air temperature in the frist layer, accounting for vapor flottability
     132      REAL mu_h2o                ! Molecular mass of water vapor (kg/mol)
     133      REAL tol_frost             ! Tolerance to consider the effect of frost on the surface
     134      REAL qsat(ngrid)           ! saturated mixing ratio of water vapor
    134135           
    135136!    Code:
    … …  
    145146      pz0tcomp(:) = 0.1*pz0(:)
    146147      rib(:) = 0.
    147       pcdv(:) = 0.
    148       pcdh(:) = 0.
     148      pcdv(:,:) = 0.
     149      pcdh(:,:) = 0.
    149150     
    150151! this formulation assumes alphah=1., implying betah=betam
    … …  
    160161     
    161162      IF(virtual) then
    162          DO ig = 1,ngrid
    163             temp_v(ig) = ph(ig,1)*(1.+pqvap(ig,1)/(mu_h2o/mumean(ig)))/(1.+pqvap(ig,1))
    164             IF(pqsurf(ig).gt.tol_frost) then
    165                call watersat(1,pts(ig),pp(ig,1),qsat(ig))
    166                tsurf_v(ig) = pts(ig)*(1.+qsat(ig)/(mu_h2o/mumean(ig)))/(1.+qsat(ig))
    167             ELSE
    168                tsurf_v(ig) = pts(ig)*(1.+pqvap(ig,1)/(mu_h2o/mumean(ig)))/(1.+pqvap(ig,1))
    169             ENDIF
     163         temp_v(:) = ph(:,1)*(1.+pqvap(:,1)/(mu_h2o/mumean(:)))/(1.+pqvap(:,1))
     164         DO islope = 1,nslope
     165            DO ig = 1,ngrid
     166               IF(pqsurf(ig,islope).gt.tol_frost) then
     167                  call watersat(1,pts(ig,islope),pp(ig,1),qsat(ig))
     168                  tsurf_v(ig,islope) = pts(ig,islope)*(1.+qsat(ig)/(mu_h2o/mumean(ig)))/(1.+qsat(ig))
     169               ELSE
     170                  tsurf_v(ig,islope) = pts(ig,islope)*(1.+pqvap(ig,1)/(mu_h2o/mumean(ig)))/(1.+pqvap(ig,1))
     171               ENDIF
     172            ENDDO
    170173         ENDDO
    171174      ELSE
    172          tsurf_v(:) = pts(:)
     175         tsurf_v(:,:) = pts(:,:)
    173176         temp_v(:) =  ph(:,1)
    174177      ENDIF
    … …  
    177180! Original formulation as in LMDZ Earth:  Cd = Ch = (kappa/(ln(1+z1/z0)))^2
    178181! NB: In forget et al., 1999, it's Cd = Ch = (kappa/(ln(z1/z0)))^2
    179          DO ig=1,ngrid
    180             z1=1.E+0 + pz(ig,1)/pz0(ig)
    181             zcd0=karman/log(z1)
    182             zcd0=zcd0*zcd0
    183             pcdv(ig)=zcd0
    184             pcdh(ig)=zcd0
     182         DO ig = 1,ngrid
     183            z1 = 1.E+0 + pz(ig,1)/pz0(ig)
     184            zcd0 = karman/log(z1)
     185            zcd0 = zcd0*zcd0
     186            DO islope = 1,nslope
     187               pcdv(ig,islope) = zcd0
     188               pcdh(ig,islope) = zcd0
     189            ENDDO
    185190         ENDDO
    186191      ELSE 
    187192! We follow the parameterization from Colaitis et al., 2013 (supplementary material)
    188          DO ig=1,ngrid
    189             ite = 0.
    190             residual = abs(pz0tcomp(ig)-pz0t)
    191             z1z0=pz(ig,1)/pz0(ig)
    192             cdn(ig)=karman/log(z1z0)
    193             cdn(ig)=cdn(ig)*cdn(ig)
     193         DO islope = 1,nslope
     194            DO ig=1,ngrid
     195               ite = 0.
     196               residual = abs(pz0tcomp(ig)-pz0t)
     197               z1z0=pz(ig,1)/pz0(ig)
     198               cdn(ig)=karman/log(z1z0)
     199               cdn(ig)=cdn(ig)*cdn(ig)
    194200           
    195             DO WHILE((residual .gt. tol_iter*pz0(ig)) .and.  (ite .lt. itemax))
     201               DO WHILE((residual .gt. tol_iter*pz0(ig)) .and.  (ite .lt. itemax))
    196202! Computations of z0T; iterated until z0T converges 
    197                pz0t = pz0tcomp(ig)
    198                z1z0t=pz(ig,1)/pz0t
    199                chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t)
    200                IF ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then                             
    201                   zu2(ig) = pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) + (log(1.+0.7*wstar(ig) + 2.3*wstar(ig)**2))**2 ! Near surface winds account for buoyancy
    202                   IF(turb_resolved) zu2(ig)=MAX(zu2(ig),1.)
     203                  pz0t = pz0tcomp(ig)
     204                  z1z0t=pz(ig,1)/pz0t
     205                  chn(ig)=cdn(ig)*log(z1z0)/log(z1z0t)
     206                  IF ((pu(ig,1) .ne. 0.) .or. (pv(ig,1) .ne. 0.)) then                             
     207                     zu2(ig) = pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) + (log(1.+0.7*wstar(ig) + 2.3*wstar(ig)**2))**2 ! Near surface winds account for buoyancy
     208                     IF(turb_resolved) zu2(ig)=MAX(zu2(ig),1.)
    203209! Richardson number formulation proposed by D.E. England et al. (1995)
    204                   rib(ig) = (pg/tsurf_v(ig))*sqrt(pz(ig,1)*pz0(ig))*(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))*(temp_v(ig)-tsurf_v(ig))/zu2(ig)
    205                 ELSE
    206                    print*,'warning, infinite Richardson at surface'
    207                    print*,pu(ig,1),pv(ig,1)
    208                    rib(ig) = ric         
    209                 ENDIF
     210                     rib(ig) = (pg/tsurf_v(ig,islope))*sqrt(pz(ig,1)*pz0(ig))*(((log(pz(ig,1)/pz0(ig)))**2)/(log(pz(ig,1)/pz0t)))*(temp_v(ig)-tsurf_v(ig,islope))/zu2(ig)
     211                   ELSE
     212                      print*,'warning, infinite Richardson at surface'
     213                      print*,pu(ig,1),pv(ig,1)
     214                      rib(ig) = ric         
     215                   ENDIF
    210216       
    211217! Compute the stability functions fm; fh depending on the stability of the surface layer, from D.E. England et al. (95)
    212218
    213219               ! STABLE BOUNDARY LAYER :
    214                IF (rib(ig) .gt. 0.) THEN
    215                   prandtl(ig) = 1.
    216                   IF(rib(ig) .lt. ric) THEN
     220                  IF (rib(ig) .gt. 0.) THEN
     221                     prandtl(ig) = 1.
     222                     IF(rib(ig) .lt. ric) THEN
    217223               ! Assuming alphah=1. and bh=bm for stable conditions :
    218                      fm(ig) = ((ric-rib(ig))/ric)**2
    219                      fh(ig) = ((ric-rib(ig))/ric)**2
    220                    ELSE
     224                        fm(ig) = ((ric-rib(ig))/ric)**2
     225                        fh(ig) = ((ric-rib(ig))/ric)**2
     226                     ELSE
    221227               ! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
    222                      fm(ig) = 1.
    223                      fh(ig) = 1.
    224                    ENDIF
     228                        fm(ig) = 1.
     229                        fh(ig) = 1.
     230                      ENDIF
    225231               ! UNSTABLE BOUNDARY LAYER :
    226                ELSE
    227                   fm(ig) = sqrt(1.-lambda*bm*rib(ig))
    228                   fh(ig) = (1./alphah)*((1.-lambda*bh*rib(ig))**0.5)*(1.-lambda*bm*rib(ig))**0.25
    229                   prandtl(ig) = alphah*((1.-lambda*bm*rib(ig))**0.25)/((1.-lambda*bh*rib(ig))**0.5)
    230                ENDIF
     232                  ELSE
     233                     fm(ig) = sqrt(1.-lambda*bm*rib(ig))
     234                     fh(ig) = (1./alphah)*((1.-lambda*bh*rib(ig))**0.5)*(1.-lambda*bm*rib(ig))**0.25
     235                     prandtl(ig) = alphah*((1.-lambda*bm*rib(ig))**0.25)/((1.-lambda*bh*rib(ig))**0.5)
     236                  ENDIF
    231237              ! Recompute the reynolds and z0t
    232                reynolds(ig) = karman*sqrt(fm(ig))*sqrt(zu2(ig))*pz0(ig)/(log(z1z0)*nu)
    233                pz0tcomp(ig) = pz0(ig)*exp(-karman*7.3*(reynolds(ig)**0.25)*(prandtl(ig)**0.5)+5*karman)
    234                residual = abs(pz0t-pz0tcomp(ig))
    235                ite = ite+1
    236             ENDDO  ! of while
     238                  reynolds(ig) = karman*sqrt(fm(ig))*sqrt(zu2(ig))*pz0(ig)/(log(z1z0)*nu)
     239                  pz0tcomp(ig) = pz0(ig)*exp(-karman*7.3*(reynolds(ig)**0.25)*(prandtl(ig)**0.5)+5*karman)
     240                  residual = abs(pz0t-pz0tcomp(ig))
     241                  ite = ite+1
     242               ENDDO  ! of while
    237243           
    238244! Compute the coefficients Cdv; Cdh : neutral coefficient x stability functions       
    239          pcdv(ig)=cdn(ig)*fm(ig)
    240          pcdh(ig)=chn(ig)*fh(ig)         
    241          pz0t = 0. ! for next grid point
    242          ENDDO ! of ngrid
    243 
     245            pcdv(ig,islope)=cdn(ig)*fm(ig)
     246            pcdh(ig,islope)=chn(ig)*fh(ig)         
     247            pz0t = 0. ! for next grid point
     248            ENDDO ! of ngrid
     249         enddo
    244250      ENDIF !of if call richardson surface layer
    245251

  • trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F

    r3164 r3165  
    112112      REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
    113113      REAL zkq(ngrid,nlay+1)
    114       REAL zcdv(ngrid),zcdh(ngrid)
    115       REAL, INTENT(OUT) :: zcdv_true(ngrid),zcdh_true(ngrid)    ! drag coeff are used by the LES to recompute u* and hfx
     114      REAL zcdv(ngrid,nslope),zcdh(ngrid,nslope)
     115      REAL, INTENT(OUT) :: zcdv_true(ngrid,nslope)
     116      REAL, INTENT(OUT) :: zcdh_true(ngrid,nslope)    ! drag coeff are used by the LES to recompute u* and hfx
     117      REAL :: zcdv_tmp(ngrid),zcdh_tmp(ngrid)    ! drag coeffs for the major sub-grid surface
     118      REAL :: zcdv_true_tmp(ngrid),zcdh_true_tmp(ngrid)    ! drag coeffs (computed with wind gustiness for the major sub-grid surface
    116119      REAL zu(ngrid,nlay),zv(ngrid,nlay)
    117120      REAL zh(ngrid,nlay)
    … …  
    302305      zdqsdif_tot(1:ngrid)=0
    303306      h2o_ice_depth(1:ngrid,1:nslope)=1
    304      
    305307      virtual = .false.
    306308     
    … …  
    420422c       ---------------------
    421423
    422       CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pplay,pu,pv,wstar,ptsrf_tmp,
    423      &          ph,virtual,mmean(:,1),zq(:,:,igcm_h2o_vap),
    424      &          pqsurf(:,igcm_h2o_ice,major_slope(:)),
     424      CALL vdif_cd(ngrid,nlay,nslope,pz0,g,pzlay,pplay,pu,pv,wstar,
     425     &          ptsrf,ph,virtual,mmean(:,1),zq(:,:,igcm_h2o_vap),
     426     &          pqsurf(:,igcm_h2o_ice,:),
    425427     &          zcdv_true,zcdh_true)
    426428
    427         zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1)
    428 
     429      zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1)
     430
     431      DO islope = 1,nslope
    429432        IF (callrichsl) THEN
    430           zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+
     433          zcdv(:,islope)=zcdv_true(:,islope)*sqrt(zu2(:)+
    431434     &     (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2)
    432           zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+
     435          zcdh(:,islope)=zcdh_true(:,islope)*sqrt(zu2(:)+
    433436     &     (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2)
    434 
    435            ustar(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+
    436      &     (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2)
    437 
    438            tstar(:)=0.
    439            DO ig=1,ngrid
    440               IF (zcdh_true(ig) .ne. 0.) THEN      ! When Cd=Ch=0, u*=t*=0
    441                  tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig))
    442      &                     *zcdh(ig)/ustar(ig)
    443               ENDIF
    444            ENDDO
    445 
    446437        ELSE
    447            zcdv(:)=zcdv_true(:)*sqrt(zu2(:))     ! 1 / bulk aerodynamic momentum conductance
    448            zcdh(:)=zcdh_true(:)*sqrt(zu2(:))     ! 1 / bulk aerodynamic heat conductance
    449            ustar(:)=sqrt(zcdv_true(:))*sqrt(zu2(:))
    450            tstar(:)=(ph(:,1)-ptsrf_tmp(:))
    451      &      *zcdh_true(:)/sqrt(zcdv_true(:))
     438           zcdv(:,islope)=zcdv_true(:,islope)*sqrt(zu2(:))     ! 1 / bulk aerodynamic momentum conductance
     439           zcdh(:,islope)=zcdh_true(:,islope)*sqrt(zu2(:))     ! 1 / bulk aerodynamic heat conductance
    452440        ENDIF
    453 
    454 ! Some usefull diagnostics for the new surface layer parametrization :
    455 
    456 !        call write_output('vdifc_zcdv_true',
    457 !     &              'momentum drag','no units',
    458 !     &                         zcdv_true(:))
    459 !        call write_output('vdifc_zcdh_true',
    460 !     &              'heat drag','no units',
    461 !     &                         zcdh_true(:))
    462 !        call write_output('vdifc_ust',
    463 !     &              'friction velocity','m/s',ust(:))
    464 !       call write_output('vdifc_tst',
    465 !     &              'friction temperature','K',tst(:))
    466 !        call write_output('vdifc_zcdv',
    467 !     &              'aerodyn momentum conductance','m/s',
    468 !     &                         zcdv(:))
    469 !        call write_output('vdifc_zcdh',
    470 !     &              'aerodyn heat conductance','m/s',
    471 !     &                         zcdh(:))
     441      ENDDO
     442      ustar(:) = 0
     443      tstar(:) = 0
     444      DO ig = 1,ngrid
     445          zcdv_tmp(ig) = zcdv(ig,major_slope(ig))
     446          zcdh_tmp(ig) = zcdh(ig,major_slope(ig))
     447          zcdv_true_tmp(ig) = zcdv_true(ig,major_slope(ig))
     448          zcdh_true_tmp(ig) = zcdh_true(ig,major_slope(ig))
     449          IF (callrichsl) THEN
     450             ustar(ig)=sqrt(zcdv_true(ig,major_slope(ig)))
     451     &                *sqrt(zu2(ig)+(log(1.+0.7*wstar(ig) +
     452     &                 2.3*wstar(ig)**2))**2)
     453             IF (zcdh_true(ig,major_slope(ig)) .ne. 0.) THEN      ! When Cd=Ch=0, u*=t*=0
     454                tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig))
     455     &                     *zcdh_tmp(ig)/ustar(ig)
     456             ENDIF
     457          ELSE
     458                ustar(ig)=sqrt(zcdv_true(ig,major_slope(ig)))
     459     &                    *sqrt(zu2(ig))
     460                tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig))
     461     &                     *zcdh_true(ig,major_slope(ig))
     462     &                     /sqrt(zcdv_true(ig,major_slope(ig)))
     463          ENDIF
     464      ENDDO
    472465
    473466c    ** schema de diffusion turbulente dans la couche limite
    … …  
    476469
    477470       CALL vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay
    478      &              ,pu,pv,ph,zcdv_true
     471     &              ,pu,pv,ph,zcdv_true_tmp,
    479472     &              ,pq2,zkv,zkh,zq)
    480473
    … …  
    483476      pt(:,:)=ph(:,:)*ppopsk(:,:)
    484477      CALL yamada4(ngrid,nlay,nq,ptimestep,g,r,pplev,pt
    485      s   ,pzlev,pzlay,pu,pv,ph,pq,zcdv_true,pq2,zkv,zkh,zkq,ustar
     478     s   ,pzlev,pzlay,pu,pv,ph,pq,zcdv_true_tmp,pq2,zkv,zkh,zkq,ustar
    486479     s   ,9)
    487480      ENDIF
    … …  
    505498         PRINT*,'Cd for momentum and potential temperature'
    506499
    507          PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
     500         PRINT*,zcdv_tmp(ngrid/2+1),zcdh_tmp(ngrid/2+1)
    508501         PRINT*,'Mixing coefficient for momentum and pot.temp.'
    509502         DO ilev=1,nlay
    … …  
    526519
    527520      zb(1:ngrid,2:nlay)=zkv(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
    528       zb(1:ngrid,1)=zcdv(1:ngrid)*zb0(1:ngrid,1)
     521      zb(1:ngrid,1)=zcdv_tmp(1:ngrid)*zb0(1:ngrid,1)
    529522
    530523      DO ig=1,ngrid
    … …  
    601594c -------------------------------------------
    602595
    603       CALL vdif_cd(ngrid,nlay,pz0,g,pzlay,pplay,zu,zv,wstar,ptsrf_tmp,
    604      &          ph,virtual,mmean(:,1),zq(:,:,igcm_h2o_vap),
    605      &          pqsurf(:,igcm_h2o_ice,major_slope(:)),
     596      CALL vdif_cd(ngrid,nlay,nslope,pz0,g,pzlay,pplay,zu,zv,wstar,
     597     &          ptsrf,ph,virtual,mmean(:,1),zq(:,:,igcm_h2o_vap),
     598     &          pqsurf(:,igcm_h2o_ice,:),
    606599     &          zcdv_true,zcdh_true)
    607600
    608       zu2(:)=zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1)         
    609          
    610       IF (callrichsl) THEN
    611           zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+
     601      zu2(:)=zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1)
     602
     603      DO islope = 1,nslope
     604        IF (callrichsl) THEN
     605          zcdv(:,islope)=zcdv_true(:,islope)*sqrt(zu2(:)+
    612606     &     (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2)
    613           zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+
     607          zcdh(:,islope)=zcdh_true(:,islope)*sqrt(zu2(:)+
    614608     &     (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2)
    615 
    616            ustar(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+
    617      &     (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2)
    618 
    619609        ELSE
    620            zcdv(:)=zcdv_true(:)*sqrt(zu2(:))     ! 1 / bulk aerodynamic momentum conductance
    621            zcdh(:)=zcdh_true(:)*sqrt(zu2(:))     ! 1 / bulk aerodynamic heat conductance
    622            ustar(:)=sqrt(zcdv_true(:))*sqrt(zu2(:))
    623         ENDIF         
     610           zcdv(:,islope)=zcdv_true(:,islope)*sqrt(zu2(:))     ! 1 / bulk aerodynamic momentum conductance
     611           zcdh(:,islope)=zcdh_true(:,islope)*sqrt(zu2(:))     ! 1 / bulk aerodynamic heat conductance
     612        ENDIF
     613      ENDDO
     614      ustar(:) = 0
     615      tstar(:) = 0
     616      DO ig = 1,ngrid
     617          zcdv_tmp(ig) = zcdv(ig,major_slope(ig))
     618          zcdh_tmp(ig) = zcdh(ig,major_slope(ig))
     619          zcdv_true_tmp(ig) = zcdv_true(ig,major_slope(ig))
     620          zcdh_true_tmp(ig) = zcdh_true(ig,major_slope(ig))
     621          IF (callrichsl) THEN
     622             ustar(ig)=sqrt(zcdv_true(ig,major_slope(ig)))
     623     &                *sqrt(zu2(ig)+(log(1.+0.7*wstar(ig) +
     624     &                 2.3*wstar(ig)**2))**2)
     625             IF (zcdh_true(ig,major_slope(ig)) .ne. 0.) THEN      ! When Cd=Ch=0, u*=t*=0
     626                tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig))
     627     &                     *zcdh_tmp(ig)/ustar(ig)
     628             ENDIF
     629          ELSE
     630                ustar(ig)=sqrt(zcdv_true(ig,major_slope(ig)))
     631     &                    *sqrt(zu2(ig))
     632                tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig))
     633     &                     *zcdh_true(ig,major_slope(ig))
     634     &                     /sqrt(zcdv_true(ig,major_slope(ig)))
     635          ENDIF
     636      ENDDO
    624637         
    625638
    … …  
    641654c Mass variation scheme:
    642655      zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
    643       zb(1:ngrid,1)=zcdh(1:ngrid)*zb0(1:ngrid,1)
     656      zb(1:ngrid,1)=zcdh_tmp(1:ngrid)*zb0(1:ngrid,1)
    644657
    645658c on initialise dm c
    … …  
    684697
    685698      ENDIF
    686 
    687 
    688699
    689700      DO ig=1,ngrid
    … …  
    793804              zdplanck(ig) = 0.
    794805            ENDIF
    795             z1(ig)=pcapcal(ig,islope)*ptsrf(ig,islope)
     806             zb(ig,1)=zcdh(ig,islope)*zb0(ig,1)
     807             z1(ig)=pcapcal(ig,islope)*ptsrf(ig,islope)
    796808     s       + cpp*zb(ig,1)*zc(ig,1)
    797809     s       + zdplanck(ig)*ptsrf(ig,islope)
    … …  
    884896           else
    885897#endif
    886             call dustlift(ngrid,nlay,nq,rho,zcdh_true,zcdh,
     898            call dustlift(ngrid,nlay,nq,rho,zcdh_true_tmp,zcdh_tmp,
    887899     &                    pqsurf_tmp(:,igcm_co2),pdqsdif_tmp)
    888900#ifndef MESOSCALE
    … …  
    10061018     &                     /float(nsubtimestep(ig))
    10071019             if(old_wsublimation_scheme) then
    1008                zb(1:ngrid,1)=zcdv(1:ngrid)*zb0(1:ngrid,1)
     1020               zb(1:ngrid,1)=zcdv(1:ngrid,islope)*zb0(1:ngrid,1)
    10091021     &                     /float(nsubtimestep(ig))
    10101022             else
    1011                zb(1:ngrid,1)=zcdh(1:ngrid)*zb0(1:ngrid,1)
     1023               zb(1:ngrid,1)=zcdh(1:ngrid,islope)*zb0(1:ngrid,1)
    10121024     &                     /float(nsubtimestep(ig))
    10131025             endif
    … …  
    10341046             zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
    10351047             if(old_wsublimation_scheme) then
    1036                 zdqsdif_surf(ig)=rho(ig)*dryness(ig)*zcdv(ig)
     1048                zdqsdif_surf(ig)=rho(ig)*dryness(ig)*zcdv(ig,islope)
    10371049     &                           *(zq1temp(ig)-qsat(ig))
    10381050             else
    1039                 zdqsdif_surf(ig)=rho(ig)*dryness(ig)*zcdh(ig)
     1051                zdqsdif_surf(ig)=rho(ig)*dryness(ig)*zcdh(ig,islope)
    10401052     &                           *(zq1temp(ig)-qsat(ig))
    10411053             endif
    10421054
    1043             !zdqsdif_tot(ig) = zdqsdif_surf(ig)
     1055             zdqsdif_tot(ig) = zdqsdif_surf(ig)
    10441056!!! Subsurface exchange
    10451057! Check for subsurface exchanges
    … …  
    11201132              qeq(ig,1)=(ztsrf(ig)/Tice(ig,1))
    11211133     &                        *qsat2(ig,1)
    1122               resist(ig,1)=(1+(h2o_ice_depth(ig,islope)*zcdh(ig)
     1134              resist(ig,1)=(1+(h2o_ice_depth(ig,islope)*zcdh(ig,islope)
    11231135     &                           /d_coef(ig,islope)))
    11241136              !write(*,*)'R=',resist(ig,islope)
    … …  
    11451157             zd(ig,1)=zb(ig,1)*z1(ig)
    11461158             zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig)
    1147              zdqsdif_ssi(ig,1)=rho(ig)*dryness(ig)*zcdv(ig)   
     1159             zdqsdif_ssi(ig,1)=rho(ig)*dryness(ig)*zcdv(ig,islope)   
    11481160     &                      *(zq1temp(ig)-qeq(ig,islope))
    11491161                      call write_output('zdq_ssi',
Note: See TracChangeset for help on using the changeset viewer.