! !$Id: cdrag.F90 2778 2017-01-26 14:49:44Z oboucher $ ! SUBROUTINE cdrag( knon, nsrf, & speed, t1, q1, zgeop1, & psol, tsurf, qsurf, z0m, z0h, & pcfm, pcfh, zri, pref ) USE dimphy USE indice_sol_mod USE print_control_mod, ONLY: lunout, prt_level 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. ! ! References: ! Louis, J. F., 1979: A parametric model of vertical eddy fluxes in the ! atmosphere. Boundary-Layer Meteorology. 01/1979; 17(2):187-202. ! Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short history of the ! operational PBL parametrization at ECMWF'. Workshop on boundary layer ! parametrization, November 1981, ECMWF, Reading, England. ! Page: 19. Equations in Table 1. ! Anton Beljaars. May 1992. The parametrization of the planetary boundary layer. ! European Centre for Medium-Range Weather Forecasts. ! Equations: 110-113. Page 40. ! Miller,M.J., A.C.M.Beljaars, T.N.Palmer. 1992. The sensitivity of the ECMWF ! model to the parameterization of evaporation from the tropical oceans. J. ! Climate, 5:418-434. ! ! ================================================================= 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 ! ! ! Parametres d'entree !***************************************************************** INTEGER, INTENT(IN) :: knon, nsrf 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 ! 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) :: 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 ! 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 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) ! ================================================================= 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 ! ! Fuxing WANG, 04/03/2015 ! To check if there are negative q1, qsurf values. ng_q1 = 0 ! Initialization ng_qsurf = 0 ! Initialization DO i = 1, knon IF (q1(i).LT.0.0) ng_q1 = ng_q1 + 1 IF (qsurf(i).LT.0.0) ng_qsurf = ng_qsurf + 1 ENDDO IF (ng_q1.GT.0 .and. prt_level > 5) THEN WRITE(lunout,*)" *** Warning: Negative q1(humidity at 1st level) values in cdrag.F90 !" WRITE(lunout,*)" The total number of the grids is: ", ng_q1 WRITE(lunout,*)" The negative q1 is set to zero " ! abort_message="voir ci-dessus" ! CALL abort_physic(modname,abort_message,1) ENDIF IF (ng_qsurf.GT.0 .and. prt_level > 5) THEN WRITE(lunout,*)" *** Warning: Negative qsurf(humidity at surface) values in cdrag.F90 !" WRITE(lunout,*)" The total number of the grids is: ", ng_qsurf WRITE(lunout,*)" The negative qsurf is set to zero " ! abort_message="voir ci-dessus" ! CALL abort_physic(modname,abort_message,1) 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) !------------------------------------------------------------ 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 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)) 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 zri(i) = MIN(20.,zri(i)) IF (.NOT.zxli) THEN 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 ! ================================================================= 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