| 1 | subroutine thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr, & |
|---|
| 2 | & masse,q,dq,qa,lev_out) |
|---|
| 3 | USE print_control_mod, ONLY: prt_level |
|---|
| 4 | implicit none |
|---|
| 5 | |
|---|
| 6 | !======================================================================= |
|---|
| 7 | ! |
|---|
| 8 | ! Calcul du transport verticale dans la couche limite en presence |
|---|
| 9 | ! de "thermiques" explicitement representes |
|---|
| 10 | ! calcul du dq/dt une fois qu'on connait les ascendances |
|---|
| 11 | ! |
|---|
| 12 | ! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr) |
|---|
| 13 | ! Introduction of an implicit computation of vertical advection in |
|---|
| 14 | ! the environment of thermal plumes in thermcell_dq |
|---|
| 15 | ! impl = 0 : explicit, 1 : implicit, -1 : old version |
|---|
| 16 | ! |
|---|
| 17 | !======================================================================= |
|---|
| 18 | |
|---|
| 19 | integer ngrid,nlay,impl |
|---|
| 20 | |
|---|
| 21 | real ptimestep |
|---|
| 22 | real masse(ngrid,nlay),fm(ngrid,nlay+1) |
|---|
| 23 | real entr(ngrid,nlay) |
|---|
| 24 | real q(ngrid,nlay) |
|---|
| 25 | real dq(ngrid,nlay) |
|---|
| 26 | integer lev_out ! niveau pour les print |
|---|
| 27 | |
|---|
| 28 | real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1) |
|---|
| 29 | |
|---|
| 30 | real zzm |
|---|
| 31 | |
|---|
| 32 | integer ig,k |
|---|
| 33 | real cfl |
|---|
| 34 | |
|---|
| 35 | real qold(ngrid,nlay),fqa(ngrid,nlay+1) |
|---|
| 36 | integer niter,iter |
|---|
| 37 | CHARACTER (LEN=20) :: modname='thermcell_dq' |
|---|
| 38 | CHARACTER (LEN=80) :: abort_message |
|---|
| 39 | |
|---|
| 40 | |
|---|
| 41 | ! Old explicite scheme |
|---|
| 42 | if (impl<=-1) then |
|---|
| 43 | call thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr, & |
|---|
| 44 | & masse,q,dq,qa,lev_out) |
|---|
| 45 | return |
|---|
| 46 | endif |
|---|
| 47 | |
|---|
| 48 | ! Calcul du critere CFL pour l'advection dans la subsidence |
|---|
| 49 | cfl = 0. |
|---|
| 50 | do k=1,nlay |
|---|
| 51 | do ig=1,ngrid |
|---|
| 52 | zzm=masse(ig,k)/ptimestep |
|---|
| 53 | cfl=max(cfl,fm(ig,k)/zzm) |
|---|
| 54 | if (entr(ig,k).gt.zzm) then |
|---|
| 55 | print*,'entr*dt>m,1',k,entr(ig,k)*ptimestep,masse(ig,k) |
|---|
| 56 | abort_message = 'entr dt > m, 1st' |
|---|
| 57 | CALL abort_physic (modname,abort_message,1) |
|---|
| 58 | endif |
|---|
| 59 | enddo |
|---|
| 60 | enddo |
|---|
| 61 | |
|---|
| 62 | qold=q |
|---|
| 63 | |
|---|
| 64 | |
|---|
| 65 | if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0' |
|---|
| 66 | |
|---|
| 67 | ! calcul du detrainement |
|---|
| 68 | do k=1,nlay |
|---|
| 69 | do ig=1,ngrid |
|---|
| 70 | detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k) |
|---|
| 71 | ! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k) |
|---|
| 72 | !test |
|---|
| 73 | if (detr(ig,k).lt.0.) then |
|---|
| 74 | entr(ig,k)=entr(ig,k)-detr(ig,k) |
|---|
| 75 | detr(ig,k)=0. |
|---|
| 76 | ! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k), |
|---|
| 77 | ! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k) |
|---|
| 78 | endif |
|---|
| 79 | if (fm(ig,k+1).lt.0.) then |
|---|
| 80 | ! print*,'fm2<0!!!' |
|---|
| 81 | endif |
|---|
| 82 | if (entr(ig,k).lt.0.) then |
|---|
| 83 | ! print*,'entr2<0!!!' |
|---|
| 84 | endif |
|---|
| 85 | enddo |
|---|
| 86 | enddo |
|---|
| 87 | |
|---|
| 88 | ! Computation of tracer concentrations in the ascending plume |
|---|
| 89 | do ig=1,ngrid |
|---|
| 90 | qa(ig,1)=q(ig,1) |
|---|
| 91 | enddo |
|---|
| 92 | |
|---|
| 93 | do k=2,nlay |
|---|
| 94 | do ig=1,ngrid |
|---|
| 95 | if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. & |
|---|
| 96 | & 1.e-5*masse(ig,k)) then |
|---|
| 97 | qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k)) & |
|---|
| 98 | & /(fm(ig,k+1)+detr(ig,k)) |
|---|
| 99 | else |
|---|
| 100 | qa(ig,k)=q(ig,k) |
|---|
| 101 | endif |
|---|
| 102 | if (qa(ig,k).lt.0.) then |
|---|
| 103 | ! print*,'qa<0!!!' |
|---|
| 104 | endif |
|---|
| 105 | if (q(ig,k).lt.0.) then |
|---|
| 106 | ! print*,'q<0!!!' |
|---|
| 107 | endif |
|---|
| 108 | enddo |
|---|
| 109 | enddo |
|---|
| 110 | |
|---|
| 111 | ! Plume vertical flux |
|---|
| 112 | do k=2,nlay-1 |
|---|
| 113 | fqa(:,k)=fm(:,k)*qa(:,k-1) |
|---|
| 114 | enddo |
|---|
| 115 | fqa(:,1)=0. ; fqa(:,nlay)=0. |
|---|
| 116 | |
|---|
| 117 | |
|---|
| 118 | ! Trace species evolution |
|---|
| 119 | if (impl==0) then |
|---|
| 120 | do k=1,nlay-1 |
|---|
| 121 | q(:,k)=q(:,k)+(fqa(:,k)-fqa(:,k+1)-fm(:,k)*q(:,k)+fm(:,k+1)*q(:,k+1)) & |
|---|
| 122 | & *ptimestep/masse(:,k) |
|---|
| 123 | enddo |
|---|
| 124 | else |
|---|
| 125 | do k=nlay-1,1,-1 |
|---|
| 126 | ! FH debut de modif : le calcul ci dessous modifiait numériquement |
|---|
| 127 | ! la concentration quand le flux de masse etait nul car on divisait |
|---|
| 128 | ! puis multipliait par masse/ptimestep. |
|---|
| 129 | ! q(:,k)=(masse(:,k)*q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1)) & |
|---|
| 130 | ! & /(fm(:,k)+masse(:,k)/ptimestep) |
|---|
| 131 | q(:,k)=(q(:,k)+ptimestep/masse(:,k)*(fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1))) & |
|---|
| 132 | & /(1.+fm(:,k)*ptimestep/masse(:,k)) |
|---|
| 133 | ! FH fin de modif. |
|---|
| 134 | enddo |
|---|
| 135 | endif |
|---|
| 136 | |
|---|
| 137 | ! Tendencies |
|---|
| 138 | do k=1,nlay |
|---|
| 139 | do ig=1,ngrid |
|---|
| 140 | dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep |
|---|
| 141 | q(ig,k)=qold(ig,k) |
|---|
| 142 | enddo |
|---|
| 143 | enddo |
|---|
| 144 | |
|---|
| 145 | return |
|---|
| 146 | end |
|---|
| 147 | |
|---|
| 148 | #ifdef ISO |
|---|
| 149 | subroutine thermcell_dq_iso(ngrid,nlay,impl,ptimestep,fm,entr, & |
|---|
| 150 | & masse,q,dq,qa,lev_out,xt,dxt,xta) |
|---|
| 151 | USE infotrac_phy, ONLY : ntraciso |
|---|
| 152 | #ifdef ISOVERIF |
|---|
| 153 | USE isotopes_mod, ONLY: iso_eau,iso_HDO |
|---|
| 154 | USE isotopes_verif_mod, ONLY: iso_verif_egalite, & |
|---|
| 155 | iso_verif_aberrant_enc_vect2D,iso_verif_egalite_vect2D |
|---|
| 156 | #endif |
|---|
| 157 | implicit none |
|---|
| 158 | |
|---|
| 159 | |
|---|
| 160 | #include "iniprint.h" |
|---|
| 161 | !======================================================================= |
|---|
| 162 | ! |
|---|
| 163 | ! Calcul du transport verticale dans la couche limite en presence |
|---|
| 164 | ! de "thermiques" explicitement representes |
|---|
| 165 | ! calcul du dq/dt une fois qu'on connait les ascendances |
|---|
| 166 | ! |
|---|
| 167 | ! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr) |
|---|
| 168 | ! Introduction of an implicit computation of vertical advection in |
|---|
| 169 | ! the environment of thermal plumes in thermcell_dq |
|---|
| 170 | ! impl = 0 : explicit, 1 : implicit, -1 : old version |
|---|
| 171 | ! |
|---|
| 172 | !======================================================================= |
|---|
| 173 | |
|---|
| 174 | integer ngrid,nlay,impl |
|---|
| 175 | |
|---|
| 176 | real ptimestep |
|---|
| 177 | real masse(ngrid,nlay),fm(ngrid,nlay+1) |
|---|
| 178 | real entr(ngrid,nlay) |
|---|
| 179 | real q(ngrid,nlay) |
|---|
| 180 | real dq(ngrid,nlay) |
|---|
| 181 | integer lev_out ! niveau pour les print |
|---|
| 182 | |
|---|
| 183 | real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1) |
|---|
| 184 | |
|---|
| 185 | real zzm |
|---|
| 186 | |
|---|
| 187 | integer ig,k |
|---|
| 188 | real cfl |
|---|
| 189 | |
|---|
| 190 | real qold(ngrid,nlay),fqa(ngrid,nlay+1) |
|---|
| 191 | integer niter,iter |
|---|
| 192 | CHARACTER (LEN=20) :: modname='thermcell_dq' |
|---|
| 193 | CHARACTER (LEN=80) :: abort_message |
|---|
| 194 | |
|---|
| 195 | ! ifdef ISO |
|---|
| 196 | real xt(ntraciso,ngrid,nlay) |
|---|
| 197 | real dxt(ntraciso,ngrid,nlay) |
|---|
| 198 | real xta(ntraciso,ngrid,nlay) |
|---|
| 199 | real q_avant_entr(ngrid,nlay) |
|---|
| 200 | integer ixt |
|---|
| 201 | real wxtd(ntraciso,ngrid,nlay+1) |
|---|
| 202 | real xtold(ntraciso,ngrid,nlay) |
|---|
| 203 | real fxta(ntraciso,ngrid,nlay+1) |
|---|
| 204 | ! endif |
|---|
| 205 | |
|---|
| 206 | ! Old explicite scheme |
|---|
| 207 | if (impl==-1) then |
|---|
| 208 | call thermcell_dq_o_iso(ngrid,nlay,ptimestep,fm,entr, & |
|---|
| 209 | & masse,q,dq,qa,lev_out,xt,dxt,xta) |
|---|
| 210 | return |
|---|
| 211 | endif |
|---|
| 212 | |
|---|
| 213 | ! Calcul du critere CFL pour l'advection dans la subsidence |
|---|
| 214 | cfl = 0. |
|---|
| 215 | do k=1,nlay |
|---|
| 216 | do ig=1,ngrid |
|---|
| 217 | zzm=masse(ig,k)/ptimestep |
|---|
| 218 | cfl=max(cfl,fm(ig,k)/zzm) |
|---|
| 219 | if (entr(ig,k).gt.zzm) then |
|---|
| 220 | print*,'entr dt > m ',entr(ig,k)*ptimestep,masse(ig,k) |
|---|
| 221 | abort_message = '' |
|---|
| 222 | CALL abort_gcm (modname,abort_message,1) |
|---|
| 223 | endif |
|---|
| 224 | enddo |
|---|
| 225 | enddo |
|---|
| 226 | |
|---|
| 227 | qold=q |
|---|
| 228 | !#ifdef ISO |
|---|
| 229 | xtold=xt |
|---|
| 230 | !#endif |
|---|
| 231 | |
|---|
| 232 | if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0' |
|---|
| 233 | |
|---|
| 234 | ! calcul du detrainement |
|---|
| 235 | do k=1,nlay |
|---|
| 236 | do ig=1,ngrid |
|---|
| 237 | detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k) |
|---|
| 238 | ! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k) |
|---|
| 239 | !test |
|---|
| 240 | if (detr(ig,k).lt.0.) then |
|---|
| 241 | entr(ig,k)=entr(ig,k)-detr(ig,k) |
|---|
| 242 | detr(ig,k)=0. |
|---|
| 243 | ! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k), |
|---|
| 244 | ! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k) |
|---|
| 245 | endif |
|---|
| 246 | if (fm(ig,k+1).lt.0.) then |
|---|
| 247 | ! print*,'fm2<0!!!' |
|---|
| 248 | endif |
|---|
| 249 | if (entr(ig,k).lt.0.) then |
|---|
| 250 | ! print*,'entr2<0!!!' |
|---|
| 251 | endif |
|---|
| 252 | enddo |
|---|
| 253 | enddo |
|---|
| 254 | ! bla90 |
|---|
| 255 | |
|---|
| 256 | ! Computation of tracer concentrations in the ascending plume |
|---|
| 257 | do ig=1,ngrid |
|---|
| 258 | qa(ig,1)=q(ig,1) |
|---|
| 259 | enddo |
|---|
| 260 | !#ifdef ISO |
|---|
| 261 | do ig=1,ngrid |
|---|
| 262 | do ixt=1,ntraciso |
|---|
| 263 | xta(ixt,ig,1)=xt(ixt,ig,1) |
|---|
| 264 | enddo |
|---|
| 265 | enddo |
|---|
| 266 | !#endif |
|---|
| 267 | |
|---|
| 268 | do k=2,nlay |
|---|
| 269 | do ig=1,ngrid |
|---|
| 270 | if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. & |
|---|
| 271 | & 1.e-5*masse(ig,k)) then |
|---|
| 272 | qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k)) & |
|---|
| 273 | & /(fm(ig,k+1)+detr(ig,k)) |
|---|
| 274 | !#ifdef ISO |
|---|
| 275 | do ixt=1,ntraciso |
|---|
| 276 | xta(ixt,ig,k)=(fm(ig,k)*xta(ixt,ig,k-1) & |
|---|
| 277 | & +entr(ig,k)*xt(ixt,ig,k)) & |
|---|
| 278 | & /(fm(ig,k+1)+detr(ig,k)) |
|---|
| 279 | enddo |
|---|
| 280 | !#endif |
|---|
| 281 | else |
|---|
| 282 | qa(ig,k)=q(ig,k) |
|---|
| 283 | !#ifdef ISO |
|---|
| 284 | do ixt=1,ntraciso |
|---|
| 285 | xta(ixt,ig,k)=xt(ixt,ig,k) |
|---|
| 286 | enddo |
|---|
| 287 | !#endif |
|---|
| 288 | endif |
|---|
| 289 | if (qa(ig,k).lt.0.) then |
|---|
| 290 | ! print*,'qa<0!!!' |
|---|
| 291 | endif |
|---|
| 292 | if (q(ig,k).lt.0.) then |
|---|
| 293 | ! print*,'q<0!!!' |
|---|
| 294 | endif |
|---|
| 295 | enddo |
|---|
| 296 | enddo |
|---|
| 297 | |
|---|
| 298 | !#ifdef ISO |
|---|
| 299 | #ifdef ISOVERIF |
|---|
| 300 | if (iso_HDO.gt.0) then |
|---|
| 301 | call iso_verif_aberrant_enc_vect2D( & |
|---|
| 302 | & xta,qa, & |
|---|
| 303 | & 'thermcell_dq_iso 320, qa',ntraciso,ngrid,nlay) |
|---|
| 304 | endif |
|---|
| 305 | #endif |
|---|
| 306 | !#endif |
|---|
| 307 | ! |
|---|
| 308 | ! Plume vertical flux |
|---|
| 309 | do k=2,nlay-1 |
|---|
| 310 | fqa(:,k)=fm(:,k)*qa(:,k-1) |
|---|
| 311 | enddo |
|---|
| 312 | fqa(:,1)=0. ; fqa(:,nlay)=0. |
|---|
| 313 | !#ifdef ISO |
|---|
| 314 | do k=2,nlay-1 |
|---|
| 315 | do ixt=1,ntraciso |
|---|
| 316 | fxta(ixt,:,k)=fm(:,k)*xta(ixt,:,k-1) |
|---|
| 317 | enddo |
|---|
| 318 | enddo |
|---|
| 319 | fxta(ixt,:,1)=0. ; fxta(ixt,:,nlay)=0. |
|---|
| 320 | !#endif |
|---|
| 321 | |
|---|
| 322 | |
|---|
| 323 | ! Trace species evolution |
|---|
| 324 | if (impl==0) then |
|---|
| 325 | do k=1,nlay-1 |
|---|
| 326 | q(:,k)=q(:,k)+(fqa(:,k)-fqa(:,k+1)-fm(:,k)*q(:,k)+fm(:,k+1)*q(:,k+1)) & |
|---|
| 327 | & *ptimestep/masse(:,k) |
|---|
| 328 | !#ifdef ISO |
|---|
| 329 | do ixt=1,ntraciso |
|---|
| 330 | xt(ixt,:,k)=xt(ixt,:,k)+(fxta(ixt,:,k)-fxta(ixt,:,k+1) & |
|---|
| 331 | & -fm(:,k)*xt(ixt,:,k)+fm(:,k+1)*xt(ixt,:,k+1)) & |
|---|
| 332 | & *ptimestep/masse(:,k) |
|---|
| 333 | enddo ! do ixt=1,ntraciso |
|---|
| 334 | !#endif |
|---|
| 335 | enddo !do k=1,nlay-1 |
|---|
| 336 | else !if (impl==0) then |
|---|
| 337 | do k=nlay-1,1,-1 |
|---|
| 338 | q(:,k)=(masse(:,k)*q(:,k)/ptimestep+fqa(:,k)-fqa(:,k+1)+fm(:,k+1)*q(:,k+1)) & |
|---|
| 339 | & /(fm(:,k)+masse(:,k)/ptimestep) |
|---|
| 340 | !#ifdef ISO |
|---|
| 341 | do ixt=1,ntraciso |
|---|
| 342 | xt(ixt,:,k)=xt(ixt,:,k)+(fxta(ixt,:,k)-fxta(ixt,:,k+1) & |
|---|
| 343 | & -fm(:,k)*xt(ixt,:,k)+fm(:,k+1)*xt(ixt,:,k+1)) & |
|---|
| 344 | & *ptimestep/masse(:,k) |
|---|
| 345 | enddo ! do ixt=1,ntraciso |
|---|
| 346 | !#endif |
|---|
| 347 | enddo !do k=nlay-1,1,-1 |
|---|
| 348 | endif !if (impl==0) then |
|---|
| 349 | |
|---|
| 350 | ! Tendencies |
|---|
| 351 | do k=1,nlay |
|---|
| 352 | do ig=1,ngrid |
|---|
| 353 | dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep |
|---|
| 354 | q(ig,k)=qold(ig,k) |
|---|
| 355 | !#ifdef ISO |
|---|
| 356 | do ixt=1,ntraciso |
|---|
| 357 | dxt(ixt,ig,k)=(xt(ixt,ig,k)-xtold(ixt,ig,k))/ptimestep |
|---|
| 358 | xt(ixt,ig,k)=xtold(ixt,ig,k) |
|---|
| 359 | enddo ! do ixt=1,ntraciso |
|---|
| 360 | !#endif |
|---|
| 361 | enddo |
|---|
| 362 | enddo |
|---|
| 363 | |
|---|
| 364 | #ifdef ISOVERIF |
|---|
| 365 | if (iso_HDO.gt.0) then |
|---|
| 366 | call iso_verif_aberrant_enc_vect2D( & |
|---|
| 367 | & xt,q, & |
|---|
| 368 | & 'thermcell_dq_iso 219, q',ntraciso,ngrid,nlay) |
|---|
| 369 | endif |
|---|
| 370 | if (iso_eau.gt.0) then |
|---|
| 371 | call iso_verif_egalite_vect2D(xt,q, & |
|---|
| 372 | & 'thermcell_dq_iso 223, q',ntraciso,ngrid,nlay) |
|---|
| 373 | call iso_verif_egalite_vect2D(dxt,dq, & |
|---|
| 374 | & 'thermcell_dq_iso 224, dq',ntraciso,ngrid,nlay) |
|---|
| 375 | endif |
|---|
| 376 | #endif |
|---|
| 377 | |
|---|
| 378 | return |
|---|
| 379 | end |
|---|
| 380 | |
|---|
| 381 | #endif |
|---|
| 382 | ! ifdef iso |
|---|
| 383 | ! |
|---|
| 384 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 385 | ! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations |
|---|
| 386 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
|---|
| 387 | |
|---|
| 388 | subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr, & |
|---|
| 389 | & masse,q,dq,qa,lev_out) |
|---|
| 390 | USE print_control_mod, ONLY: prt_level |
|---|
| 391 | implicit none |
|---|
| 392 | |
|---|
| 393 | !======================================================================= |
|---|
| 394 | ! |
|---|
| 395 | ! Calcul du transport verticale dans la couche limite en presence |
|---|
| 396 | ! de "thermiques" explicitement representes |
|---|
| 397 | ! calcul du dq/dt une fois qu'on connait les ascendances |
|---|
| 398 | ! |
|---|
| 399 | !======================================================================= |
|---|
| 400 | |
|---|
| 401 | integer ngrid,nlay,impl |
|---|
| 402 | |
|---|
| 403 | real ptimestep |
|---|
| 404 | real masse(ngrid,nlay),fm(ngrid,nlay+1) |
|---|
| 405 | real entr(ngrid,nlay) |
|---|
| 406 | real q(ngrid,nlay) |
|---|
| 407 | real dq(ngrid,nlay) |
|---|
| 408 | integer lev_out ! niveau pour les print |
|---|
| 409 | |
|---|
| 410 | real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1) |
|---|
| 411 | |
|---|
| 412 | real zzm |
|---|
| 413 | |
|---|
| 414 | integer ig,k |
|---|
| 415 | real cfl |
|---|
| 416 | |
|---|
| 417 | real qold(ngrid,nlay) |
|---|
| 418 | real ztimestep |
|---|
| 419 | integer niter,iter |
|---|
| 420 | CHARACTER (LEN=20) :: modname='thermcell_dq' |
|---|
| 421 | CHARACTER (LEN=80) :: abort_message |
|---|
| 422 | |
|---|
| 423 | |
|---|
| 424 | |
|---|
| 425 | ! Calcul du critere CFL pour l'advection dans la subsidence |
|---|
| 426 | cfl = 0. |
|---|
| 427 | do k=1,nlay |
|---|
| 428 | do ig=1,ngrid |
|---|
| 429 | zzm=masse(ig,k)/ptimestep |
|---|
| 430 | cfl=max(cfl,fm(ig,k)/zzm) |
|---|
| 431 | if (entr(ig,k).gt.zzm) then |
|---|
| 432 | print*,'entr*dt>m,2',k,entr(ig,k)*ptimestep,masse(ig,k) |
|---|
| 433 | abort_message = 'entr dt > m, 2nd' |
|---|
| 434 | CALL abort_physic (modname,abort_message,1) |
|---|
| 435 | endif |
|---|
| 436 | enddo |
|---|
| 437 | enddo |
|---|
| 438 | |
|---|
| 439 | !IM 090508 print*,'CFL CFL CFL CFL ',cfl |
|---|
| 440 | |
|---|
| 441 | #undef CFL |
|---|
| 442 | #ifdef CFL |
|---|
| 443 | ! On subdivise le calcul en niter pas de temps. |
|---|
| 444 | niter=int(cfl)+1 |
|---|
| 445 | #else |
|---|
| 446 | niter=1 |
|---|
| 447 | #endif |
|---|
| 448 | |
|---|
| 449 | ztimestep=ptimestep/niter |
|---|
| 450 | qold=q |
|---|
| 451 | |
|---|
| 452 | |
|---|
| 453 | do iter=1,niter |
|---|
| 454 | if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0' |
|---|
| 455 | |
|---|
| 456 | ! calcul du detrainement |
|---|
| 457 | do k=1,nlay |
|---|
| 458 | do ig=1,ngrid |
|---|
| 459 | detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k) |
|---|
| 460 | ! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k) |
|---|
| 461 | !test |
|---|
| 462 | if (detr(ig,k).lt.0.) then |
|---|
| 463 | entr(ig,k)=entr(ig,k)-detr(ig,k) |
|---|
| 464 | detr(ig,k)=0. |
|---|
| 465 | ! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k), |
|---|
| 466 | ! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k) |
|---|
| 467 | endif |
|---|
| 468 | if (fm(ig,k+1).lt.0.) then |
|---|
| 469 | ! print*,'fm2<0!!!' |
|---|
| 470 | endif |
|---|
| 471 | if (entr(ig,k).lt.0.) then |
|---|
| 472 | ! print*,'entr2<0!!!' |
|---|
| 473 | endif |
|---|
| 474 | enddo |
|---|
| 475 | enddo |
|---|
| 476 | |
|---|
| 477 | ! calcul de la valeur dans les ascendances |
|---|
| 478 | do ig=1,ngrid |
|---|
| 479 | qa(ig,1)=q(ig,1) |
|---|
| 480 | enddo |
|---|
| 481 | |
|---|
| 482 | do k=2,nlay |
|---|
| 483 | do ig=1,ngrid |
|---|
| 484 | if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt. & |
|---|
| 485 | & 1.e-5*masse(ig,k)) then |
|---|
| 486 | qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k)) & |
|---|
| 487 | & /(fm(ig,k+1)+detr(ig,k)) |
|---|
| 488 | else |
|---|
| 489 | qa(ig,k)=q(ig,k) |
|---|
| 490 | endif |
|---|
| 491 | if (qa(ig,k).lt.0.) then |
|---|
| 492 | ! print*,'qa<0!!!' |
|---|
| 493 | endif |
|---|
| 494 | if (q(ig,k).lt.0.) then |
|---|
| 495 | ! print*,'q<0!!!' |
|---|
| 496 | endif |
|---|
| 497 | enddo |
|---|
| 498 | enddo |
|---|
| 499 | |
|---|
| 500 | ! Calcul du flux subsident |
|---|
| 501 | |
|---|
| 502 | do k=2,nlay |
|---|
| 503 | do ig=1,ngrid |
|---|
| 504 | #undef centre |
|---|
| 505 | #ifdef centre |
|---|
| 506 | wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k)) |
|---|
| 507 | #else |
|---|
| 508 | |
|---|
| 509 | #define plusqueun |
|---|
| 510 | #ifdef plusqueun |
|---|
| 511 | ! Schema avec advection sur plus qu'une maille. |
|---|
| 512 | zzm=masse(ig,k)/ztimestep |
|---|
| 513 | if (fm(ig,k)>zzm) then |
|---|
| 514 | wqd(ig,k)=zzm*q(ig,k)+(fm(ig,k)-zzm)*q(ig,k+1) |
|---|
| 515 | else |
|---|
| 516 | wqd(ig,k)=fm(ig,k)*q(ig,k) |
|---|
| 517 | endif |
|---|
| 518 | #else |
|---|
| 519 | wqd(ig,k)=fm(ig,k)*q(ig,k) |
|---|
| 520 | #endif |
|---|
| 521 | #endif |
|---|
| 522 | |
|---|
| 523 | if (wqd(ig,k).lt.0.) then |
|---|
| 524 | ! print*,'wqd<0!!!' |
|---|
| 525 | endif |
|---|
| 526 | enddo |
|---|
| 527 | enddo |
|---|
| 528 | do ig=1,ngrid |
|---|
| 529 | wqd(ig,1)=0. |
|---|
| 530 | wqd(ig,nlay+1)=0. |
|---|
| 531 | enddo |
|---|
| 532 | |
|---|
| 533 | |
|---|
| 534 | ! Calcul des tendances |
|---|
| 535 | do k=1,nlay |
|---|
| 536 | do ig=1,ngrid |
|---|
| 537 | q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k) & |
|---|
| 538 | & -wqd(ig,k)+wqd(ig,k+1)) & |
|---|
| 539 | & *ztimestep/masse(ig,k) |
|---|
| 540 | ! if (dq(ig,k).lt.0.) then |
|---|
| 541 | ! print*,'dq<0!!!' |
|---|
| 542 | ! endif |
|---|
| 543 | enddo |
|---|
| 544 | enddo |
|---|
| 545 | |
|---|
| 546 | |
|---|
| 547 | enddo |
|---|
| 548 | |
|---|
| 549 | |
|---|
| 550 | ! Calcul des tendances |
|---|
| 551 | do k=1,nlay |
|---|
| 552 | do ig=1,ngrid |
|---|
| 553 | dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep |
|---|
| 554 | q(ig,k)=qold(ig,k) |
|---|
| 555 | enddo |
|---|
| 556 | enddo |
|---|
| 557 | |
|---|
| 558 | return |
|---|
| 559 | end |
|---|
| 560 | |
|---|
| 561 | |
|---|
| 562 | #ifdef ISO |
|---|
| 563 | subroutine thermcell_dq_o_iso(ngrid,nlay,ptimestep,fm,entr, & |
|---|
| 564 | & masse,q,dq,qa,lev_out,xt,dxt,xta) |
|---|
| 565 | #ifdef ISO |
|---|
| 566 | use infotrac_phy, ONLY: ntraciso |
|---|
| 567 | #ifdef ISOVERIF |
|---|
| 568 | USE isotopes_mod, ONLY: iso_eau,iso_HDO |
|---|
| 569 | USE isotopes_verif_mod, ONLY: iso_verif_egalite, & |
|---|
| 570 | iso_verif_aberrant_enc_vect2D,iso_verif_egalite_vect2D |
|---|
| 571 | #endif |
|---|
| 572 | #endif |
|---|
| 573 | implicit none |
|---|
| 574 | |
|---|
| 575 | #include "iniprint.h" |
|---|
| 576 | !======================================================================= |
|---|
| 577 | ! |
|---|
| 578 | ! Calcul du transport verticale dans la couche limite en presence |
|---|
| 579 | ! de "thermiques" explicitement representes |
|---|
| 580 | ! calcul du dq/dt une fois qu'on connait les ascendances |
|---|
| 581 | ! |
|---|
| 582 | !======================================================================= |
|---|
| 583 | |
|---|
| 584 | integer ngrid,nlay |
|---|
| 585 | |
|---|
| 586 | real ptimestep |
|---|
| 587 | real masse(ngrid,nlay),fm(ngrid,nlay+1) |
|---|
| 588 | real entr(ngrid,nlay) |
|---|
| 589 | real q(ngrid,nlay) |
|---|
| 590 | real dq(ngrid,nlay) |
|---|
| 591 | integer lev_out ! niveau pour les print |
|---|
| 592 | |
|---|
| 593 | real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1) |
|---|
| 594 | |
|---|
| 595 | real zzm |
|---|
| 596 | |
|---|
| 597 | integer ig,k |
|---|
| 598 | real cfl |
|---|
| 599 | |
|---|
| 600 | real qold(ngrid,nlay) |
|---|
| 601 | real ztimestep |
|---|
| 602 | integer niter,iter |
|---|
| 603 | CHARACTER (LEN=20) :: modname='thermcell_dq' |
|---|
| 604 | CHARACTER (LEN=80) :: abort_message |
|---|
| 605 | |
|---|
| 606 | ! ifdef ISO |
|---|
| 607 | real xt(ntraciso,ngrid,nlay) |
|---|
| 608 | real dxt(ntraciso,ngrid,nlay) |
|---|
| 609 | real xta(ntraciso,ngrid,nlay) |
|---|
| 610 | real xtold(ntraciso,ngrid,nlay) |
|---|
| 611 | real q_avant_entr(ngrid,nlay) |
|---|
| 612 | integer ixt |
|---|
| 613 | real wxtd(ntraciso,ngrid,nlay+1) |
|---|
| 614 | ! endif |
|---|
| 615 | |
|---|
| 616 | ! Calcul du critere CFL pour l'advection dans la subsidence |
|---|
| 617 | cfl = 0. |
|---|
| 618 | do k=1,nlay |
|---|
| 619 | do ig=1,ngrid |
|---|
| 620 | zzm=masse(ig,k)/ptimestep |
|---|
| 621 | cfl=max(cfl,fm(ig,k)/zzm) |
|---|
| 622 | if (entr(ig,k).gt.zzm) then |
|---|
| 623 | print*,'entr dt > m ',entr(ig,k)*ptimestep,masse(ig,k) |
|---|
| 624 | abort_message = '' |
|---|
| 625 | CALL abort_gcm (modname,abort_message,1) |
|---|
| 626 | endif |
|---|
| 627 | enddo |
|---|
| 628 | enddo |
|---|
| 629 | |
|---|
| 630 | !IM 090508 print*,'CFL CFL CFL CFL ',cfl |
|---|
| 631 | |
|---|
| 632 | #undef CFL |
|---|
| 633 | #ifdef CFL |
|---|
| 634 | ! On subdivise le calcul en niter pas de temps. |
|---|
| 635 | niter=int(cfl)+1 |
|---|
| 636 | #else |
|---|
| 637 | niter=1 |
|---|
| 638 | #endif |
|---|
| 639 | |
|---|
| 640 | ztimestep=ptimestep/niter |
|---|
| 641 | qold=q |
|---|
| 642 | !#ifdef ISO |
|---|
| 643 | xtold=xt |
|---|
| 644 | !#endif |
|---|
| 645 | |
|---|
| 646 | do iter=1,niter |
|---|
| 647 | if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0' |
|---|
| 648 | |
|---|
| 649 | ! calcul du detrainement |
|---|
| 650 | do k=1,nlay |
|---|
| 651 | do ig=1,ngrid |
|---|
| 652 | detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k) |
|---|
| 653 | ! print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k) |
|---|
| 654 | !test |
|---|
| 655 | if (detr(ig,k).lt.0.) then |
|---|
| 656 | entr(ig,k)=entr(ig,k)-detr(ig,k) |
|---|
| 657 | detr(ig,k)=0. |
|---|
| 658 | ! print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k), |
|---|
| 659 | ! s 'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k) |
|---|
| 660 | endif |
|---|
| 661 | if (fm(ig,k+1).lt.0.) then |
|---|
| 662 | ! print*,'fm2<0!!!' |
|---|
| 663 | endif |
|---|
| 664 | if (entr(ig,k).lt.0.) then |
|---|
| 665 | ! print*,'entr2<0!!!' |
|---|
| 666 | endif |
|---|
| 667 | enddo |
|---|
| 668 | enddo |
|---|
| 669 | |
|---|
| 670 | ! calcul de la valeur dans les ascendances |
|---|
| 671 | do ig=1,ngrid |
|---|
| 672 | qa(ig,1)=q(ig,1) |
|---|
| 673 | enddo |
|---|
| 674 | !#ifdef ISO |
|---|
| 675 | do ig=1,ngrid |
|---|
| 676 | do ixt=1,ntraciso |
|---|
| 677 | xta(ixt,ig,1)=xt(ixt,ig,1) |
|---|
| 678 | enddo |
|---|
| 679 | enddo |
|---|
| 680 | !#endif |
|---|
| 681 | |
|---|
| 682 | do k=2,nlay |
|---|
| 683 | do ig=1,ngrid |
|---|
| 684 | if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt. & |
|---|
| 685 | & 1.e-5*masse(ig,k)) then |
|---|
| 686 | qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k)) & |
|---|
| 687 | & /(fm(ig,k+1)+detr(ig,k)) |
|---|
| 688 | !#ifdef ISO |
|---|
| 689 | do ixt=1,ntraciso |
|---|
| 690 | xta(ixt,ig,k)=(fm(ig,k)*xta(ixt,ig,k-1) & |
|---|
| 691 | & +entr(ig,k)*xt(ixt,ig,k)) & |
|---|
| 692 | & /(fm(ig,k+1)+detr(ig,k)) |
|---|
| 693 | enddo |
|---|
| 694 | !#endif |
|---|
| 695 | else |
|---|
| 696 | qa(ig,k)=q(ig,k) |
|---|
| 697 | !#ifdef ISO |
|---|
| 698 | do ixt=1,ntraciso |
|---|
| 699 | xta(ixt,ig,k)=xt(ixt,ig,k) |
|---|
| 700 | enddo |
|---|
| 701 | !#endif |
|---|
| 702 | endif |
|---|
| 703 | if (qa(ig,k).lt.0.) then |
|---|
| 704 | ! print*,'qa<0!!!' |
|---|
| 705 | endif |
|---|
| 706 | if (q(ig,k).lt.0.) then |
|---|
| 707 | ! print*,'q<0!!!' |
|---|
| 708 | endif |
|---|
| 709 | enddo |
|---|
| 710 | enddo |
|---|
| 711 | |
|---|
| 712 | !#ifdef ISO |
|---|
| 713 | #ifdef ISOVERIF |
|---|
| 714 | if (iso_HDO.gt.0) then |
|---|
| 715 | call iso_verif_aberrant_enc_vect2D( & |
|---|
| 716 | & xta,qa, & |
|---|
| 717 | & 'thermcell_dq_iso 320, qa',ntraciso,ngrid,nlay) |
|---|
| 718 | endif |
|---|
| 719 | #endif |
|---|
| 720 | !#endif |
|---|
| 721 | |
|---|
| 722 | ! Calcul du flux subsident |
|---|
| 723 | |
|---|
| 724 | do k=2,nlay |
|---|
| 725 | do ig=1,ngrid |
|---|
| 726 | #undef centre |
|---|
| 727 | #ifdef centre |
|---|
| 728 | wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k)) |
|---|
| 729 | !#ifdef ISO |
|---|
| 730 | ! attention à l'advection: on suppose que R varie linéairement entre les |
|---|
| 731 | ! couches, comme q. |
|---|
| 732 | do ixt=1,ntraciso |
|---|
| 733 | wxtd(ixt,ig,k)=wqd(ig,k)*0.5*( & |
|---|
| 734 | & xt(ixt,ig,k-1)/q(ig,k-1)+xt(ixt,ig,k)/q(ig,k)) |
|---|
| 735 | enddo !do ixt=1,ntraciso |
|---|
| 736 | !#endif |
|---|
| 737 | #else |
|---|
| 738 | |
|---|
| 739 | #define plusqueun |
|---|
| 740 | #ifdef plusqueun |
|---|
| 741 | ! Schema avec advection sur plus qu'une maille. |
|---|
| 742 | zzm=masse(ig,k)/ztimestep |
|---|
| 743 | if (fm(ig,k)>zzm) then |
|---|
| 744 | wqd(ig,k)=zzm*q(ig,k)+(fm(ig,k)-zzm)*q(ig,k+1) |
|---|
| 745 | !#ifdef ISO |
|---|
| 746 | do ixt=1,ntraciso |
|---|
| 747 | wxtd(ixt,ig,k)=zzm*xt(ixt,ig,k)+(fm(ig,k)-zzm)*xt(ixt,ig,k+1) |
|---|
| 748 | enddo !do ixt=1,ntraciso |
|---|
| 749 | !#endif |
|---|
| 750 | else |
|---|
| 751 | wqd(ig,k)=fm(ig,k)*q(ig,k) |
|---|
| 752 | !#ifdef ISO |
|---|
| 753 | do ixt=1,ntraciso |
|---|
| 754 | wxtd(ixt,ig,k)=fm(ig,k)*xt(ixt,ig,k) |
|---|
| 755 | enddo !do ixt=1,ntraciso |
|---|
| 756 | !#endif |
|---|
| 757 | endif |
|---|
| 758 | #else |
|---|
| 759 | wqd(ig,k)=fm(ig,k)*q(ig,k) |
|---|
| 760 | !#ifdef ISO |
|---|
| 761 | do ixt=1,ntraciso |
|---|
| 762 | wxtd(ixt,ig,k)=fm(ig,k)*xt(ixt,ig,k) |
|---|
| 763 | enddo !do ixt=1,ntraciso |
|---|
| 764 | !#endif |
|---|
| 765 | #endif |
|---|
| 766 | #endif |
|---|
| 767 | |
|---|
| 768 | if (wqd(ig,k).lt.0.) then |
|---|
| 769 | ! print*,'wqd<0!!!' |
|---|
| 770 | endif |
|---|
| 771 | enddo |
|---|
| 772 | enddo |
|---|
| 773 | do ig=1,ngrid |
|---|
| 774 | wqd(ig,1)=0. |
|---|
| 775 | wqd(ig,nlay+1)=0. |
|---|
| 776 | !#ifdef ISO |
|---|
| 777 | do ixt=1,ntraciso |
|---|
| 778 | wxtd(ixt,ig,1)=0.0 |
|---|
| 779 | wxtd(ixt,ig,nlay+1)=0. |
|---|
| 780 | enddo !do ixt=1,ntraciso |
|---|
| 781 | !#endif |
|---|
| 782 | enddo |
|---|
| 783 | |
|---|
| 784 | |
|---|
| 785 | ! Calcul des tendances |
|---|
| 786 | do k=1,nlay |
|---|
| 787 | do ig=1,ngrid |
|---|
| 788 | !#ifdef ISO |
|---|
| 789 | q_avant_entr(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k) & |
|---|
| 790 | & -wqd(ig,k)+wqd(ig,k+1)) & |
|---|
| 791 | & *ztimestep/masse(ig,k) |
|---|
| 792 | !#endif |
|---|
| 793 | q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k) & |
|---|
| 794 | & -wqd(ig,k)+wqd(ig,k+1)) & |
|---|
| 795 | & *ztimestep/masse(ig,k) |
|---|
| 796 | ! if (dq(ig,k).lt.0.) then |
|---|
| 797 | ! print*,'dq<0!!!' |
|---|
| 798 | ! endif |
|---|
| 799 | enddo |
|---|
| 800 | enddo |
|---|
| 801 | |
|---|
| 802 | !#ifdef ISO |
|---|
| 803 | ! Hypothèse: on ajoute le détrainement, et on fait la subsidence. Ensuite, on |
|---|
| 804 | ! retranche la perte vers l'entrainement en supposant que c'est un puit non |
|---|
| 805 | ! fractionnant |
|---|
| 806 | do k=1,nlay |
|---|
| 807 | do ig=1,ngrid |
|---|
| 808 | do ixt=1,ntraciso |
|---|
| 809 | xt(ixt,ig,k)=xt(ixt,ig,k)+(detr(ig,k)*xta(ixt,ig,k) & |
|---|
| 810 | & -wxtd(ixt,ig,k)+wxtd(ixt,ig,k+1)) & |
|---|
| 811 | & *ztimestep/masse(ig,k) |
|---|
| 812 | xt(ixt,ig,k)=xt(ixt,ig,k)/q_avant_entr(ig,k)*q(ig,k) |
|---|
| 813 | enddo !do ixt=1,ntraciso |
|---|
| 814 | enddo !do ig=1,ngrid |
|---|
| 815 | enddo !do k=1,nlay |
|---|
| 816 | #ifdef ISOVERIF |
|---|
| 817 | if (iso_eau.gt.0) then |
|---|
| 818 | call iso_verif_egalite_vect2D( & |
|---|
| 819 | & xt,q, & |
|---|
| 820 | & 'thermcell_dq_iso 421',ntraciso,ngrid,nlay) |
|---|
| 821 | endif |
|---|
| 822 | if (iso_HDO.gt.0) then |
|---|
| 823 | call iso_verif_aberrant_enc_vect2D( & |
|---|
| 824 | & xt,q, & |
|---|
| 825 | & 'thermcell_dq_iso 426',ntraciso,ngrid,nlay) |
|---|
| 826 | endif |
|---|
| 827 | #endif |
|---|
| 828 | !#endif |
|---|
| 829 | |
|---|
| 830 | enddo |
|---|
| 831 | |
|---|
| 832 | |
|---|
| 833 | ! Calcul des tendances |
|---|
| 834 | do k=1,nlay |
|---|
| 835 | do ig=1,ngrid |
|---|
| 836 | dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep |
|---|
| 837 | q(ig,k)=qold(ig,k) |
|---|
| 838 | !#ifdef ISO |
|---|
| 839 | do ixt=1,ntraciso |
|---|
| 840 | dxt(ixt,ig,k)=(xt(ixt,ig,k)-xtold(ixt,ig,k))/ptimestep |
|---|
| 841 | xt(ixt,ig,k)=xtold(ixt,ig,k) |
|---|
| 842 | enddo ! do ixt=1,ntraciso |
|---|
| 843 | !#endif |
|---|
| 844 | enddo |
|---|
| 845 | enddo |
|---|
| 846 | |
|---|
| 847 | return |
|---|
| 848 | end |
|---|
| 849 | |
|---|
| 850 | #endif |
|---|
| 851 | ! ifdef iso |
|---|
