Changeset 3167 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Dec 30, 2023, 6:37:47 PM (12 months ago)
Author:
llange
Message:

Mars PCM
Introducing the scheme from ATKE workshop to solve turbulent diffusion + surface layer parameterization.
Works only if callatke = .true. in the callphys.def. By default, it is false and the model runs as usual with the yamada 2.5 closure
Tuning of the several parameters for the ATKE in progress
LL

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/callkeys.h

    r3062 r3167  
    1818     &   ,latentheat_surfwater,gwd_convective_source,startphy_file      &
    1919     &   ,hdo,hdofrac,cst_cap_albedo,temp_dependent_m,refill_watercap   &
    20      &   ,cloud_adapt_ts
     20     &   ,cloud_adapt_ts, callatke
    2121!$OMP THREADPRIVATE(/callkeys_l/)
    2222
    … …  
    3636     &   ,callnirco2,callnlte,callthermos,callconduct,                  &
    3737     &    calleuv,callmolvis,callmoldiff,thermochem,thermoswater        &
    38      &   ,calltherm,callrichsl,callslope,tituscap,callyamada4
     38     &   ,calltherm,callrichsl,callslope,tituscap,callyamada4,callatke
    3939
    4040      COMMON/aeroutput/dustiropacity

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

    r3144 r3167  
    243243         call getin_p("callyamada4",callyamada4)
    244244         write(*,*) " callyamada4 = ",callyamada4
    245 
     245         
     246         write(*,*) "used latest version of ATKE scheme?"
     247         callatke =.false. ! default value
     248         call getin_p("callatke",callatke)
     249         write(*,*) " callatke = ",callatke
     250
     251         if(callyamada4.and.callatke) then
     252          print*,"You can't use both yamada and atke"
     253          print*,"Choose either Yamada or ATKE"
     254          call abort_physic(modname,
     255     &     "Can't use both yamada and ATKE",1)
     256         endif
     257     
    246258         if (calltherm .and. .not.callyamada4) then
    247259          print*,'!!!! WARNING WARNING WARNING !!!!'

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

    r3151 r3167  
    2525      use write_output_mod, only: write_output
    2626      use turb_mod, only: turb_resolved
    27      
     27      use lmdz_atke_turbulence_ini, only : smmin, cinf, cepsilon, pr_slope, pr_asym, pr_neut, ri0, ri1, cn, rpi
     28
    2829      IMPLICIT NONE
    2930     
    … …  
    136137      REAL vhf(ngrid),vvv(ngrid)                                    ! vertical heat flux (W/m^2) and vertical velocity variance (m/s)
    137138      REAL Teta_out(ngrid,n_out)                                    ! Temporary variable to compute interpolated potential temperature (K)
     139      REAL sm                                                       ! Stability function computed with the ATKE scheme
     140      REAL f_ri_cd_min                                              ! minimum of the stability function with the ATKE scheme
    138141
    139142!    Code:
    … …  
    160163      itemax= 10
    161164      tol_iter =  0.01
    162      
     165      f_ri_cd_min = 0.01
     166
    163167! this formulation assumes alphah=1., implying betah=betam
    164168! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers :
    … …  
    194198            ENDIF
    195199
    196 ! Compute the stability functions fm; fh depending on the stability of the surface layer, from D.E. England et al. (95)
    197                ! STABLE BOUNDARY LAYER :
    198             IF (rib(ig) .gt. 0.) THEN
    199                prandtl(ig) = 1.
    200                IF(rib(ig) .lt. ric) THEN
    201                ! Assuming alphah=1. and bh=bm for stable conditions :
    202                   fm(ig) = ((ric-rib(ig))/ric)**2
    203                   fh(ig) = ((ric-rib(ig))/ric)**2
     200! Compute the stability functions fm; fh depending on the stability of the surface layer
     201
     202            IF(callatke) THEN
     203               ! Computation following parameterizaiton from ATKE
     204               IF(rib(ig) .gt. 0) THEN
     205                  ! STABLE BOUNDARY LAYER :
     206                  sm = MAX(smmin,cn*(1.-rib(ig)/ric))
     207                  ! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019
     208                  prandtl(ig) = pr_neut*exp(-pr_slope/pr_neut*rib(ig)+rib(ig)/pr_neut) + rib(ig) * pr_slope
     209                  fm(ig) = max((sm**(3./2.) * sqrt(cepsilon) * (1 - rib(ig) / prandtl(ig))**(1./2.)),f_ri_cd_min)
     210                  fh(ig) = max((fm(ig) / prandtl(ig)), f_ri_cd_min)
    204211               ELSE
    205                ! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
    206                   fm(ig) = 1.
    207                   fh(ig) = 1.
    208                ! should be 0 no?
    209                ENDIF
    210                ! UNSTABLE BOUNDARY LAYER :
    211             ELSE
    212                fm(ig) = sqrt(1.-lambda*bm*rib(ig))
    213                fh(ig) = (1./alphah)*((1.-lambda*bh*rib(ig))**0.5)*(1.-lambda*bm*rib(ig))**0.25
    214                prandtl(ig) = alphah*((1.-lambda*bm*rib(ig))**0.25)/((1.-lambda*bh*rib(ig))**0.5)
    215             ENDIF
    216               ! Recompute the reynolds and z0t
    217             reynolds(ig) = karman*sqrt(fm(ig))*sqrt(zu2(ig))*pz0(ig)/(log(z1z0)*nu)
    218             pz0tcomp(ig) = pz0(ig)*exp(-karman*7.3*(reynolds(ig)**0.25)*(prandtl(ig)**0.5)+5*karman)
    219             residual = abs(pz0t-pz0tcomp(ig))
    220             ite = ite+1
     212                  ! UNSTABLE BOUNDARY LAYER :
     213                  sm = 2./rpi * (cinf - cn) * atan(-rib(ig)/ri0) + cn
     214                  prandtl(ig) = -2./rpi * (pr_asym - pr_neut) * atan(rib(ig)/ri1) + pr_neut
     215                  fm(ig) = MAX((sm**(3./2.) * sqrt(cepsilon) * (1 - rib(ig) / prandtl(ig))**(1./2.)),f_ri_cd_min)
     216                  fh(ig) = MAX((fm(ig) / prandtl(ig)), f_ri_cd_min)
     217                ENDIF ! Rib < 0
     218             ELSE
     219                ! Computation following parameterizaiton from from D.E. England et al. (95)
     220                IF (rib(ig) .gt. 0.) THEN
     221                   ! STABLE BOUNDARY LAYER :
     222                   prandtl(ig) = 1.
     223                   IF(rib(ig) .lt. ric) THEN
     224                      ! Assuming alphah=1. and bh=bm for stable conditions :
     225                      fm(ig) = ((ric-rib(ig))/ric)**2
     226                      fh(ig) = ((ric-rib(ig))/ric)**2
     227                   ELSE
     228                      ! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
     229                      fm(ig) = 1.
     230                      fh(ig) = 1.
     231                   ENDIF
     232                ELSE
     233                   ! UNSTABLE BOUNDARY LAYER :
     234                   fm(ig) = sqrt(1.-lambda*bm*rib(ig))
     235                   fh(ig) = (1./alphah)*((1.-lambda*bh*rib(ig))**0.5)*(1.-lambda*bm*rib(ig))**0.25
     236                   prandtl(ig) = alphah*((1.-lambda*bm*rib(ig))**0.25)/((1.-lambda*bh*rib(ig))**0.5)
     237                ENDIF ! Rib < 0
     238             ENDIF ! atke
     239             ! Recompute the reynolds and z0t
     240             reynolds(ig) = karman*sqrt(fm(ig))*sqrt(zu2(ig))*pz0(ig)/(log(z1z0)*nu)
     241             pz0tcomp(ig) = pz0(ig)*exp(-karman*7.3*(reynolds(ig)**0.25)*(prandtl(ig)**0.5)+5*karman)
     242             residual = abs(pz0t-pz0tcomp(ig))
     243             ite = ite+1
    221244         ENDDO  ! of while
    222245! Compute the coefficients Cdv; Cdh : neutral coefficient x stability functions       

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

    r3159 r3167  
    124124      use write_output_mod, only: write_output
    125125      use pbl_parameters_mod, only: pbl_parameters
     126      use lmdz_atke_turbulence_ini, only : atke_ini
    126127      IMPLICIT NONE
    127128c=======================================================================
    … …  
    574575      REAL :: ztmp1,ztmp2                         ! intermediate variables to compute the mean molar mass of the layer
    575576
     577! Variable for the computation of the TKE with parameterization from ATKE working group
     578      REAL :: viscom                              ! kinematic molecular viscosity for momentum
     579      REAL :: viscoh                              ! kinematic molecular viscosity for heat
     580
     581
    576582c=======================================================================
    577583      pdq(:,:,:) = 0.
    … …  
    737743         ENDIF
    738744         icount=1
     745         
     746         
    739747
    740748#ifndef MESOSCALE
    … …  
    794802c        ~~~~~~~~~~~~~~~
    795803        if (topflows) call topmons_setup(ngrid)
     804       
     805c        Parameterization of the ATKE routine
     806c        ~~~~~~~~~~~~~~~
     807        if (callatke) then
     808           viscom = 0.001
     809           viscoh = 0.001
     810           CALL atke_ini(g, r, pi, cpp, 0., viscom, viscoh)
     811        endif
    796812
    797813#endif
    … …  
    24872503
    24882504      ENDIF ! refill_watercap
    2489 
    24902505
    24912506c-----------------------------------------------------------------------

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

    r3128 r3167  
    14431443
    14441444!       call write_output("coeff_diffusion_soil",'interlayer diffusion coefficient','m^2/s',D_inter)                         
    1445      
    1446      
    1447 
     1445           
    14481446!endif
    14491447endif

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

    r3165 r3167  
    2121   use turb_mod, only: turb_resolved
    2222   use watersat_mod, only: watersat
    23 
     23   use lmdz_atke_turbulence_ini, only : smmin, cinf, cepsilon, pr_slope, pr_asym, pr_neut, ri0, ri1, cn, rpi
     24   
    2425   IMPLICIT NONE
    25       
     26   include "callkeys.h"      
    2627!=======================================================================
    2728!
    … …  
    6364!   -------------
    6465
    65 #include "callkeys.h"
    66 
    6766
    6867!   Arguments:
    … …  
    134133      REAL qsat(ngrid)           ! saturated mixing ratio of water vapor
    135134           
     135      REAL sm                    ! Stability function computed with the ATKE scheme
     136      REAL f_ri_cd_min           ! minimum of the stability function with the ATKE scheme
     137     
    136138!    Code:
    137139!    --------
    … …  
    148150      pcdv(:,:) = 0.
    149151      pcdh(:,:) = 0.
    150      
     152      f_ri_cd_min = 0.01
    151153! this formulation assumes alphah=1., implying betah=betam
    152154! We use Dyer et al. parameters, as they cover a broad range of Richardson numbers :
    … …  
    215217                   ENDIF
    216218       
    217 ! Compute the stability functions fm; fh depending on the stability of the surface layer, from D.E. England et al. (95)
    218 
    219                ! STABLE BOUNDARY LAYER :
    220                   IF (rib(ig) .gt. 0.) THEN
    221                      prandtl(ig) = 1.
    222                      IF(rib(ig) .lt. ric) THEN
     219! Compute the stability functions fm; fh depending on the stability of the surface layer
     220
     221                   IF(callatke) THEN
     222                   ! Computation following parameterizaiton from ATKE
     223                      IF(rib(ig) .gt. 0) THEN
     224                         ! STABLE BOUNDARY LAYER :
     225                         sm = MAX(smmin,cn*(1.-rib(ig)/ric))
     226                         ! prandlt expression from venayagamoorthy and stretch 2010, Li et al 2019
     227                         prandtl(ig) = pr_neut*exp(-pr_slope/pr_neut*rib(ig)+rib(ig)/pr_neut) + rib(ig) * pr_slope
     228                         fm(ig) = max((sm**(3./2.) * sqrt(cepsilon) * (1 - rib(ig) / prandtl(ig))**(1./2.)),f_ri_cd_min)
     229                         fh(ig) = max((fm(ig) / prandtl(ig)), f_ri_cd_min)
     230                      ELSE
     231                         ! UNSTABLE BOUNDARY LAYER :
     232                         sm = 2./rpi * (cinf - cn) * atan(-rib(ig)/ri0) + cn
     233                         prandtl(ig) = -2./rpi * (pr_asym - pr_neut) * atan(rib(ig)/ri1) + pr_neut
     234                         fm(ig) = MAX((sm**(3./2.) * sqrt(cepsilon) * (1 - rib(ig) / prandtl(ig))**(1./2.)),f_ri_cd_min)
     235                         fh(ig) = MAX((fm(ig) / prandtl(ig)), f_ri_cd_min)
     236                      ENDIF ! Rib < 0
     237                   ELSE
     238                   ! Computation following parameterizaiton from from D.E. England et al. (95)
     239                      IF (rib(ig) .gt. 0.) THEN
     240                        ! STABLE BOUNDARY LAYER :
     241                        prandtl(ig) = 1.
     242                        IF(rib(ig) .lt. ric) THEN
    223243               ! Assuming alphah=1. and bh=bm for stable conditions :
    224                         fm(ig) = ((ric-rib(ig))/ric)**2
    225                         fh(ig) = ((ric-rib(ig))/ric)**2
     244                           fm(ig) = ((ric-rib(ig))/ric)**2
     245                           fh(ig) = ((ric-rib(ig))/ric)**2
     246                        ELSE
     247               ! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
     248                           fm(ig) = 1.
     249                           fh(ig) = 1.
     250                         ENDIF
    226251                     ELSE
    227                ! For Ri>Ric, we consider Ri->Infinity => no turbulent mixing at surface
    228                         fm(ig) = 1.
    229                         fh(ig) = 1.
    230                       ENDIF
    231                ! UNSTABLE BOUNDARY LAYER :
    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
     252                        ! UNSTABLE BOUNDARY LAYER :
     253                        fm(ig) = sqrt(1.-lambda*bm*rib(ig))
     254                        fh(ig) = (1./alphah)*((1.-lambda*bh*rib(ig))**0.5)*(1.-lambda*bm*rib(ig))**0.25
     255                        prandtl(ig) = alphah*((1.-lambda*bm*rib(ig))**0.25)/((1.-lambda*bh*rib(ig))**0.5)
     256                     ENDIF ! Rib < 0
     257                  ENDIF ! atke
    237258              ! Recompute the reynolds and z0t
    238259                  reynolds(ig) = karman*sqrt(fm(ig))*sqrt(zu2(ig))*pz0(ig)/(log(z1z0)*nu)

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

    r3166 r3167  
    3636      use comsoil_h, only: layer, mlayer,adsorption_soil
    3737      use vdif_cd_mod, only: vdif_cd
     38      use lmdz_call_atke, only: call_atke
    3839      IMPLICIT NONE
    3940
    … …  
    466467c    ** schema de diffusion turbulente dans la couche limite
    467468c       ----------------------------------------------------
    468        IF (.not. callyamada4) THEN
    469 
    470        CALL vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay
    471      &              ,pu,pv,ph,zcdv_true_tmp
    472      &              ,pq2,zkv,zkh,zq)
    473 
    474       ELSE
    475 
    476469      pt(:,:)=ph(:,:)*ppopsk(:,:)
    477       CALL yamada4(ngrid,nlay,nq,ptimestep,g,r,pplev,pt
     470      if (callyamada4) then
     471         call yamada4(ngrid,nlay,nq,ptimestep,g,r,pplev,pt
    478472     s   ,pzlev,pzlay,pu,pv,ph,pq,zcdv_true_tmp,pq2,zkv,zkh,zkq,ustar
    479      s   ,9)
    480       ENDIF
    481 
     473     s   ,9)     
     474     
     475      elseif (callatke) then
     476         call call_atke(ptimestep,ngrid,nlay,zcdv_true_tmp,
     477     s               zcdh_true_tmp,pu(:,1),pv(:,1),ptsrf_tmp,
     478     s               pu,pv,pt,zq(:,1,igcm_h2o_vap),pplay,pplev,
     479     s               pzlay,pzlev,pq2,zkv(:,1:nlay),zkh(:,1:nlay))
     480
     481         zkv(:,nlay+1) = zkv(:,nlay)
     482         zkh(:,nlay+1) = zkh(:,nlay)   
     483      else
     484        call vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay
     485     s              ,pu,pv,ph,zcdv_true_tmp
     486     s              ,pq2,zkv,zkh,zq)
     487     
     488      endif
    482489      if ((doubleq).and.(ngrid.eq.1)) then
    483490        kmixmin = 80. !80.! minimum eddy mix coeff in 1D
    484         do ilev=1,nlay
     491        do ilev=2,nlay
    485492          do ig=1,ngrid
    486493           zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev))
    … …  
    10131020            DO tsub=1,nsubtimestep(ig)
    10141021              if(tsub.eq.nsubtimestep(ig)) writeoutput = .true.
    1015 
    10161022c           C'est parti !
    10171023             zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay)
Note: See TracChangeset for help on using the changeset viewer.