Ignore:
Timestamp:
Jul 26, 2023, 10:59:36 PM (15 months ago)
Author:
evignon
Message:

modifications en vue de la reprise de l'atelier tke a la rentree

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/atke_exchange_coeff_mod.F90

    r4632 r4644  
    66
    77subroutine atke_compute_km_kh(ngrid,nlay,dtime, &
    8                         wind_u,wind_v,temp,play,pinterf, &
     8                        wind_u,wind_v,temp,play,pinterf,cdrag_uv, &
    99                        tke,Km_out,Kh_out)
    1010
    … …  
    2626
    2727
    28 USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd
    29 USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, rg, rd
    30 USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, iflag_atke_lmix, lmin
     28USE atke_turbulence_ini_mod, ONLY : iflag_atke, kappa, l0, ric, cinf, rpi, rcpd, ok_vdiff_tke
     29USE atke_turbulence_ini_mod, ONLY : cepsilon, pr_slope, pr_asym, pr_neut, ctkes,rg, rd
     30USE atke_turbulence_ini_mod, ONLY : viscom, viscoh, clmix, iflag_atke_lmix, lmin, smmin
    3131
    3232implicit none
    … …  
    4545REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
    4646REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: pinterf   ! pressure at interfaces(Pa)
    47 
     47REAL, DIMENSION(ngrid), INTENT(IN)            :: cdrag_uv   ! surface drag coefficient for momentum
    4848
    4949REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke  ! turbulent kinetic energy at interface between layers
    … …  
    7878REAL    :: lstrat     ! mixing length depending on local stratification
    7979REAL    :: taustrat   ! caracteristic timescale for turbulence in very stable conditions
     80REAL    :: netloss    ! net loss term of tke
     81REAL    :: netsource  ! net source term of tke
     82REAL    :: ustar      ! friction velocity estimation
    8083
    8184! Initializations:
    … …  
    142145            Prandtl(igrid,ilay) = -2./rpi * (pr_asym - pr_neut) * atan(Ri(igrid,ilay)/Ri1) + pr_neut
    143146        ELSE ! stable cases
    144             Sm(igrid,ilay) = max(0.,cn*(1.-Ri(igrid,ilay)/Ric))
     147            Sm(igrid,ilay) = max(smmin,cn*(1.-Ri(igrid,ilay)/Ric))
    145148            Prandtl(igrid,ilay) = pr_neut + Ri(igrid,ilay) * pr_slope
    146149            IF (Ri(igrid,ilay) .GE. Prandtl(igrid,ilay)) THEN
    … …  
    182185          l_exchange(igrid,ilay) = kappa*l0*z_interf(igrid,ilay) / (kappa*z_interf(igrid,ilay) + l0)
    183186          IF (N2(igrid,ilay) .GT. 0.) THEN
    184              lstrat=clmix*sqrt(tke(igrid,ilay))/(2.*max(sqrt(shear2(igrid,ilay)),1E-10)*(1.+sqrt(Ri(igrid,ilay))/2.))
     187             lstrat=clmix*sqrt(tke(igrid,ilay))/(2.*max(sqrt(shear2(igrid,ilay)),1.E-10)*(1.+sqrt(Ri(igrid,ilay))/2.))
    185188             IF (lstrat .LT. l_exchange(igrid,ilay)) THEN
    186189                l_exchange(igrid,ilay)=max(lstrat,lmin)
    … …  
    204207
    205208
    206 ! Computing the TKE:
    207 !===================
     209! Computing the TKE k>=2:
     210!========================
    208211IF (iflag_atke == 0) THEN
     212
     213! stationary solution (dtke/dt=0)
    209214
    210215   DO ilay=2,nlay
    … …  
    221226    DO ilay=2,nlay
    222227        DO igrid=1,ngrid
    223            qq=sqrt(2.*tke(igrid,ilay))
     228           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
    224229           delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * &
    225230                (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*shear2(igrid,ilay) &
    … …  
    230235    ENDDO
    231236
     237
    232238ELSE IF (iflag_atke == 2) THEN
    233239
    234 ! full implicit scheme resolved with a second order polynomial equation
    235 ! + exact resolution for very stable cases (iflag_atke_lmix=1)
     240! semi implicit scheme when l does not depend on tke
     241! positive-guaranteed if pr slope in stable condition >1
     242
     243   DO ilay=2,nlay
     244        DO igrid=1,ngrid
     245           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
     246           qq=(qq+l_exchange(igrid,ilay)*Sm(igrid,ilay)*dtime/sqrt(2.)      &
     247               *shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) &
     248               /(1.+qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
     249           tke(igrid,ilay)=0.5*(qq**2)
     250        ENDDO
     251    ENDDO
     252
     253
     254ELSE IF (iflag_atke == 3) THEN
     255! numerical resolution adapted from that in MAR (Deleersnijder 1992)
     256! positively defined by construction
     257
    236258    DO ilay=2,nlay
    237259        DO igrid=1,ngrid
    238            qq=sqrt(2.*tke(igrid,ilay))
    239            IF (switch_num(igrid,ilay)) THEN
    240            taustrat=clmix/2./sqrt(N2(igrid,ilay))*Sm(igrid,ilay)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &
    241                     - sqrt(N2(igrid,ilay))/2./cepsilon/clmix
    242            taustrat=-1./taustrat
    243            qq=qq*exp(-dtime/taustrat)
     260           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
     261           IF (Ri(igrid,ilay) .LT. 0.) THEN
     262              netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))
     263              netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))
    244264           ELSE
    245            delta=1.+4.*dtime/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)) * &
    246                 (qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay) * shear2(igrid,ilay) &
    247                 *(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)))
    248            qq=(-1. + sqrt(delta))/dtime*cepsilon*sqrt(2.)*l_exchange(igrid,ilay)
     265              netloss=qq/(2.*sqrt(2.)*cepsilon*l_exchange(igrid,ilay))+ &
     266                      l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*N2(igrid,ilay)/Prandtl(igrid,ilay)
     267              netsource=l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)
    249268           ENDIF
     269           qq=((qq**2)/dtime+qq*netsource)/(qq/dtime+netloss)
    250270           tke(igrid,ilay)=0.5*(qq**2)
    251271        ENDDO
    252272    ENDDO
    253273
    254 
    255 
    256 ELSE IF (iflag_atke == 3) THEN
    257 
    258 ! semi implicit scheme when l does not depend on tke
    259 ! positive-guaranteed if pr slope in stable condition >1
     274ELSE IF (iflag_atke == 4) THEN
     275! semi implicit scheme from Arpege (V. Masson methodology with
     276! Taylor expansion of the dissipation term)
    260277    DO ilay=2,nlay
    261278        DO igrid=1,ngrid
    262            qq=sqrt(2.*tke(igrid,ilay))
    263            qq=(qq+dtime*l_exchange(igrid,ilay)/sqrt(2.)*Sm(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay))) &
    264              /(1.+dtime*qq/cepsilon/l_exchange(igrid,ilay)/(2.**(3/2)))     
    265            tke(igrid,ilay)=0.5*(qq**2)
     279           qq=max(sqrt(2.*tke(igrid,ilay)),1.e-10)
     280           qq=(l_exchange(igrid,ilay)*Sm(igrid,ilay)/sqrt(2.)*shear2(igrid,ilay)*(1.-Ri(igrid,ilay)/Prandtl(igrid,ilay)) &
     281             +qq*(1.+dtime*qq/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))) &
     282             /(1.+2.*qq*dtime/(cepsilon*l_exchange(igrid,ilay)*2.*sqrt(2.)))
     283           qq=max(0.,qq)
     284           tke(igrid,ilay)=0.5*(qq**2)
    266285        ENDDO
    267286    ENDDO
    … …  
    273292
    274293END IF
     294
     295! We impose a 0 tke at nlay+1
     296!==============================
     297
     298DO igrid=1,ngrid
     299 tke(igrid,nlay+1)=0.
     300END DO
     301
     302
     303! Calculation of surface TKE (k=1)
     304!=================================
     305! surface TKE calculation inspired from what is done in Arpege (see E. Bazile note)
     306DO igrid=1,ngrid
     307 ustar=sqrt(cdrag_uv(igrid)*(wind_u(igrid,1)**2+wind_v(igrid,1)**2))
     308 tke(igrid,1)=ctkes*(ustar**2)
     309END DO
     310
     311
     312! vertical diffusion of TKE
     313!==========================
     314IF (ok_vdiff_tke) THEN
     315   CALL atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
     316ENDIF
    275317
    276318
    … …  
    292334end subroutine atke_compute_km_kh
    293335
     336!===============================================================================================
     337subroutine atke_vdiff_tke(ngrid,nlay,dtime,z_lay,z_interf,temp,play,l_exchange,Sm,tke)
     338
     339! routine that computes the vertical diffusion of TKE by the turbulence
     340! using an implicit resolution (See note by Dufresne and Ghattas
     341! E Vignon, July 2023
     342
     343USE atke_turbulence_ini_mod, ONLY : rd, cke, viscom
     344
     345
     346INTEGER, INTENT(IN) :: ngrid ! number of horizontal index (flat grid)
     347INTEGER, INTENT(IN) :: nlay  ! number of vertical index 
     348
     349REAL, INTENT(IN)    :: dtime ! physics time step (s)
     350REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: z_lay   ! altitude of mid-layers (m)
     351REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)       :: z_interf   ! altitude of bottom interfaces (m)
     352REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: temp   ! temperature (K)
     353REAL, DIMENSION(ngrid,nlay), INTENT(IN)       :: play   ! pressure (Pa)
     354REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: l_exchange     ! mixing length at interfaces between layers
     355REAL, DIMENSION(ngrid,nlay+1), INTENT(IN)     :: Sm     ! stability function for eddy diffusivity for momentum at interface between layers
     356
     357REAL, DIMENSION(ngrid,nlay+1), INTENT(INOUT)  :: tke    ! turbulent kinetic energy at interface between layers
     358
     359
     360
     361INTEGER                                       :: igrid,ilay
     362REAL, DIMENSION(ngrid,nlay+1)                 :: Ke     ! eddy diffusivity for TKE
     363REAL, DIMENSION(ngrid,nlay+1)                 :: dtke
     364REAL, DIMENSION(ngrid,nlay+1)                 :: ak, bk, ck, CCK, DDK
     365REAL                                          :: gammak,Kem,KKb,KKt
     366
     367
     368! Few initialisations
     369CCK(:,:)=0.
     370DDK(:,:)=0.
     371dtke(:,:)=0.
     372
     373
     374! Eddy diffusivity for TKE
     375
     376DO ilay=2,nlay
     377    DO igrid=1,ngrid
     378       Ke(igrid,ilay)=(viscom+l_exchange(igrid,ilay)*Sm(igrid,ilay)*sqrt(tke(igrid,ilay)))*cke
     379    ENDDO
     380ENDDO
     381! at the top of the atmosphere set to 0
     382Ke(:,nlay+1)=0.
     383! at the surface, set it equal to that at the first model level
     384Ke(:,1)=Ke(:,2)
     385
     386
     387! calculate intermediary variables
     388
     389DO ilay=2,nlay
     390    DO igrid=1,ngrid
     391    Kem=0.5*(Ke(igrid,ilay+1)+Ke(igrid,ilay))   
     392    KKt=Kem*play(igrid,ilay)/rd/temp(igrid,ilay)/(z_interf(igrid,ilay+1)-z_interf(igrid,ilay))
     393    Kem=0.5*(Ke(igrid,ilay)+Ke(igrid,ilay-1))
     394    KKb=Kem*play(igrid,ilay-1)/rd/temp(igrid,ilay-1)/(z_interf(igrid,ilay)-z_interf(igrid,ilay-1))
     395    gammak=1./(z_lay(igrid,ilay)-z_lay(igrid,ilay-1))
     396    ak(igrid,ilay)=-gammak*dtime*KKb
     397    ck(igrid,ilay)=-gammak*dtime*KKt
     398    bk(igrid,ilay)=1.+gammak*dtime*(KKt+KKb)
     399    ENDDO
     400ENDDO
     401
     402! calculate CCK and DDK coefficients
     403! downhill phase
     404
     405DO igrid=1,ngrid
     406  CCK(igrid,nlay)=tke(igrid,nlay)/bk(igrid,nlay)
     407  DDK(igrid,nlay)=-ak(igrid,nlay)/bk(igrid,nlay)
     408ENDDO
     409
     410
     411DO ilay=nlay-1,2,-1
     412    DO igrid=1,ngrid
     413        CCK(igrid,ilay)=(tke(igrid,ilay)/bk(igrid,ilay)-ck(igrid,ilay)/bk(igrid,ilay)*CCK(igrid,ilay+1)) &
     414                       / (1.+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
     415        DDK(igrid,ilay)=-ak(igrid,ilay)/bk(igrid,ilay)/(1+ck(igrid,ilay)/bk(igrid,ilay)*DDK(igrid,ilay+1))
     416    ENDDO
     417ENDDO
     418
     419! calculate TKE
     420! uphill phase
     421
     422DO ilay=2,nlay+1
     423    DO igrid=1,ngrid
     424        dtke(igrid,ilay)=CCK(igrid,ilay)+DDK(igrid,ilay)*tke(igrid,ilay-1)-tke(igrid,ilay)
     425    ENDDO
     426ENDDO
     427
     428! update TKE
     429tke(:,:)=tke(:,:)+dtke(:,:)
     430
     431
     432end subroutine atke_vdiff_tke
     433
     434
    294435
    295436end module atke_exchange_coeff_mod
Note: See TracChangeset for help on using the changeset viewer.