| [2232] | 1 | ! |
|---|
| 2 | !$Id: cdrag.F90 ???? 2015-??-?? ??:??:??Z ? $ |
|---|
| 3 | ! |
|---|
| 4 | SUBROUTINE cdrag( knon, nsrf, & |
|---|
| 5 | speed, t1, q1, zgeop1, & |
|---|
| [2243] | 6 | psol, tsurf, qsurf, z0m, z0h, & |
|---|
| [2232] | 7 | pcfm, pcfh, zri, pref ) |
|---|
| 8 | |
|---|
| 9 | USE dimphy |
|---|
| 10 | USE indice_sol_mod |
|---|
| [2311] | 11 | USE print_control_mod, ONLY: lunout |
|---|
| [2232] | 12 | IMPLICIT NONE |
|---|
| 13 | ! ================================================================= c |
|---|
| 14 | ! |
|---|
| 15 | ! Objet : calcul des cdrags pour le moment (pcfm) et |
|---|
| 16 | ! les flux de chaleur sensible et latente (pcfh). |
|---|
| 17 | ! |
|---|
| 18 | ! Modified histroy: |
|---|
| 19 | ! 27-Jan-2014: richardson number inconsistant between |
|---|
| 20 | ! coefcdrag.F90 and clcdrag.F90, Fuxing WANG wrote this subroutine |
|---|
| 21 | ! by merging coefcdrag and clcdrag. |
|---|
| 22 | ! |
|---|
| 23 | ! References: |
|---|
| 24 | ! Louis, J. F., 1979: A parametric model of vertical eddy fluxes in the |
|---|
| 25 | ! atmosphere. Boundary-Layer Meteorology. 01/1979; 17(2):187-202. |
|---|
| 26 | ! Louis, J. F., Tiedtke, M. and Geleyn, J. F., 1982: `A short history of the |
|---|
| 27 | ! operational PBL parametrization at ECMWF'. Workshop on boundary layer |
|---|
| 28 | ! parametrization, November 1981, ECMWF, Reading, England. |
|---|
| 29 | ! Page: 19. Equations in Table 1. |
|---|
| 30 | ! Anton Beljaars. May 1992. The parametrization of the planetary boundary layer. |
|---|
| 31 | ! European Centre for Medium-Range Weather Forecasts. |
|---|
| 32 | ! Equations: 110-113. Page 40. |
|---|
| 33 | ! Miller,M.J., A.C.M.Beljaars, T.N.Palmer. 1992. The sensitivity of the ECMWF |
|---|
| 34 | ! model to the parameterization of evaporation from the tropical oceans. J. |
|---|
| 35 | ! Climate, 5:418-434. |
|---|
| 36 | ! |
|---|
| 37 | ! ================================================================= c |
|---|
| 38 | ! |
|---|
| 39 | ! knon----input-I- nombre de points pour un type de surface |
|---|
| 40 | ! nsrf----input-I- indice pour le type de surface; voir indicesol.h |
|---|
| 41 | ! speed---input-R- module du vent au 1er niveau du modele |
|---|
| 42 | ! t1------input-R- temperature de l'air au 1er niveau du modele |
|---|
| 43 | ! q1------input-R- humidite de l'air au 1er niveau du modele |
|---|
| 44 | ! zgeop---input-R- geopotentiel au 1er niveau du modele |
|---|
| 45 | ! psol----input-R- pression au sol |
|---|
| 46 | ! tsurf---input-R- temperature de l'air a la surface |
|---|
| 47 | ! qsurf---input-R- humidite de l'air a la surface |
|---|
| [2243] | 48 | ! z0m, z0h---input-R- rugosite |
|---|
| [2232] | 49 | !! u1, v1 are removed, speed is used. Fuxing WANG, 04/03/2015, |
|---|
| 50 | !! u1------input-R- vent zonal au 1er niveau du modele |
|---|
| 51 | !! v1------input-R- vent meridien au 1er niveau du modele |
|---|
| 52 | ! |
|---|
| 53 | ! pcfm---output-R- cdrag pour le moment |
|---|
| 54 | ! pcfh---output-R- cdrag pour les flux de chaleur latente et sensible |
|---|
| 55 | ! zri----output-R- Richardson number |
|---|
| 56 | ! pref---output-R- pression au niveau zgeop/RG |
|---|
| 57 | ! |
|---|
| 58 | ! Parameters: |
|---|
| 59 | ! ckap-----Karman constant |
|---|
| 60 | ! cb,cc,cd-parameters in Louis et al., 1982 |
|---|
| 61 | ! ================================================================= c |
|---|
| 62 | ! |
|---|
| 63 | ! |
|---|
| 64 | ! Parametres d'entree |
|---|
| 65 | !***************************************************************** |
|---|
| 66 | INTEGER, INTENT(IN) :: knon, nsrf |
|---|
| 67 | REAL, DIMENSION(klon), INTENT(IN) :: speed ! module du vent au 1er niveau du modele |
|---|
| 68 | REAL, DIMENSION(klon), INTENT(IN) :: zgeop1! geopotentiel au 1er niveau du modele |
|---|
| 69 | REAL, DIMENSION(klon), INTENT(IN) :: psol ! pression au sol |
|---|
| 70 | REAL, DIMENSION(klon), INTENT(IN) :: t1 ! temperature at 1st level |
|---|
| 71 | REAL, DIMENSION(klon), INTENT(IN) :: q1 ! humidity at 1st level |
|---|
| 72 | REAL, DIMENSION(klon), INTENT(IN) :: tsurf ! Surface temperature (K) |
|---|
| 73 | REAL, DIMENSION(klon), INTENT(IN) :: qsurf ! Surface humidity (Kg/Kg) |
|---|
| [2243] | 74 | REAL, DIMENSION(klon), INTENT(IN) :: z0m, z0h ! Rugosity at surface (m) |
|---|
| [2232] | 75 | ! paprs, pplay u1, v1: to be deleted |
|---|
| 76 | ! they were in the old clcdrag. Fuxing WANG, 04/03/2015 |
|---|
| 77 | ! REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs |
|---|
| 78 | ! REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay |
|---|
| 79 | ! REAL, DIMENSION(klon), INTENT(IN) :: u1, v1 |
|---|
| 80 | |
|---|
| 81 | ! Parametres de sortie |
|---|
| 82 | !****************************************************************** |
|---|
| 83 | REAL, DIMENSION(klon), INTENT(OUT) :: pcfm ! Drag coefficient for heat flux |
|---|
| 84 | REAL, DIMENSION(klon), INTENT(OUT) :: pcfh ! Drag coefficient for momentum |
|---|
| 85 | REAL, DIMENSION(klon), INTENT(OUT) :: zri ! Richardson number |
|---|
| 86 | REAL, DIMENSION(klon), INTENT(OUT) :: pref ! Pression au niveau zgeop/RG |
|---|
| 87 | |
|---|
| 88 | ! Parametres local |
|---|
| 89 | INTEGER :: ng_q1 ! Number of grids that q1 < 0.0 |
|---|
| 90 | INTEGER :: ng_qsurf ! Number of grids that qsurf < 0.0 |
|---|
| 91 | ! zgeop1, psol: to be deleted, they are inputs now. Fuxing WANG, 04/03/2015 |
|---|
| 92 | ! REAL, DIMENSION(klon) :: zgeop1! geopotentiel au 1er niveau du modele |
|---|
| 93 | ! REAL, DIMENSION(klon) :: psol ! pression au sol |
|---|
| 94 | ! |
|---|
| 95 | ! ================================================================= c |
|---|
| 96 | ! |
|---|
| 97 | INCLUDE "YOMCST.h" |
|---|
| 98 | INCLUDE "YOETHF.h" |
|---|
| 99 | ! INCLUDE "indicesol.h" |
|---|
| 100 | INCLUDE "clesphys.h" |
|---|
| 101 | ! |
|---|
| 102 | ! Quelques constantes et options: |
|---|
| 103 | !!$PB REAL, PARAMETER :: ckap=0.35, cb=5.0, cc=5.0, cd=5.0, cepdu2=(0.1)**2 |
|---|
| 104 | REAL, PARAMETER :: CKAP=0.40, CB=5.0, CC=5.0, CD=5.0, CEPDU2 = (0.1)**2 |
|---|
| 105 | ! |
|---|
| 106 | ! Variables locales : |
|---|
| 107 | INTEGER :: i |
|---|
| 108 | REAL :: zdu2, ztsolv |
|---|
| 109 | REAL :: ztvd, zscf |
|---|
| 110 | REAL :: zucf, zcr |
|---|
| 111 | REAL :: friv, frih |
|---|
| 112 | REAL, DIMENSION(klon) :: zcfm1, zcfm2 ! Drag coefficient for momentum |
|---|
| 113 | REAL, DIMENSION(klon) :: zcfh1, zcfh2 ! Drag coefficient for heat flux |
|---|
| 114 | LOGICAL, PARAMETER :: zxli=.FALSE. ! calcul des cdrags selon Laurent Li |
|---|
| [2243] | 115 | REAL, DIMENSION(klon) :: zcdn_m, zcdn_h ! Drag coefficient in neutral conditions |
|---|
| [2441] | 116 | REAL zzzcd |
|---|
| [2232] | 117 | ! |
|---|
| 118 | ! Fonctions thermodynamiques et fonctions d'instabilite |
|---|
| 119 | REAL :: fsta, fins, x |
|---|
| 120 | fsta(x) = 1.0 / (1.0+10.0*x*(1+8.0*x)) |
|---|
| 121 | fins(x) = SQRT(1.0-18.0*x) |
|---|
| 122 | |
|---|
| 123 | ! ================================================================= c |
|---|
| 124 | ! Fuxing WANG, 04/03/2015, delete the calculation of zgeop1 |
|---|
| 125 | ! (le geopotentiel du premier couche de modele). |
|---|
| 126 | ! zgeop1 is an input ivariable in this subroutine. |
|---|
| 127 | ! DO i = 1, knon |
|---|
| 128 | ! zgeop1(i) = RD * t1(i) / (0.5*(paprs(i,1)+pplay(i,1))) & |
|---|
| 129 | ! * (paprs(i,1)-pplay(i,1)) |
|---|
| 130 | ! END DO |
|---|
| 131 | ! ================================================================= c |
|---|
| 132 | ! |
|---|
| 133 | ! Fuxing WANG, 04/03/2015 |
|---|
| 134 | ! To check if there are negative q1, qsurf values. |
|---|
| 135 | ng_q1 = 0 ! Initialization |
|---|
| 136 | ng_qsurf = 0 ! Initialization |
|---|
| 137 | DO i = 1, knon |
|---|
| 138 | IF (q1(i).LT.0.0) ng_q1 = ng_q1 + 1 |
|---|
| 139 | IF (qsurf(i).LT.0.0) ng_qsurf = ng_qsurf + 1 |
|---|
| 140 | ENDDO |
|---|
| 141 | IF (ng_q1.GT.0) THEN |
|---|
| 142 | WRITE(lunout,*)" *** Warning: Negative q1(humidity at 1st level) values in cdrag.F90 !" |
|---|
| 143 | WRITE(lunout,*)" The total number of the grids is: ", ng_q1 |
|---|
| 144 | WRITE(lunout,*)" The negative q1 is set to zero " |
|---|
| 145 | ! abort_message="voir ci-dessus" |
|---|
| [2311] | 146 | ! CALL abort_physic(modname,abort_message,1) |
|---|
| [2232] | 147 | ENDIF |
|---|
| 148 | IF (ng_qsurf.GT.0) THEN |
|---|
| 149 | WRITE(lunout,*)" *** Warning: Negative qsurf(humidity at surface) values in cdrag.F90 !" |
|---|
| 150 | WRITE(lunout,*)" The total number of the grids is: ", ng_qsurf |
|---|
| 151 | WRITE(lunout,*)" The negative qsurf is set to zero " |
|---|
| 152 | ! abort_message="voir ci-dessus" |
|---|
| [2311] | 153 | ! CALL abort_physic(modname,abort_message,1) |
|---|
| [2232] | 154 | ENDIF |
|---|
| 155 | |
|---|
| 156 | ! Calculer le frottement au sol (Cdrag) |
|---|
| 157 | DO i = 1, knon |
|---|
| 158 | !------------------------------------------------------------ |
|---|
| 159 | ! u1, v1 are replaced by speed. Fuxing WANG, 04/03/2015, |
|---|
| 160 | ! zdu2 = MAX(CEPDU2,u1(i)**2+v1(i)**2) |
|---|
| 161 | !------------------------------------------------------------ |
|---|
| 162 | zdu2 = MAX(CEPDU2, speed(i)**2) |
|---|
| 163 | ! psol(i) = paprs(i,1) |
|---|
| 164 | pref(i) = EXP(LOG(psol(i)) - zgeop1(i)/(RD*t1(i)* & |
|---|
| 165 | (1.+ RETV * max(q1(i),0.0)))) ! negative q1 set to zero |
|---|
| 166 | !------------ the old calculations in clcdrag---------------- |
|---|
| 167 | ! ztsolv = tsurf(i) * (1.0+RETV*qsurf(i)) |
|---|
| 168 | ! ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) & |
|---|
| 169 | ! *(1.+RETV*q1(i)) |
|---|
| 170 | !------------------------------------------------------------ |
|---|
| 171 | ! Fuxing WANG, 04/03/2015, in this revised version, |
|---|
| 172 | ! the negative qsurf and q1 are set to zero (as in coefcdrag) |
|---|
| 173 | ztsolv = tsurf(i) * (1.0+RETV*max(qsurf(i),0.0)) ! negative qsurf set to zero |
|---|
| 174 | ztvd = (t1(i)+zgeop1(i)/RCPD/(1.+RVTMP2*q1(i))) & |
|---|
| 175 | *(1.+RETV*max(q1(i),0.0)) ! negative q1 set to zero |
|---|
| 176 | zri(i) = zgeop1(i)*(ztvd-ztsolv)/(zdu2*ztvd) |
|---|
| 177 | |
|---|
| [2243] | 178 | |
|---|
| [2441] | 179 | ! Coefficients CD neutres pour m et h : k^2/ln(z/z0) et k^2/(ln(z/z0)*ln(z/z0h)) |
|---|
| 180 | zzzcd=CKAP/LOG(1.+zgeop1(i)/(RG*z0m(i))) |
|---|
| 181 | zcdn_m(i) = zzzcd*zzzcd |
|---|
| [2460] | 182 | zcdn_h(i) = zzzcd*(CKAP/LOG(1.+zgeop1(i)/(RG*z0h(i)))) |
|---|
| [2243] | 183 | |
|---|
| [2232] | 184 | IF (zri(i) .GT. 0.) THEN ! situation stable |
|---|
| 185 | zri(i) = MIN(20.,zri(i)) |
|---|
| 186 | IF (.NOT.zxli) THEN |
|---|
| 187 | zscf = SQRT(1.+CD*ABS(zri(i))) |
|---|
| 188 | friv = AMAX1(1. / (1.+2.*CB*zri(i)/ZSCF), f_ri_cd_min) |
|---|
| [2243] | 189 | zcfm1(i) = zcdn_m(i) * friv |
|---|
| [2232] | 190 | frih = AMAX1(1./ (1.+3.*CB*zri(i)*ZSCF), f_ri_cd_min ) |
|---|
| 191 | !!$ PB zcfh1(i) = zcdn(i) * frih |
|---|
| 192 | !!$ PB zcfh1(i) = f_cdrag_stable * zcdn(i) * frih |
|---|
| [2243] | 193 | zcfh1(i) = f_cdrag_ter * zcdn_h(i) * frih |
|---|
| 194 | IF(nsrf.EQ.is_oce) zcfh1(i) = f_cdrag_oce * zcdn_h(i) * frih |
|---|
| [2232] | 195 | !!$ PB |
|---|
| 196 | pcfm(i) = zcfm1(i) |
|---|
| 197 | pcfh(i) = zcfh1(i) |
|---|
| 198 | ELSE |
|---|
| [2243] | 199 | pcfm(i) = zcdn_m(i)* fsta(zri(i)) |
|---|
| 200 | pcfh(i) = zcdn_h(i)* fsta(zri(i)) |
|---|
| [2232] | 201 | ENDIF |
|---|
| 202 | ELSE ! situation instable |
|---|
| 203 | IF (.NOT.zxli) THEN |
|---|
| [2243] | 204 | zucf = 1./(1.+3.0*CB*CC*zcdn_m(i)*SQRT(ABS(zri(i)) & |
|---|
| 205 | *(1.0+zgeop1(i)/(RG*z0m(i))))) |
|---|
| 206 | zcfm2(i) = zcdn_m(i)*amax1((1.-2.0*CB*zri(i)*zucf),f_ri_cd_min) |
|---|
| 207 | !!$ PB zcfh2(i) = zcdn_h(i)*amax1((1.-3.0*cb*zri(i)*zucf),f_ri_cd_min) |
|---|
| 208 | zcfh2(i) = f_cdrag_ter*zcdn_h(i)*amax1((1.-3.0*CB*zri(i)*zucf),f_ri_cd_min) |
|---|
| [2232] | 209 | pcfm(i) = zcfm2(i) |
|---|
| 210 | pcfh(i) = zcfh2(i) |
|---|
| 211 | ELSE |
|---|
| [2243] | 212 | pcfm(i) = zcdn_m(i)* fins(zri(i)) |
|---|
| 213 | pcfh(i) = zcdn_h(i)* fins(zri(i)) |
|---|
| [2232] | 214 | ENDIF |
|---|
| [2278] | 215 | IF(iflag_gusts==0) THEN |
|---|
| 216 | ! cdrah sur l'ocean cf. Miller et al. (1992) - only active when gustiness parameterization is not active |
|---|
| 217 | zcr = (0.0016/(zcdn_m(i)*SQRT(zdu2)))*ABS(ztvd-ztsolv)**(1./3.) |
|---|
| 218 | IF(nsrf.EQ.is_oce) pcfh(i) =f_cdrag_oce* zcdn_h(i)*(1.0+zcr**1.25)**(1./1.25) |
|---|
| 219 | ENDIF |
|---|
| [2232] | 220 | ENDIF |
|---|
| 221 | END DO |
|---|
| 222 | |
|---|
| 223 | ! ================================================================= c |
|---|
| 224 | |
|---|
| 225 | ! IM cf JLD : on seuille cdrag_m et cdrag_h |
|---|
| 226 | IF (nsrf == is_oce) THEN |
|---|
| 227 | DO i=1,knon |
|---|
| 228 | pcfm(i)=MIN(pcfm(i),cdmmax) |
|---|
| 229 | pcfh(i)=MIN(pcfh(i),cdhmax) |
|---|
| 230 | END DO |
|---|
| 231 | END IF |
|---|
| 232 | |
|---|
| 233 | END SUBROUTINE cdrag |
|---|
