Index: LMDZ5/trunk/libf/phylmd/cdrag.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/cdrag.F90 (revision 2867) +++ LMDZ5/trunk/libf/phylmd/cdrag.F90 (revision 2868) @@ -2,22 +2,26 @@ !$Id$ ! - SUBROUTINE cdrag( knon, nsrf, & + SUBROUTINE cdrag(knon, nsrf, & speed, t1, q1, zgeop1, & psol, tsurf, qsurf, z0m, z0h, & - pcfm, pcfh, zri, pref ) + cdm, cdh, zri, pref) USE dimphy USE indice_sol_mod USE print_control_mod, ONLY: lunout, prt_level + USE ioipsl_getin_p_mod, ONLY : getin_p + IMPLICIT NONE ! ================================================================= c ! ! Objet : calcul des cdrags pour le moment (pcfm) et -! les flux de chaleur sensible et latente (pcfh). -! -! Modified histroy: -! 27-Jan-2014: richardson number inconsistant between -! coefcdrag.F90 and clcdrag.F90, Fuxing WANG wrote this subroutine -! by merging coefcdrag and clcdrag. +! les flux de chaleur sensible et latente (pcfh) d'apr??s +! Louis 1982, Louis 1979, King et al 2001 +! ou Zilitinkevich et al 2002 pour les cas stables, Louis 1979 +! et 1982 pour les cas instables +! +! Modified history: +! writting on the 20/05/2016 +! modified on the 13/12/2016 to be adapted to LMDZ6 ! ! References: @@ -28,4 +32,7 @@ ! parametrization, November 1981, ECMWF, Reading, England. ! Page: 19. Equations in Table 1. +! Mascart P, Noilhan J, Giordani H 1995.A MODIFIED PARAMETERIZATION OF FLUX-PROFILE RELATIONSHIPS +! IN THE SURFACE LAYER USING DIFFERENT ROUGHNESS LENGTH VALUES FOR HEAT AND MOMENTUM +! Boundary-Layer Meteorology 72: 331-344 ! Anton Beljaars. May 1992. The parametrization of the planetary boundary layer. ! European Centre for Medium-Range Weather Forecasts. @@ -34,103 +41,121 @@ ! model to the parameterization of evaporation from the tropical oceans. J. ! Climate, 5:418-434. +! King J.C, Connolley, W.M ad Derbyshire S.H. 2001, Sensitivity of Modelled Antarctic climate +! to surface and boundary-layer flux parametrizations +! QJRMS, 127, pp 779-794 ! ! ================================================================= c -! -! knon----input-I- nombre de points pour un type de surface -! nsrf----input-I- indice pour le type de surface; voir indicesol.h -! speed---input-R- module du vent au 1er niveau du modele -! t1------input-R- temperature de l'air au 1er niveau du modele -! q1------input-R- humidite de l'air au 1er niveau du modele -! zgeop---input-R- geopotentiel au 1er niveau du modele -! psol----input-R- pression au sol -! tsurf---input-R- temperature de l'air a la surface -! qsurf---input-R- humidite de l'air a la surface -! z0m, z0h---input-R- rugosite -!! u1, v1 are removed, speed is used. Fuxing WANG, 04/03/2015, -!! u1------input-R- vent zonal au 1er niveau du modele -!! v1------input-R- vent meridien au 1er niveau du modele -! -! pcfm---output-R- cdrag pour le moment -! pcfh---output-R- cdrag pour les flux de chaleur latente et sensible -! zri----output-R- Richardson number -! pref---output-R- pression au niveau zgeop/RG -! -! Parameters: -! ckap-----Karman constant -! cb,cc,cd-parameters in Louis et al., 1982 ! ================================================================= c -! -! +! On choisit le couple de fonctions de correction avec deux flags: +! Un pour les cas instables, un autre pour les cas stables +! +! iflag_corr_insta: +! 1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h) +! 2: Louis 1982 +! 3: Laurent Li +! +! iflag_corr_sta: +! 1: Louis 1979 avec les modifications de Mascart 1995 (z0/= z0h) +! 2: Louis 1982 +! 3: Laurent Li +! 4: King 2001 (SHARP) +! 5: MO 1st order theory (allow collapse of turbulence) +! +! +!***************************************************************** ! Parametres d'entree !***************************************************************** - INTEGER, INTENT(IN) :: knon, nsrf + + INTEGER, INTENT(IN) :: knon, nsrf ! nombre de points de grille sur l'horizontal + type de surface REAL, DIMENSION(klon), INTENT(IN) :: speed ! module du vent au 1er niveau du modele REAL, DIMENSION(klon), INTENT(IN) :: zgeop1! geopotentiel au 1er niveau du modele - REAL, DIMENSION(klon), INTENT(IN) :: psol ! pression au sol - REAL, DIMENSION(klon), INTENT(IN) :: t1 ! temperature at 1st level - REAL, DIMENSION(klon), INTENT(IN) :: q1 ! humidity at 1st level REAL, DIMENSION(klon), INTENT(IN) :: tsurf ! Surface temperature (K) REAL, DIMENSION(klon), INTENT(IN) :: qsurf ! Surface humidity (Kg/Kg) REAL, DIMENSION(klon), INTENT(IN) :: z0m, z0h ! Rugosity at surface (m) -! paprs, pplay u1, v1: to be deleted -! they were in the old clcdrag. Fuxing WANG, 04/03/2015 -! REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs -! REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay -! REAL, DIMENSION(klon), INTENT(IN) :: u1, v1 + REAL, DIMENSION(klon), INTENT(IN) :: t1 ! Temperature au premier niveau (K) + REAL, DIMENSION(klon), INTENT(IN) :: q1 ! humidite specifique au premier niveau (kg/kg) + REAL, DIMENSION(klon), INTENT(IN) :: psol ! pression au sol + + ! Parametres de sortie !****************************************************************** - REAL, DIMENSION(klon), INTENT(OUT) :: pcfm ! Drag coefficient for heat flux - REAL, DIMENSION(klon), INTENT(OUT) :: pcfh ! Drag coefficient for momentum + REAL, DIMENSION(klon), INTENT(OUT) :: cdm ! Drag coefficient for heat flux + REAL, DIMENSION(klon), INTENT(OUT) :: cdh ! Drag coefficient for momentum REAL, DIMENSION(klon), INTENT(OUT) :: zri ! Richardson number REAL, DIMENSION(klon), INTENT(OUT) :: pref ! Pression au niveau zgeop/RG -! Parametres local - INTEGER :: ng_q1 ! Number of grids that q1 < 0.0 - INTEGER :: ng_qsurf ! Number of grids that qsurf < 0.0 -! zgeop1, psol: to be deleted, they are inputs now. Fuxing WANG, 04/03/2015 -! REAL, DIMENSION(klon) :: zgeop1! geopotentiel au 1er niveau du modele -! REAL, DIMENSION(klon) :: psol ! pression au sol -! -! ================================================================= c -! +! Variables Locales +!****************************************************************** + + INCLUDE "YOMCST.h" INCLUDE "YOETHF.h" -! INCLUDE "indicesol.h" INCLUDE "clesphys.h" -! -! Quelques constantes et options: -!!$PB REAL, PARAMETER :: ckap=0.35, cb=5.0, cc=5.0, cd=5.0, cepdu2=(0.1)**2 - REAL, PARAMETER :: CKAP=0.40, CB=5.0, CC=5.0, CD=5.0, CEPDU2 = (0.1)**2 -! -! Variables locales : - INTEGER :: i - REAL :: zdu2, ztsolv - REAL :: ztvd, zscf - REAL :: zucf, zcr - REAL :: friv, frih - REAL, DIMENSION(klon) :: zcfm1, zcfm2 ! Drag coefficient for momentum - REAL, DIMENSION(klon) :: zcfh1, zcfh2 ! Drag coefficient for heat flux - LOGICAL, PARAMETER :: zxli=.FALSE. ! calcul des cdrags selon Laurent Li - REAL, DIMENSION(klon) :: zcdn_m, zcdn_h ! Drag coefficient in neutral conditions + + + REAL, PARAMETER :: CKAP=0.40, CKAPT=0.42 + REAL CEPDU2 + REAL ALPHA + REAL CB,CC,CD,C2,C3 + REAL MU, CM, CH, B, CMstar, CHstar + REAL PM, PH, BPRIME + REAL C + INTEGER ng_q1 ! Number of grids that q1 < 0.0 + INTEGER ng_qsurf ! Number of grids that qsurf < 0.0 + INTEGER i + REAL zdu2, ztsolv + REAL ztvd, zscf + REAL zucf, zcr + REAL friv, frih + REAL, DIMENSION(klon) :: FM, FH ! stability functions + REAL, DIMENSION(klon) :: cdmn, cdhn ! Drag coefficient in neutral conditions REAL zzzcd -! -! Fonctions thermodynamiques et fonctions d'instabilite - REAL :: fsta, fins, x - fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x)) - fins(x) = SQRT(1.0-18.0*x) + + LOGICAL, SAVE :: firstcall = .TRUE. + !$OMP THREADPRIVATE(firstcall) + INTEGER, SAVE :: iflag_corr_sta + !$OMP THREADPRIVATE(iflag_corr_sta) + INTEGER, SAVE :: iflag_corr_insta + !$OMP THREADPRIVATE(iflag_corr_insta) + +!===================================================================c +! Valeurs numeriques des constantes +!===================================================================c + + +! Minimum du carre du vent + + CEPDU2 = (0.1)**2 + +! Louis 1982 + + CB=5.0 + CC=5.0 + CD=5.0 + + +! King 2001 + + C2=0.25 + C3=0.0625 + + +! Louis 1979 + + BPRIME=4.7 + B=9.4 + + +!MO + + ALPHA=5.0 + ! ================================================================= c -! Fuxing WANG, 04/03/2015, delete the calculation of zgeop1 -! (le geopotentiel du premier couche de modele). -! zgeop1 is an input ivariable in this subroutine. -! DO i = 1, knon -! zgeop1(i) = RD * t1(i) / (0.5*(paprs(i,1)+pplay(i,1))) & -! * (paprs(i,1)-pplay(i,1)) -! END DO -! ================================================================= c -! +! Tests avant de commencer ! Fuxing WANG, 04/03/2015 ! To check if there are negative q1, qsurf values. +!====================================================================c ng_q1 = 0 ! Initialization ng_qsurf = 0 ! Initialization @@ -154,80 +179,233 @@ ENDIF -! Calculer le frottement au sol (Cdrag) - DO i = 1, knon -!------------------------------------------------------------ -! u1, v1 are replaced by speed. Fuxing WANG, 04/03/2015, -! zdu2 = MAX(CEPDU2,u1(i)**2+v1(i)**2) -!------------------------------------------------------------ + + +!=============================================================================c +! Calcul du cdrag +!=============================================================================c + +! On choisit les fonctions de stabilite utilisees au premier appel +!************************************************************************** + IF (firstcall) THEN + iflag_corr_sta=2 + iflag_corr_insta=2 + + CALL getin_p('iflag_corr_sta',iflag_corr_sta) + CALL getin_p('iflag_corr_insta',iflag_corr_insta) + + firstcall = .FALSE. + ENDIF + +!xxxxxxxxxxxxxxxxxxxxxxx + DO i = 1, knon ! Boucle sur l'horizontal +!xxxxxxxxxxxxxxxxxxxxxxx + + +! calculs preliminaires: +!*********************** + + zdu2 = MAX(CEPDU2, speed(i)**2) -! psol(i) = paprs(i,1) pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* & - (1.+ RETV * max(q1(i),0.0)))) ! negative q1 set to zero -!------------ the old calculations in clcdrag---------------- -! ztsolv = tsurf(i) * (1.0+RETV*qsurf(i)) -! ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) & -! *(1.+RETV*q1(i)) -!------------------------------------------------------------ -! Fuxing WANG, 04/03/2015, in this revised version, -! the negative qsurf and q1 are set to zero (as in coefcdrag) - ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0)) ! negative qsurf set to zero - ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*max(q1(i),0.0))) & - *(1.+RETV*max(q1(i),0.0)) ! negative q1 set to zero + (1.+ RETV * max(q1(i),0.0)))) ! negative q1 set to zero + ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0)) ! negative qsurf set to zero + ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) & + *(1.+RETV*max(q1(i),0.0)) ! negative q1 set to zero zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd) -! Coefficients CD neutres pour m et h : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)) +! Coefficients CD neutres : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)): +!******************************************************************** + zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i))) - zcdn_m(i) = zzzcd*zzzcd - zcdn_h(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i)))) - - IF (zri(i) .GT. 0.) THEN ! situation stable + cdmn(i) = zzzcd*zzzcd + cdhn(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i)))) + + +! Calcul des fonctions de stabilit?? FMs, FHs, FMi, FHi : +!******************************************************* + +!'''''''''''''' +! Cas instables +!'''''''''''''' + + IF (zri(i) .LT. 0.) THEN + + + SELECT CASE (iflag_corr_insta) + + CASE (1) ! Louis 1979 + Mascart 1995 + + MU=LOG(MAX(z0m(i)/z0h(i),0.01)) + CMstar=6.8741+2.6933*MU-0.3601*(MU**2)+0.0154*(MU**3) + PM=0.5233-0.0815*MU+0.0135*(MU**2)-0.001*(MU**3) + CHstar=3.2165+4.3431*MU+0.536*(MU**2)-0.0781*(MU**3) + PH=0.5802-0.1571*MU+0.0327*(MU**2)-0.0026*(MU**3) + CH=CHstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & + & * CKAPT/LOG(z0h(i)+zgeop1(i)/(RG*z0h(i))) & + & * ((zgeop1(i)/(RG*z0h(i)))**PH) + CM=CMstar*B*CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & + & *CKAP/LOG(z0m(i)+zgeop1(i)/(RG*z0m(i))) & + & * ((zgeop1(i)/(RG*z0m(i)))**PM) + + + + + FM(i)=1.-B*zri(i)/(1.+CM*SQRT(ABS(zri(i)))) + FH(i)=1.-B*zri(i)/(1.+CH*SQRT(ABS(zri(i)))) + + CASE (2) ! Louis 1982 + + zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) & + *(1.0+zgeop1(i)/(RG*z0m(i))))) + FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min) + FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min) + + + CASE (3) ! Laurent Li + + + FM(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min) + FH(i) = MAX(SQRT(1.0-18.0*zri(i)),f_ri_cd_min) + + + + CASE default ! Louis 1982 + + zucf = 1./(1.+3.0*CB*CC*cdmn(i)*SQRT(ABS(zri(i)) & + *(1.0+zgeop1(i)/(RG*z0m(i))))) + FM(i) = AMAX1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min) + FH(i) = AMAX1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min) + + + END SELECT + + + +! Calcul des drags + + + cdm(i)=cdmn(i)*FM(i) + cdh(i)=f_cdrag_ter*cdhn(i)*FH(i) + + +! Traitement particulier des cas oceaniques +! on applique Miller et al 1992 en l'absence de gustiness + + IF (nsrf == is_oce) THEN +! cdh(i)=f_cdrag_oce*cdhn(i)*FH(i) + + IF(iflag_gusts==0) THEN + zcr = (0.0016/(cdmn(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.) + cdh(i) =f_cdrag_oce* cdhn(i)*(1.0+zcr**1.25)**(1./1.25) + ENDIF + + + cdm(i)=MIN(cdm(i),cdmmax) + cdh(i)=MIN(cdh(i),cdhmax) + + END IF + + + + ELSE + +!''''''''''''''' +! Cas stables : +!''''''''''''''' zri(i) = MIN(20.,zri(i)) - IF (.NOT.zxli) THEN + + SELECT CASE (iflag_corr_sta) + + CASE (1) ! Louis 1979 + Mascart 1995 + + FM(i)=MAX(1./((1+BPRIME*zri(i))**2),f_ri_cd_min) + FH(i)=FM(i) + + + CASE (2) ! Louis 1982 + zscf = SQRT(1.+CD*ABS(zri(i))) - friv = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), f_ri_cd_min) - zcfm1(i) = zcdn_m(i) * friv - frih = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), f_ri_cd_min ) -!!$ PB zcfh1(i) = zcdn(i) * frih -!!$ PB zcfh1(i) = f_cdrag_stable * zcdn(i) * frih - zcfh1(i) = f_cdrag_ter * zcdn_h(i) * frih - IF(nsrf.EQ.is_oce) zcfh1(i) = f_cdrag_oce * zcdn_h(i) * frih -!!$ PB - pcfm(i) = zcfm1(i) - pcfh(i) = zcfh1(i) - ELSE - pcfm(i) = zcdn_m(i)* fsta(zri(i)) - pcfh(i) = zcdn_h(i)* fsta(zri(i)) - ENDIF - ELSE ! situation instable - IF (.NOT.zxli) THEN - zucf = 1./(1.+3.0*CB*CC*zcdn_m(i)*SQRT(ABS(zri(i)) & - *(1.0+zgeop1(i)/(RG*z0m(i))))) - zcfm2(i) = zcdn_m(i)*amax1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min) -!!$ PB zcfh2(i) = zcdn_h(i)*amax1((1.-3.0*cb*zri(i)*zucf),f_ri_cd_min) - zcfh2(i) = f_cdrag_ter*zcdn_h(i)*amax1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min) - pcfm(i) = zcfm2(i) - pcfh(i) = zcfh2(i) - ELSE - pcfm(i) = zcdn_m(i)* fins(zri(i)) - pcfh(i) = zcdn_h(i)* fins(zri(i)) - ENDIF - IF(iflag_gusts==0) THEN -! cdrah sur l'ocean cf. Miller et al. (1992) - only active when gustiness parameterization is not active - zcr = (0.0016/(zcdn_m(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.) - IF(nsrf.EQ.is_oce) pcfh(i) =f_cdrag_oce* zcdn_h(i)*(1.0+zcr**1.25)**(1./1.25) - ENDIF - ENDIF - END DO - + FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min) + FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min ) + + + CASE (3) ! Laurent Li + + FM(i)=MAX(1.0 / (1.0+10.0*zri(i)*(1+8.0*zri(i))),f_ri_cd_min) + FH(i)=FM(i) + + + CASE (4) ! King 2001 + + if (zri(i) .LT. C2/2.) then + FM(i)=MAX((1.-zri(i)/C2)**2,f_ri_cd_min) + FH(i)= FM(i) + + + else + FM(i)=MAX(C3*((C2/zri(i))**2),f_ri_cd_min) + FH(i)= FM(i) + endif + + + CASE (5) ! MO + + if (zri(i) .LT. 1./alpha) then + + FM(i)=MAX((1.-alpha*zri(i))**2,f_ri_cd_min) + FH(i)=FM(i) + + else + + + FM(i)=MAX(1E-7,f_ri_cd_min) + FH(i)=FM(i) + + endif + + + + + + CASE default ! Louis 1982 + + zscf = SQRT(1.+CD*ABS(zri(i))) + FM(i)= AMAX1(1. / (1.+2.*CB*zri(i)/zscf), f_ri_cd_min) + FH(i)= AMAX1(1./ (1.+3.*CB*zri(i)*zscf), f_ri_cd_min ) + + + + END SELECT + +! Calcul des drags + + + cdm(i)=cdmn(i)*FM(i) + cdh(i)=f_cdrag_ter*cdhn(i)*FH(i) + + IF(nsrf.EQ.is_oce) THEN + + cdh(i)=f_cdrag_oce*cdhn(i)*FH(i) + cdm(i)=MIN(cdm(i),cdmmax) + cdh(i)=MIN(cdh(i),cdhmax) + + ENDIF + + + ENDIF + + + + +!xxxxxxxxxxx + END DO ! Fin de la boucle sur l'horizontal +!xxxxxxxxxxx ! ================================================================= c - ! IM cf JLD : on seuille cdrag_m et cdrag_h - IF (nsrf == is_oce) THEN - DO i=1,knon - pcfm(i)=MIN(pcfm(i),cdmmax) - pcfh(i)=MIN(pcfh(i),cdhmax) - END DO - END IF + END SUBROUTINE cdrag + + +!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++