[787] | 1 | subroutine moistadj(ngrid, nq, pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb) |
---|
[135] | 2 | |
---|
[728] | 3 | use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD, Psat_water, Lcpdqsat_water |
---|
[787] | 4 | USE tracer_h |
---|
[135] | 5 | |
---|
| 6 | implicit none |
---|
| 7 | |
---|
| 8 | |
---|
| 9 | !===================================================================== |
---|
| 10 | ! |
---|
| 11 | ! Purpose |
---|
| 12 | ! ------- |
---|
| 13 | ! Calculates moist convective adjustment by the method of Manabe. |
---|
| 14 | ! |
---|
| 15 | ! Authors |
---|
| 16 | ! ------- |
---|
| 17 | ! Adapted from the LMDTERRE code by R. Wordsworth (2010) |
---|
| 18 | ! Original author Z. X. Li (1993) |
---|
| 19 | ! |
---|
| 20 | !===================================================================== |
---|
| 21 | |
---|
| 22 | #include "dimensions.h" |
---|
| 23 | #include "dimphys.h" |
---|
| 24 | #include "comcstfi.h" |
---|
| 25 | |
---|
[787] | 26 | INTEGER ngrid, nq |
---|
[135] | 27 | |
---|
[787] | 28 | REAL pt(ngrid,nlayermx) ! temperature (K) |
---|
| 29 | REAL pq(ngrid,nlayermx,nq) ! tracer (kg/kg) |
---|
| 30 | REAL pdq(ngrid,nlayermx,nq) |
---|
[135] | 31 | |
---|
[787] | 32 | REAL pdqmana(ngrid,nlayermx,nq) ! tendency of tracers (kg/kg.s-1) |
---|
| 33 | REAL pdtmana(ngrid,nlayermx) ! temperature increment |
---|
| 34 | |
---|
[728] | 35 | ! local variables |
---|
[787] | 36 | REAL zt(ngrid,nlayermx) ! temperature (K) |
---|
| 37 | REAL zq(ngrid,nlayermx) ! humidite specifique (kg/kg) |
---|
| 38 | REAL pplev(ngrid,nlayermx+1) ! pression a inter-couche (Pa) |
---|
| 39 | REAL pplay(ngrid,nlayermx) ! pression au milieu de couche (Pa) |
---|
[135] | 40 | |
---|
[787] | 41 | REAL d_t(ngrid,nlayermx) ! temperature increment |
---|
| 42 | REAL d_q(ngrid,nlayermx) ! incrementation pour vapeur d'eau |
---|
| 43 | REAL d_ql(ngrid,nlayermx) ! incrementation pour l'eau liquide |
---|
| 44 | REAL rneb(ngrid,nlayermx) ! cloud fraction |
---|
[135] | 45 | REAL ptimestep |
---|
| 46 | |
---|
| 47 | ! REAL t_coup |
---|
| 48 | ! PARAMETER (t_coup=234.0) |
---|
| 49 | REAL seuil_vap |
---|
| 50 | PARAMETER (seuil_vap=1.0E-10) |
---|
| 51 | |
---|
| 52 | ! Local variables |
---|
[728] | 53 | INTEGER i, k, iq, nn |
---|
[773] | 54 | INTEGER, PARAMETER :: niter = 1 |
---|
[135] | 55 | INTEGER k1, k1p, k2, k2p |
---|
[787] | 56 | LOGICAL itest(ngrid) |
---|
| 57 | REAL delta_q(ngrid, nlayermx) |
---|
[135] | 58 | REAL cp_new_t(nlayermx) |
---|
| 59 | REAL cp_delta_t(nlayermx) |
---|
| 60 | REAL v_cptj(nlayermx), v_cptjk1, v_ssig |
---|
[787] | 61 | REAL v_cptt(ngrid,nlayermx), v_p, v_t, zpsat |
---|
| 62 | REAL zqs(ngrid,nlayermx), zdqs(ngrid,nlayermx) |
---|
| 63 | REAL zq1(ngrid), zq2(ngrid) |
---|
| 64 | REAL gamcpdz(ngrid,2:nlayermx) |
---|
[135] | 65 | REAL zdp, zdpm |
---|
| 66 | |
---|
| 67 | REAL zsat ! super-saturation |
---|
| 68 | REAL zflo ! flotabilite |
---|
| 69 | |
---|
[787] | 70 | REAL local_q(ngrid,nlayermx),local_t(ngrid,nlayermx) |
---|
[135] | 71 | |
---|
| 72 | REAL zdelta, zcor, zcvm5 |
---|
| 73 | |
---|
| 74 | REAL dEtot, dqtot, masse ! conservation diagnostics |
---|
| 75 | real dL1tot, dL2tot |
---|
| 76 | |
---|
| 77 | ! Indices of water vapour and water ice tracers |
---|
| 78 | INTEGER,SAVE :: i_h2o=0 ! water vapour |
---|
| 79 | INTEGER,SAVE :: i_ice=0 ! water ice |
---|
| 80 | |
---|
| 81 | LOGICAL firstcall |
---|
| 82 | SAVE firstcall |
---|
| 83 | |
---|
| 84 | DATA firstcall /.TRUE./ |
---|
| 85 | |
---|
| 86 | IF (firstcall) THEN |
---|
| 87 | |
---|
| 88 | i_h2o=igcm_h2o_vap |
---|
| 89 | i_ice=igcm_h2o_ice |
---|
| 90 | |
---|
| 91 | write(*,*) "rain: i_ice=",i_ice |
---|
| 92 | write(*,*) " i_h2o=",i_h2o |
---|
| 93 | |
---|
| 94 | firstcall = .FALSE. |
---|
| 95 | ENDIF |
---|
| 96 | |
---|
| 97 | ! GCM -----> subroutine variables |
---|
[787] | 98 | zq(1:ngrid,1:nlayermx) = pq(1:ngrid,1:nlayermx,i_h2o)+ pdq(1:ngrid,1:nlayermx,i_h2o)*ptimestep |
---|
| 99 | zt(1:ngrid,1:nlayermx) = pt(1:ngrid,1:nlayermx) |
---|
| 100 | pdqmana(1:ngrid,1:nlayermx,1:nq)=0.0 |
---|
[728] | 101 | |
---|
[135] | 102 | DO k = 1, nlayermx |
---|
[787] | 103 | DO i = 1, ngrid |
---|
[728] | 104 | if(zq(i,k).lt.0.)then |
---|
| 105 | zq(i,k)=0.0 |
---|
[135] | 106 | endif |
---|
[728] | 107 | ENDDO |
---|
[135] | 108 | ENDDO |
---|
[728] | 109 | |
---|
[787] | 110 | local_q(1:ngrid,1:nlayermx) = zq(1:ngrid,1:nlayermx) |
---|
| 111 | local_t(1:ngrid,1:nlayermx) = zt(1:ngrid,1:nlayermx) |
---|
| 112 | rneb(1:ngrid,1:nlayermx) = 0.0 |
---|
| 113 | d_ql(1:ngrid,1:nlayermx) = 0.0 |
---|
| 114 | d_t(1:ngrid,1:nlayermx) = 0.0 |
---|
| 115 | d_q(1:ngrid,1:nlayermx) = 0.0 |
---|
[135] | 116 | |
---|
| 117 | ! Calculate v_cptt |
---|
| 118 | DO k = 1, nlayermx |
---|
[787] | 119 | DO i = 1, ngrid |
---|
[135] | 120 | v_cptt(i,k) = RCPD * local_t(i,k) |
---|
[728] | 121 | v_t = MAX(local_t(i,k),15.) |
---|
[135] | 122 | v_p = pplay(i,k) |
---|
| 123 | |
---|
[728] | 124 | call Psat_water(v_t,v_p,zpsat,zqs(i,k)) |
---|
| 125 | call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k)) |
---|
[135] | 126 | ENDDO |
---|
| 127 | ENDDO |
---|
| 128 | |
---|
| 129 | ! Calculate Gamma * Cp * dz: (gamma is the critical gradient) |
---|
| 130 | DO k = 2, nlayermx |
---|
[787] | 131 | DO i = 1, ngrid |
---|
[135] | 132 | zdp = pplev(i,k)-pplev(i,k+1) |
---|
| 133 | zdpm = pplev(i,k-1)-pplev(i,k) |
---|
[728] | 134 | gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp) & |
---|
| 135 | + RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) ) & |
---|
| 136 | * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) & |
---|
| 137 | / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) ) |
---|
[135] | 138 | ENDDO |
---|
| 139 | ENDDO |
---|
| 140 | |
---|
| 141 | !------------------------------------ modification of unstable profile |
---|
[787] | 142 | DO 9999 i = 1, ngrid |
---|
| 143 | |
---|
[135] | 144 | itest(i) = .FALSE. |
---|
| 145 | |
---|
| 146 | ! print*,'we in the loop' |
---|
| 147 | ! stop |
---|
| 148 | |
---|
| 149 | k1 = 0 |
---|
| 150 | k2 = 1 |
---|
| 151 | |
---|
| 152 | 810 CONTINUE ! look for k1, the base of the column |
---|
| 153 | k2 = k2 + 1 |
---|
| 154 | IF (k2 .GT. nlayermx) GOTO 9999 |
---|
| 155 | zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2) |
---|
[728] | 156 | zsat=(local_q(i,k2-1)-zqs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2)) & |
---|
| 157 | +(local_q(i,k2)-zqs(i,k2))*(pplev(i,k2)-pplev(i,k2+1)) |
---|
[135] | 158 | |
---|
| 159 | IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 ) GOTO 810 |
---|
| 160 | k1 = k2 - 1 |
---|
| 161 | itest(i) = .TRUE. |
---|
| 162 | |
---|
| 163 | 820 CONTINUE !! look for k2, the top of the column |
---|
| 164 | IF (k2 .EQ. nlayermx) GOTO 821 |
---|
| 165 | k2p = k2 + 1 |
---|
[728] | 166 | zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-zqs(i,k2p)) |
---|
[135] | 167 | zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p) |
---|
| 168 | |
---|
| 169 | IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 821 |
---|
| 170 | k2 = k2p |
---|
| 171 | GOTO 820 |
---|
| 172 | 821 CONTINUE |
---|
| 173 | |
---|
| 174 | !------------------------------------------------------ local adjustment |
---|
| 175 | 830 CONTINUE ! actual adjustment |
---|
[728] | 176 | Do nn=1,niter |
---|
[135] | 177 | v_cptj(k1) = 0.0 |
---|
| 178 | zdp = pplev(i,k1)-pplev(i,k1+1) |
---|
[728] | 179 | v_cptjk1 = ( (1.0+zdqs(i,k1))*(v_cptt(i,k1)+v_cptj(k1)) + RLVTT*(local_q(i,k1)-zqs(i,k1)) ) * zdp |
---|
| 180 | v_ssig = zdp * (1.0+zdqs(i,k1)) |
---|
[135] | 181 | |
---|
| 182 | k1p = k1 + 1 |
---|
| 183 | DO k = k1p, k2 |
---|
| 184 | zdp = pplev(i,k)-pplev(i,k+1) |
---|
| 185 | v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k) |
---|
[728] | 186 | v_cptjk1 = v_cptjk1 + zdp * ( (1.0+zdqs(i, k))*(v_cptt(i,k)+v_cptj(k)) + RLVTT*(local_q(i,k)-zqs(i,k)) ) |
---|
| 187 | v_ssig = v_ssig + zdp *(1.0+zdqs(i,k)) |
---|
[135] | 188 | ENDDO |
---|
| 189 | |
---|
| 190 | |
---|
| 191 | ! this right here is where the adjustment is done??? |
---|
| 192 | DO k = k1, k2 |
---|
| 193 | cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k) |
---|
| 194 | cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k) |
---|
[728] | 195 | v_cptt(i,k)=cp_new_t(k) |
---|
| 196 | local_q(i,k) = zqs(i,k) + zdqs(i,k)*cp_delta_t(k)/RLVTT |
---|
[135] | 197 | local_t(i,k) = cp_new_t(k) / RCPD |
---|
[253] | 198 | |
---|
[728] | 199 | v_t = MAX(local_t(i,k),15.) |
---|
| 200 | v_p = pplay(i,k) |
---|
| 201 | |
---|
| 202 | call Psat_water(v_t,v_p,zpsat,zqs(i,k)) |
---|
| 203 | call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k)) |
---|
| 204 | |
---|
| 205 | |
---|
| 206 | |
---|
| 207 | ! print*,'i,k,zqs,cpdT=',i,k,zqs(i,k),cp_delta_t(k) |
---|
[135] | 208 | ENDDO |
---|
[728] | 209 | Enddo ! nn=1,niter |
---|
[135] | 210 | |
---|
[253] | 211 | |
---|
[135] | 212 | !--------------------------------------------------- sounding downwards |
---|
| 213 | ! -- we refine the prognostic variables in |
---|
| 214 | ! -- the layer about to be adjusted |
---|
| 215 | |
---|
[728] | 216 | ! DO k = k1, k2 |
---|
| 217 | ! v_cptt(i,k) = RCPD * local_t(i,k) |
---|
| 218 | ! v_t = local_t(i,k) |
---|
| 219 | ! v_p = pplay(i,k) |
---|
| 220 | ! |
---|
| 221 | ! call Psat_water(v_t,v_p,zpsat,zqs(i,k)) |
---|
| 222 | ! call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k)) |
---|
| 223 | ! ENDDO |
---|
[135] | 224 | |
---|
| 225 | DO k = 2, nlayermx |
---|
| 226 | zdpm = pplev(i,k-1) - pplev(i,k) |
---|
| 227 | zdp = pplev(i,k) - pplev(i,k+1) |
---|
[728] | 228 | gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp) & |
---|
| 229 | + RLVTT /(zdpm+zdp) * (zqs(i,k-1)*zdpm + zqs(i,k)*zdp) ) & |
---|
| 230 | * (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) & |
---|
| 231 | / (1.0+ (zdqs(i,k-1)*zdpm + zdqs(i,k)*zdp)/(zdpm+zdp) ) |
---|
[135] | 232 | ENDDO |
---|
| 233 | |
---|
| 234 | ! Test to see if we've reached the bottom |
---|
| 235 | |
---|
| 236 | IF (k1 .EQ. 1) GOTO 841 ! yes we have! |
---|
| 237 | zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1) |
---|
[728] | 238 | zsat=(local_q(i,k1-1)-zqs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1)) & |
---|
| 239 | + (local_q(i,k1)-zqs(i,k1))*(pplev(i,k1)-pplev(i,k1+1)) |
---|
[135] | 240 | IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 841 ! yes we have! |
---|
| 241 | |
---|
| 242 | 840 CONTINUE |
---|
| 243 | k1 = k1 - 1 |
---|
| 244 | IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995) |
---|
[728] | 245 | zsat = zsat + (local_q(i,k1-1)-zqs(i,k1-1)) & |
---|
[135] | 246 | *(pplev(i,k1-1)-pplev(i,k1)) |
---|
| 247 | zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1) |
---|
| 248 | IF (zflo.GT.0.0 .AND. zsat.GT.0.0) THEN |
---|
| 249 | GOTO 840 |
---|
| 250 | ELSE |
---|
| 251 | GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995) |
---|
| 252 | ENDIF |
---|
| 253 | 841 CONTINUE |
---|
| 254 | |
---|
| 255 | GOTO 810 ! look for other layers higher up |
---|
| 256 | |
---|
| 257 | 9999 CONTINUE ! loop over all the points |
---|
| 258 | |
---|
| 259 | ! print*,'k1=',k1 |
---|
| 260 | ! print*,'k2=',k2 |
---|
| 261 | |
---|
| 262 | ! print*,'local_t=',local_t |
---|
| 263 | ! print*,'v_cptt=',v_cptt |
---|
| 264 | ! print*,'gamcpdz=',gamcpdz |
---|
| 265 | |
---|
| 266 | !----------------------------------------------------------------------- |
---|
| 267 | ! Determine the cloud fraction (hypothese: la nebulosite a lieu |
---|
| 268 | ! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement): |
---|
| 269 | |
---|
| 270 | DO k = 1, nlayermx |
---|
[787] | 271 | DO i = 1, ngrid |
---|
[135] | 272 | IF (itest(i)) THEN |
---|
[728] | 273 | delta_q(i,k) = local_q(i,k) - zq(i,k) |
---|
[135] | 274 | IF (delta_q(i,k).LT.0.) rneb(i,k) = 1.0 |
---|
| 275 | ENDIF |
---|
| 276 | ENDDO |
---|
| 277 | ENDDO |
---|
| 278 | |
---|
| 279 | ! Distribuer l'eau condensee en eau liquide nuageuse (hypothese: |
---|
| 280 | ! l'eau liquide est distribuee aux endroits ou la vapeur d'eau |
---|
| 281 | ! diminue et d'une maniere proportionnelle a cet diminution): |
---|
| 282 | |
---|
[787] | 283 | DO i = 1, ngrid |
---|
[135] | 284 | IF (itest(i)) THEN |
---|
| 285 | zq1(i) = 0.0 |
---|
| 286 | zq2(i) = 0.0 |
---|
| 287 | ENDIF |
---|
| 288 | ENDDO |
---|
| 289 | DO k = 1, nlayermx |
---|
[787] | 290 | DO i = 1, ngrid |
---|
[135] | 291 | IF (itest(i)) THEN |
---|
| 292 | zdp = pplev(i,k)-pplev(i,k+1) |
---|
| 293 | zq1(i) = zq1(i) - delta_q(i,k) * zdp |
---|
| 294 | zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp |
---|
| 295 | ENDIF |
---|
| 296 | ENDDO |
---|
| 297 | ENDDO |
---|
| 298 | DO k = 1, nlayermx |
---|
[787] | 299 | DO i = 1, ngrid |
---|
[135] | 300 | IF (itest(i)) THEN |
---|
[728] | 301 | IF (zq2(i).NE.0.0) d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i) |
---|
[135] | 302 | ENDIF |
---|
| 303 | ENDDO |
---|
| 304 | ENDDO |
---|
| 305 | |
---|
[253] | 306 | ! print*,'local_q BEFORE=',local_q |
---|
| 307 | |
---|
[135] | 308 | DO k = 1, nlayermx |
---|
[787] | 309 | DO i = 1, ngrid |
---|
[135] | 310 | local_q(i, k) = MAX(local_q(i, k), seuil_vap) |
---|
| 311 | ENDDO |
---|
| 312 | ENDDO |
---|
| 313 | |
---|
| 314 | DO k = 1, nlayermx |
---|
[787] | 315 | DO i = 1, ngrid |
---|
[728] | 316 | d_t(i,k) = local_t(i,k) - zt(i,k) |
---|
| 317 | d_q(i,k) = local_q(i,k) - zq(i,k) |
---|
[135] | 318 | ENDDO |
---|
| 319 | ENDDO |
---|
| 320 | |
---|
| 321 | ! now subroutine -----> GCM variables |
---|
| 322 | DO k = 1, nlayermx |
---|
[787] | 323 | DO i = 1, ngrid |
---|
[135] | 324 | |
---|
[728] | 325 | pdtmana(i,k) = d_t(i,k)/ptimestep |
---|
| 326 | pdqmana(i,k,i_h2o) = d_q(i,k)/ptimestep |
---|
| 327 | pdqmana(i,k,i_ice) = d_ql(i,k)/ptimestep |
---|
[135] | 328 | |
---|
| 329 | ENDDO |
---|
| 330 | ENDDO |
---|
| 331 | |
---|
[253] | 332 | |
---|
[135] | 333 | RETURN |
---|
[253] | 334 | END |
---|