[253] | 1 | subroutine convadj(ngrid,nlay,nq,ptimestep, |
---|
| 2 | & pplay,pplev,ppopsk, |
---|
| 3 | & pu,pv,ph,pq, |
---|
| 4 | & pdufi,pdvfi,pdhfi,pdqfi, |
---|
| 5 | & pduadj,pdvadj,pdhadj, |
---|
| 6 | & pdqadj) |
---|
[135] | 7 | |
---|
[253] | 8 | implicit none |
---|
[135] | 9 | |
---|
[253] | 10 | !================================================================== |
---|
| 11 | ! |
---|
| 12 | ! Purpose |
---|
| 13 | ! ------- |
---|
| 14 | ! Calculates dry convective adjustment. If one tracer is CO2, |
---|
| 15 | ! we take into account the molecular mass variation (e.g. when |
---|
| 16 | ! CO2 condenses) to trigger convection (F. Forget 01/2005) |
---|
| 17 | ! |
---|
| 18 | ! Authors |
---|
| 19 | ! ------- |
---|
| 20 | ! Original author unknown. |
---|
| 21 | ! Modif. 2005 by F. Forget. |
---|
| 22 | ! |
---|
| 23 | !================================================================== |
---|
[135] | 24 | |
---|
[253] | 25 | ! ------------ |
---|
| 26 | ! Declarations |
---|
| 27 | ! ------------ |
---|
| 28 | |
---|
[135] | 29 | #include "dimensions.h" |
---|
| 30 | #include "dimphys.h" |
---|
| 31 | #include "comcstfi.h" |
---|
| 32 | #include "callkeys.h" |
---|
| 33 | #include "tracer.h" |
---|
| 34 | |
---|
| 35 | |
---|
[253] | 36 | ! Arguments |
---|
| 37 | ! --------- |
---|
[135] | 38 | |
---|
| 39 | INTEGER ngrid,nlay |
---|
| 40 | REAL ptimestep |
---|
| 41 | REAL ph(ngrid,nlay),pdhfi(ngrid,nlay),pdhadj(ngrid,nlay) |
---|
| 42 | REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1),ppopsk(ngrid,nlay) |
---|
| 43 | REAL pu(ngrid,nlay),pdufi(ngrid,nlay),pduadj(ngrid,nlay) |
---|
| 44 | REAL pv(ngrid,nlay),pdvfi(ngrid,nlay),pdvadj(ngrid,nlay) |
---|
| 45 | |
---|
[253] | 46 | ! Tracers |
---|
[135] | 47 | integer nq |
---|
| 48 | real pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) |
---|
| 49 | real pdqadj(ngrid,nlay,nq) |
---|
| 50 | |
---|
| 51 | |
---|
[253] | 52 | ! Local |
---|
| 53 | ! ----- |
---|
[135] | 54 | |
---|
| 55 | INTEGER ig,i,l,l1,l2,jj |
---|
| 56 | INTEGER jcnt, jadrs(ngridmx) |
---|
| 57 | |
---|
| 58 | REAL sig(nlayermx+1),sdsig(nlayermx),dsig(nlayermx) |
---|
| 59 | REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx) |
---|
| 60 | REAL zh(ngridmx,nlayermx) |
---|
| 61 | REAL zu2(ngridmx,nlayermx),zv2(ngridmx,nlayermx) |
---|
| 62 | REAL zh2(ngridmx,nlayermx), zhc(ngridmx,nlayermx) |
---|
| 63 | REAL zhm,zsm,zdsm,zum,zvm,zalpha,zhmc |
---|
| 64 | |
---|
[253] | 65 | ! Tracers |
---|
[135] | 66 | INTEGER iq,ico2 |
---|
| 67 | save ico2 |
---|
| 68 | REAL zq(ngridmx,nlayermx,nqmx), zq2(ngridmx,nlayermx,nqmx) |
---|
| 69 | REAL zqm(nqmx),zqco2m |
---|
| 70 | real m_co2, m_noco2, A , B |
---|
| 71 | save A, B |
---|
| 72 | |
---|
| 73 | real mtot1, mtot2 , mm1, mm2 |
---|
| 74 | integer l1ref, l2ref |
---|
| 75 | LOGICAL vtest(ngridmx),down,firstcall |
---|
| 76 | save firstcall |
---|
| 77 | data firstcall/.true./ |
---|
| 78 | |
---|
[253] | 79 | ! for conservation test |
---|
| 80 | real masse,cadjncons |
---|
| 81 | |
---|
[135] | 82 | EXTERNAL SCOPY |
---|
[253] | 83 | |
---|
| 84 | ! -------------- |
---|
| 85 | ! Initialisation |
---|
| 86 | ! -------------- |
---|
| 87 | |
---|
[135] | 88 | IF (firstcall) THEN |
---|
| 89 | IF(ngrid.NE.ngridmx) THEN |
---|
| 90 | PRINT* |
---|
[253] | 91 | PRINT*,'STOP in convadj' |
---|
[135] | 92 | PRINT*,'ngrid =',ngrid |
---|
| 93 | PRINT*,'ngridmx =',ngridmx |
---|
| 94 | ENDIF |
---|
| 95 | ico2=0 |
---|
| 96 | if (tracer) then |
---|
[253] | 97 | ! Prepare Special treatment if one of the tracers is CO2 gas |
---|
[135] | 98 | do iq=1,nqmx |
---|
| 99 | if (noms(iq).eq."co2") then |
---|
| 100 | print*,'dont go there' |
---|
| 101 | stop |
---|
| 102 | ico2=iq |
---|
| 103 | m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol) |
---|
| 104 | m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol) |
---|
[253] | 105 | ! Compute A and B coefficient use to compute |
---|
| 106 | ! mean molecular mass Mair defined by |
---|
| 107 | ! 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2 |
---|
| 108 | ! 1/Mair = A*q(ico2) + B |
---|
[135] | 109 | A =(1/m_co2 - 1/m_noco2) |
---|
| 110 | B=1/m_noco2 |
---|
| 111 | end if |
---|
| 112 | enddo |
---|
| 113 | endif |
---|
| 114 | firstcall=.false. |
---|
| 115 | ENDIF ! of IF (firstcall) |
---|
| 116 | |
---|
| 117 | DO l=1,nlay |
---|
| 118 | DO ig=1,ngrid |
---|
| 119 | zh(ig,l)=ph(ig,l)+pdhfi(ig,l)*ptimestep |
---|
| 120 | zu(ig,l)=pu(ig,l)+pdufi(ig,l)*ptimestep |
---|
| 121 | zv(ig,l)=pv(ig,l)+pdvfi(ig,l)*ptimestep |
---|
| 122 | ENDDO |
---|
| 123 | ENDDO |
---|
| 124 | |
---|
| 125 | if(tracer) then |
---|
| 126 | DO iq =1, nq |
---|
| 127 | DO l=1,nlay |
---|
| 128 | DO ig=1,ngrid |
---|
| 129 | zq(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep |
---|
| 130 | ENDDO |
---|
| 131 | ENDDO |
---|
| 132 | ENDDO |
---|
| 133 | end if |
---|
| 134 | |
---|
| 135 | CALL scopy(ngrid*nlay,zh,1,zh2,1) |
---|
| 136 | CALL scopy(ngrid*nlay,zu,1,zu2,1) |
---|
| 137 | CALL scopy(ngrid*nlay,zv,1,zv2,1) |
---|
| 138 | CALL scopy(ngrid*nlay*nq,zq,1,zq2,1) |
---|
| 139 | |
---|
[253] | 140 | ! ----------------------------- |
---|
| 141 | ! Detection of unstable columns |
---|
| 142 | ! ----------------------------- |
---|
| 143 | ! If ph(above) < ph(below) we set vtest=.true. |
---|
| 144 | |
---|
[135] | 145 | DO ig=1,ngrid |
---|
[253] | 146 | vtest(ig)=.false. |
---|
[135] | 147 | ENDDO |
---|
[253] | 148 | |
---|
[135] | 149 | if (ico2.ne.0) then |
---|
[253] | 150 | ! Special case if one of the tracers is CO2 gas |
---|
[135] | 151 | DO l=1,nlay |
---|
| 152 | DO ig=1,ngrid |
---|
| 153 | zhc(ig,l) = zh2(ig,l)*(A*zq2(ig,l,ico2)+B) |
---|
| 154 | ENDDO |
---|
| 155 | ENDDO |
---|
| 156 | else |
---|
| 157 | CALL scopy(ngrid*nlay,zh2,1,zhc,1) |
---|
| 158 | end if |
---|
| 159 | |
---|
[253] | 160 | ! Find out which grid points are convectively unstable |
---|
[135] | 161 | DO l=2,nlay |
---|
| 162 | DO ig=1,ngrid |
---|
[253] | 163 | IF(zhc(ig,l).LT.zhc(ig,l-1)) vtest(ig)=.true. |
---|
[135] | 164 | ENDDO |
---|
| 165 | ENDDO |
---|
[253] | 166 | |
---|
| 167 | ! Make a list of them |
---|
[135] | 168 | jcnt=0 |
---|
| 169 | DO ig=1,ngrid |
---|
| 170 | IF(vtest(ig)) THEN |
---|
| 171 | jcnt=jcnt+1 |
---|
| 172 | jadrs(jcnt)=ig |
---|
| 173 | ENDIF |
---|
| 174 | ENDDO |
---|
| 175 | |
---|
| 176 | |
---|
[253] | 177 | ! --------------------------------------------------------------- |
---|
| 178 | ! Adjustment of the "jcnt" unstable profiles indicated by "jadrs" |
---|
| 179 | ! --------------------------------------------------------------- |
---|
| 180 | |
---|
[135] | 181 | DO jj = 1, jcnt ! loop on every convective grid point |
---|
[253] | 182 | |
---|
[135] | 183 | i = jadrs(jj) |
---|
| 184 | |
---|
[253] | 185 | ! Calculate sigma in this column |
---|
[135] | 186 | DO l=1,nlay+1 |
---|
| 187 | sig(l)=pplev(i,l)/pplev(i,1) |
---|
| 188 | |
---|
| 189 | ENDDO |
---|
| 190 | DO l=1,nlay |
---|
| 191 | dsig(l)=sig(l)-sig(l+1) |
---|
| 192 | sdsig(l)=ppopsk(i,l)*dsig(l) |
---|
| 193 | ENDDO |
---|
| 194 | l2 = 1 |
---|
[253] | 195 | |
---|
| 196 | ! Test loop upwards on l2 |
---|
| 197 | |
---|
[135] | 198 | DO |
---|
| 199 | l2 = l2 + 1 |
---|
| 200 | IF (l2 .GT. nlay) EXIT |
---|
| 201 | IF (zhc(i, l2) .LT. zhc(i, l2-1)) THEN |
---|
| 202 | |
---|
[253] | 203 | ! l2 is the highest level of the unstable column |
---|
[135] | 204 | |
---|
| 205 | l1 = l2 - 1 |
---|
| 206 | l = l1 |
---|
| 207 | zsm = sdsig(l2) |
---|
| 208 | zdsm = dsig(l2) |
---|
| 209 | zhm = zh2(i, l2) |
---|
| 210 | if(ico2.ne.0) zqco2m = zq2(i,l2,ico2) |
---|
[253] | 211 | |
---|
| 212 | ! Test loop downwards |
---|
| 213 | |
---|
[135] | 214 | DO |
---|
| 215 | zsm = zsm + sdsig(l) |
---|
| 216 | zdsm = zdsm + dsig(l) |
---|
| 217 | zhm = zhm + sdsig(l) * (zh2(i, l) - zhm) / zsm |
---|
| 218 | if(ico2.ne.0) then |
---|
| 219 | zqco2m = |
---|
| 220 | & zqco2m + dsig(l) * (zq2(i,l,ico2) - zqco2m) / zdsm |
---|
| 221 | zhmc = zhm*(A*zqco2m+B) |
---|
| 222 | else |
---|
| 223 | zhmc = zhm |
---|
| 224 | end if |
---|
| 225 | |
---|
[253] | 226 | ! do we have to extend the column downwards? |
---|
[135] | 227 | |
---|
[253] | 228 | down = .false. |
---|
| 229 | IF (l1 .ne. 1) then !-- and then |
---|
| 230 | IF (zhmc .lt. zhc(i, l1-1)) then |
---|
| 231 | down = .true. |
---|
[135] | 232 | END IF |
---|
| 233 | END IF |
---|
| 234 | |
---|
[253] | 235 | ! this could be a problem... |
---|
| 236 | |
---|
| 237 | if (down) then |
---|
[135] | 238 | |
---|
| 239 | l1 = l1 - 1 |
---|
| 240 | l = l1 |
---|
| 241 | |
---|
[253] | 242 | else |
---|
[135] | 243 | |
---|
[253] | 244 | ! can we extend the column upwards? |
---|
[135] | 245 | |
---|
[253] | 246 | if (l2 .eq. nlay) exit |
---|
[135] | 247 | |
---|
[253] | 248 | if (zhc(i, l2+1) .ge. zhmc) exit |
---|
| 249 | |
---|
[135] | 250 | l2 = l2 + 1 |
---|
| 251 | l = l2 |
---|
[253] | 252 | |
---|
| 253 | end if |
---|
| 254 | |
---|
| 255 | enddo |
---|
| 256 | |
---|
| 257 | ! New constant profile (average value) |
---|
| 258 | |
---|
| 259 | |
---|
[135] | 260 | zalpha=0. |
---|
| 261 | zum=0. |
---|
| 262 | zvm=0. |
---|
| 263 | do iq=1,nq |
---|
| 264 | zqm(iq) = 0. |
---|
| 265 | end do |
---|
| 266 | DO l = l1, l2 |
---|
| 267 | if(ico2.ne.0) then |
---|
| 268 | zalpha=zalpha+ |
---|
| 269 | & ABS(zhc(i,l)/(A+B*zqco2m) -zhm)*dsig(l) |
---|
| 270 | else |
---|
| 271 | zalpha=zalpha+ABS(zh2(i,l)-zhm)*dsig(l) |
---|
| 272 | endif |
---|
| 273 | zh2(i, l) = zhm |
---|
[253] | 274 | ! modifs by RDW !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 275 | zum=zum+dsig(l)*zu2(i,l) |
---|
| 276 | zvm=zvm+dsig(l)*zv2(i,l) |
---|
| 277 | ! zum=zum+dsig(l)*zu(i,l) |
---|
| 278 | ! zvm=zvm+dsig(l)*zv(i,l) |
---|
[135] | 279 | do iq=1,nq |
---|
[253] | 280 | zqm(iq) = zqm(iq)+dsig(l)*zq2(i,l,iq) |
---|
| 281 | ! zqm(iq) = zqm(iq)+dsig(l)*zq(i,l,iq) |
---|
| 282 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 283 | |
---|
| 284 | ! to conserve tracers/ KE, we must calculate zum, zvm and zqm using |
---|
| 285 | ! the up-to-date column values. If we do not do this, there are cases |
---|
| 286 | ! where convection stops at one level and starts at the next where we |
---|
| 287 | ! can break conservation of stuff (particularly tracers) significantly. |
---|
| 288 | |
---|
[135] | 289 | end do |
---|
| 290 | ENDDO |
---|
| 291 | zalpha=zalpha/(zhm*(sig(l1)-sig(l2+1))) |
---|
| 292 | zum=zum/(sig(l1)-sig(l2+1)) |
---|
| 293 | zvm=zvm/(sig(l1)-sig(l2+1)) |
---|
| 294 | do iq=1,nq |
---|
| 295 | zqm(iq) = zqm(iq)/(sig(l1)-sig(l2+1)) |
---|
| 296 | end do |
---|
| 297 | |
---|
| 298 | IF(zalpha.GT.1.) THEN |
---|
| 299 | zalpha=1. |
---|
| 300 | ELSE |
---|
[253] | 301 | ! IF(zalpha.LT.0.) STOP |
---|
[135] | 302 | IF(zalpha.LT.1.e-4) zalpha=1.e-4 |
---|
| 303 | ENDIF |
---|
| 304 | |
---|
| 305 | DO l=l1,l2 |
---|
| 306 | zu2(i,l)=zu2(i,l)+zalpha*(zum-zu2(i,l)) |
---|
| 307 | zv2(i,l)=zv2(i,l)+zalpha*(zvm-zv2(i,l)) |
---|
| 308 | do iq=1,nq |
---|
[253] | 309 | ! zq2(i,l,iq)=zq2(i,l,iq)+zalpha*(zqm(iq)-zq2(i,l,iq)) |
---|
[135] | 310 | zq2(i,l,iq)=zqm(iq) |
---|
| 311 | end do |
---|
| 312 | ENDDO |
---|
| 313 | if (ico2.ne.0) then |
---|
| 314 | DO l=l1, l2 |
---|
| 315 | zhc(i,l) = zh2(i,l)*(A*zq2(i,l,ico2)+B) |
---|
| 316 | ENDDO |
---|
| 317 | end if |
---|
| 318 | |
---|
[253] | 319 | |
---|
[135] | 320 | l2 = l2 + 1 |
---|
[253] | 321 | |
---|
| 322 | END IF ! End of l1 to l2 instability treatment |
---|
| 323 | ! We now continue to test from l2 upwards |
---|
| 324 | |
---|
| 325 | ENDDO ! End of upwards loop on l2 |
---|
| 326 | |
---|
| 327 | |
---|
| 328 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 329 | ! check conservation |
---|
| 330 | cadjncons=0.0 |
---|
| 331 | if(water)then |
---|
| 332 | do l = 1, nlay |
---|
| 333 | masse = (pplev(i,l) - pplev(i,l+1))/g |
---|
| 334 | iq = igcm_h2o_vap |
---|
| 335 | cadjncons = cadjncons + |
---|
| 336 | & masse*(zq2(i,l,iq)-zq(i,l,iq))/ptimestep |
---|
| 337 | end do |
---|
| 338 | endif |
---|
| 339 | |
---|
| 340 | if(cadjncons.lt.-1.e-6)then |
---|
| 341 | print*,'convadj has just crashed...' |
---|
| 342 | print*,'i = ',i |
---|
| 343 | print*,'l1 = ',l1 |
---|
| 344 | print*,'l2 = ',l2 |
---|
| 345 | print*,'cadjncons = ',cadjncons |
---|
| 346 | do l = 1, nlay |
---|
| 347 | print*,'dsig = ',dsig(l) |
---|
| 348 | end do |
---|
| 349 | do l = 1, nlay |
---|
| 350 | print*,'dsig = ',dsig(l) |
---|
| 351 | end do |
---|
| 352 | do l = 1, nlay |
---|
| 353 | print*,'sig = ',sig(l) |
---|
| 354 | end do |
---|
| 355 | do l = 1, nlay |
---|
| 356 | print*,'pplay(ig,:) = ',pplay(i,l) |
---|
| 357 | end do |
---|
| 358 | do l = 1, nlay+1 |
---|
| 359 | print*,'pplev(ig,:) = ',pplev(i,l) |
---|
| 360 | end do |
---|
| 361 | do l = 1, nlay |
---|
| 362 | print*,'ph(ig,:) = ',ph(i,l) |
---|
| 363 | end do |
---|
| 364 | do l = 1, nlay |
---|
| 365 | print*,'ph(ig,:) = ',ph(i,l) |
---|
| 366 | end do |
---|
| 367 | do l = 1, nlay |
---|
| 368 | print*,'ph(ig,:) = ',ph(i,l) |
---|
| 369 | end do |
---|
| 370 | do l = 1, nlay |
---|
| 371 | print*,'zh(ig,:) = ',zh(i,l) |
---|
| 372 | end do |
---|
| 373 | do l = 1, nlay |
---|
| 374 | print*,'zh2(ig,:) = ',zh2(i,l) |
---|
| 375 | end do |
---|
| 376 | do l = 1, nlay |
---|
| 377 | print*,'zq(ig,:,vap) = ',zq(i,l,igcm_h2o_vap) |
---|
| 378 | end do |
---|
| 379 | do l = 1, nlay |
---|
| 380 | print*,'zq2(ig,:,vap) = ',zq2(i,l,igcm_h2o_vap) |
---|
| 381 | end do |
---|
| 382 | print*,'zqm(vap) = ',zqm(igcm_h2o_vap) |
---|
| 383 | print*,'jadrs=',jadrs |
---|
| 384 | |
---|
| 385 | call abort |
---|
| 386 | endif |
---|
| 387 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 388 | |
---|
| 389 | |
---|
[135] | 390 | ENDDO |
---|
[253] | 391 | |
---|
[135] | 392 | DO l=1,nlay |
---|
| 393 | DO ig=1,ngrid |
---|
| 394 | pdhadj(ig,l)=(zh2(ig,l)-zh(ig,l))/ptimestep |
---|
| 395 | pduadj(ig,l)=(zu2(ig,l)-zu(ig,l))/ptimestep |
---|
| 396 | pdvadj(ig,l)=(zv2(ig,l)-zv(ig,l))/ptimestep |
---|
| 397 | ENDDO |
---|
| 398 | ENDDO |
---|
| 399 | |
---|
| 400 | if(tracer) then |
---|
| 401 | do iq=1, nq |
---|
| 402 | do l=1,nlay |
---|
| 403 | do ig=1, ngrid |
---|
| 404 | pdqadj(ig,l,iq)=(zq2(ig,l,iq)-zq(ig,l,iq))/ptimestep |
---|
| 405 | end do |
---|
| 406 | end do |
---|
| 407 | end do |
---|
| 408 | end if |
---|
| 409 | |
---|
[253] | 410 | |
---|
| 411 | ! output |
---|
[135] | 412 | ! if (ngrid.eq.1) then |
---|
| 413 | ! ig=1 |
---|
| 414 | ! iq =1 |
---|
| 415 | ! write(*,*)'**** l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq)' |
---|
| 416 | ! do l=nlay,1,-1 |
---|
| 417 | ! write(*,*) l, pq(ig,l,iq),zq(ig,l,iq),zq2(ig,l,iq) |
---|
| 418 | ! end do |
---|
| 419 | ! end if |
---|
| 420 | |
---|
| 421 | |
---|
[253] | 422 | return |
---|
| 423 | end |
---|