[185] | 1 | subroutine thermcell_dqupdown(ngrid,nlayer,ptimestep,fm0,entr0, & |
---|
[161] | 2 | & detr0,masse0,q_therm,dq_therm,ztvd,fm_down,ztv,charvar,lmax) |
---|
| 3 | implicit none |
---|
| 4 | |
---|
| 5 | ! #include "iniprint.h" |
---|
| 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 | ! Version modifiee pour prendre les downdrafts a la place de la |
---|
| 12 | ! subsidence compensatoire |
---|
| 13 | !======================================================================= |
---|
| 14 | |
---|
[185] | 15 | #include "dimensions.h" |
---|
| 16 | #include "dimphys.h" |
---|
| 17 | |
---|
[161] | 18 | ! ============================ INPUTS ============================ |
---|
| 19 | |
---|
[185] | 20 | INTEGER, INTENT(IN) :: ngrid,nlayer |
---|
[161] | 21 | REAL, INTENT(IN) :: ptimestep |
---|
[185] | 22 | REAL, INTENT(IN) :: fm0(ngridmx,nlayermx+1) |
---|
| 23 | REAL, INTENT(IN) :: entr0(ngridmx,nlayermx),detr0(ngridmx,nlayermx) |
---|
| 24 | REAL, INTENT(IN) :: q_therm(ngridmx,nlayermx) |
---|
| 25 | REAL, INTENT(IN) :: fm_down(ngridmx,nlayermx+1) |
---|
| 26 | REAL, INTENT(IN) :: ztvd(ngridmx,nlayermx) |
---|
| 27 | REAL, INTENT(IN) :: ztv(ngridmx,nlayermx) |
---|
[161] | 28 | CHARACTER (LEN=20), INTENT(IN) :: charvar |
---|
[185] | 29 | REAL, INTENT(IN) :: masse0(ngridmx,nlayermx) |
---|
| 30 | INTEGER, INTENT(IN) :: lmax(ngridmx) |
---|
[161] | 31 | |
---|
| 32 | ! ============================ OUTPUTS =========================== |
---|
| 33 | |
---|
[185] | 34 | REAL, INTENT(OUT) :: dq_therm(ngridmx,nlayermx) |
---|
[161] | 35 | |
---|
| 36 | ! ============================ LOCAL ============================= |
---|
| 37 | |
---|
[185] | 38 | ! REAL detr0(ngridmx,nlayermx) |
---|
| 39 | REAL detrd(ngridmx,nlayermx) |
---|
| 40 | REAL entrd(ngridmx,nlayermx) |
---|
| 41 | REAL fmd(ngridmx,nlayermx+1) |
---|
| 42 | REAL q(ngridmx,nlayermx) |
---|
| 43 | REAL qa(ngridmx,nlayermx) |
---|
| 44 | REAL qd(ngridmx,nlayermx) |
---|
[161] | 45 | INTEGER ig,k |
---|
[185] | 46 | LOGICAL active(ngridmx,nlayermx) |
---|
| 47 | INTEGER lmax_down(ngridmx),lmin_down(ngridmx) |
---|
[161] | 48 | INTEGER ncorec |
---|
| 49 | |
---|
| 50 | ! =========== Init ============================================== |
---|
| 51 | |
---|
| 52 | entrd(:,:)=0. |
---|
| 53 | detrd(:,:)=0. |
---|
| 54 | qa(:,:)=q_therm(:,:) |
---|
| 55 | q(:,:)=q_therm(:,:) |
---|
| 56 | qd(:,:)=q_therm(:,:) |
---|
| 57 | active(:,:)=.false. |
---|
| 58 | |
---|
| 59 | ! previous calculation of zdthl_down uses the divergence of fmd |
---|
| 60 | ! so it can be negative without problem. Here we include the sign |
---|
| 61 | ! of fmd in the equations, so it has to be positive |
---|
| 62 | ! |
---|
[256] | 63 | fmd(:,:)=-fm_down(:,:) |
---|
[161] | 64 | ! |
---|
| 65 | !! ========== Entrainment, Detrainement and Mass ================= |
---|
| 66 | ! |
---|
| 67 | !! ========== DOWNDRAFT TRANSPORT DISABLED FOR NOW =============== |
---|
| 68 | ! |
---|
[185] | 69 | ! do ig=1,ngridmx |
---|
| 70 | ! if (ztv(ig,nlayermx)-ztvd(ig,nlayermx) .gt. 0.5) then |
---|
[161] | 71 | ! print*,"downdraft non nul derniere couche !!! (dqupdown)" |
---|
| 72 | ! endif |
---|
[185] | 73 | ! detrd(ig,nlayermx)=0. |
---|
| 74 | ! entrd(ig,nlayermx)=0. |
---|
[161] | 75 | ! enddo |
---|
| 76 | ! |
---|
[185] | 77 | ! do k=nlayermx-1,1,-1 |
---|
| 78 | ! do ig=1,ngridmx |
---|
[161] | 79 | ! |
---|
| 80 | ! if (ztv(ig,k)-ztvd(ig,k) .gt. 0.0001) then |
---|
| 81 | ! detrd(ig,k)=MAX(0.,(fmd(ig,k+1)*(ztv(ig,k)-ztvd(ig,k+1))) & |
---|
| 82 | ! & /(ztv(ig,k)-ztvd(ig,k)) - fmd(ig,k)) |
---|
| 83 | ! entrd(ig,k)=MAX(0.,(fmd(ig,k+1)*(ztvd(ig,k)-ztvd(ig,k+1))) & |
---|
| 84 | ! & /(ztv(ig,k)-ztvd(ig,k))) |
---|
| 85 | ! |
---|
| 86 | ! endif |
---|
| 87 | ! enddo |
---|
| 88 | ! enddo |
---|
| 89 | ! |
---|
| 90 | !! ======= We have computed entrainment and detrainment from a prescribed |
---|
| 91 | !! mass flux and potential temp profile. Due to the way downdraft are parametrized, |
---|
| 92 | !! this can yield negative entr and detr. We force it to be positive, but in |
---|
| 93 | !! order to conserve tracers, we need to recompute an adequate mass flux |
---|
| 94 | !! and modify interface rates, to preserve consistency. |
---|
| 95 | ! lmax_down(:)=1 |
---|
| 96 | ! lmin_down(:)=1 |
---|
| 97 | ! |
---|
[185] | 98 | ! do k=1,nlayermx |
---|
| 99 | ! do ig=1,ngridmx |
---|
[161] | 100 | ! if ((entrd(ig,k).gt.0.) .or. (detrd(ig,k).gt.0.)) then |
---|
| 101 | !! if (entrd(ig,k).gt.detrd(ig,k)) then |
---|
| 102 | ! lmax_down(ig)=min(k,lmax(ig)) |
---|
| 103 | ! endif |
---|
| 104 | ! enddo |
---|
| 105 | ! enddo |
---|
[185] | 106 | ! do k=nlayermx,1,-1 |
---|
| 107 | ! do ig=1,ngridmx |
---|
[161] | 108 | ! if ((entrd(ig,k).gt.0.) .or. (detrd(ig,k).gt.0.)) then |
---|
| 109 | !! if (detrd(ig,k).gt.entrd(ig,k)) then |
---|
| 110 | ! lmin_down(ig)=k |
---|
| 111 | ! endif |
---|
| 112 | ! enddo |
---|
| 113 | ! enddo |
---|
| 114 | ! |
---|
| 115 | ! |
---|
| 116 | ! fmd(:,:)=0. |
---|
| 117 | ! |
---|
[185] | 118 | ! do ig=1,ngridmx |
---|
[161] | 119 | ! if ((lmax_down(ig).gt.1) .and. ((lmax_down(ig)-lmin_down(ig)).gt.1)) then |
---|
| 120 | !! fmd(ig,lmax_down(ig))=0. |
---|
| 121 | !! entrd(ig,lmax_down(ig))=detrd(ig,lmax_down(ig)) |
---|
| 122 | !! detrd(ig,lmax_down(ig))=0. |
---|
| 123 | !! print*,lmin_down(ig),lmax_down(ig),lmax(ig) |
---|
| 124 | ! |
---|
| 125 | ! fmd(ig,lmax_down(ig)+1)=0. |
---|
| 126 | ! |
---|
| 127 | ! do k=lmax_down(ig),lmin_down(ig)+1,-1 |
---|
| 128 | ! fmd(ig,k)=fmd(ig,k+1)+entrd(ig,k)-detrd(ig,k) |
---|
| 129 | ! enddo |
---|
| 130 | ! |
---|
| 131 | ! fmd(ig,lmin_down(ig))=0. |
---|
| 132 | ! detrd(ig,lmin_down(ig))=fmd(ig,lmin_down(ig)+1)+entrd(ig,lmin_down(ig)) |
---|
| 133 | ! |
---|
| 134 | ! else |
---|
| 135 | ! entrd(ig,:)=0. |
---|
| 136 | ! detrd(ig,:)=0. |
---|
| 137 | ! active(ig,:)=.false. |
---|
| 138 | ! endif |
---|
| 139 | ! |
---|
| 140 | ! enddo |
---|
| 141 | ! ncorec=0 |
---|
[185] | 142 | ! do k=nlayermx,2,-1 |
---|
| 143 | ! do ig=1,ngridmx |
---|
[161] | 144 | ! if (fmd(ig,k).lt.0.) then |
---|
| 145 | !! detrd(ig,k)=max(0.,detrd(ig,k)+fmd(ig,k-1)) |
---|
| 146 | !! fmd(ig,k-1)=0. |
---|
| 147 | !! entrd(ig,k-1)=0. |
---|
| 148 | !! detrd(ig,k-1)=0. |
---|
| 149 | !! lmin_down(ig)=k-1 |
---|
| 150 | ! fmd(ig,k)=fmd(ig,k+1) |
---|
| 151 | ! detrd(ig,k)=entrd(ig,k) |
---|
| 152 | ! ncorec=ncorec+1 |
---|
| 153 | !! fmd(ig,k)=0. |
---|
| 154 | !! detrd(ig,k)=entrd(ig,k)+fmd(ig,k+1) |
---|
| 155 | ! |
---|
| 156 | ! endif |
---|
| 157 | ! enddo |
---|
| 158 | ! enddo |
---|
| 159 | ! |
---|
| 160 | ! if (ncorec .ne. 0) then |
---|
| 161 | ! print*, 'corrections for negative downward mass flux :',ncorec |
---|
| 162 | ! endif |
---|
| 163 | ! print*, lmin_down(:),lmax_down(:) |
---|
| 164 | ! |
---|
[185] | 165 | ! do k=2,nlayermx |
---|
| 166 | ! do ig=1,ngridmx |
---|
[161] | 167 | ! active(ig,k)=(k.ge.lmin_down(ig)).and.(k.le.lmax_down(ig)) & |
---|
| 168 | ! & .and.(((fmd(ig,k)+detrd(ig,k))*ptimestep).gt.1.e-6*masse0(ig,k)) |
---|
| 169 | ! enddo |
---|
| 170 | ! enddo |
---|
| 171 | ! |
---|
| 172 | ! |
---|
[185] | 173 | ! do ig=1,ngridmx |
---|
[161] | 174 | ! do k=lmin_down(ig),lmax_down(ig) |
---|
| 175 | ! if(.not.active(ig,k)) then |
---|
| 176 | ! active(ig,:)=.false. |
---|
| 177 | ! endif |
---|
| 178 | ! enddo |
---|
| 179 | ! enddo |
---|
| 180 | ! |
---|
| 181 | ! if(charvar .eq. 'tke') then |
---|
| 182 | ! active(:,:)=.false. |
---|
| 183 | ! endif |
---|
| 184 | ! |
---|
[185] | 185 | !! do ig=1,ngridmx |
---|
[161] | 186 | !! active(ig,lmax_down(ig))=(((fmd(ig,lmax_down(ig))+detrd(ig,lmax_down(ig)))* & |
---|
| 187 | !! & ptimestep).gt.1.e-8*masse0(ig,lmax_down(ig))) |
---|
| 188 | !! enddo |
---|
| 189 | !! |
---|
[185] | 190 | !! do ig=1,ngridmx |
---|
[161] | 191 | !! if (lmax_down(ig).gt.1) then |
---|
| 192 | !! do k=lmax_down(ig)-1,lmin_down(ig),-1 |
---|
| 193 | !! active(ig,k)=(((fmd(ig,k)+detrd(ig,k))* & |
---|
| 194 | !! & ptimestep).gt.1.e-8*masse0(ig,k)) & |
---|
| 195 | !! & .and. active(ig,k+1) |
---|
| 196 | !! enddo |
---|
| 197 | !! else |
---|
| 198 | !! active(ig,:)=.false. |
---|
| 199 | !! endif |
---|
| 200 | !! enddo |
---|
| 201 | !! |
---|
| 202 | !! ========== qa : q in updraft ================================== |
---|
| 203 | ! |
---|
[185] | 204 | do k=2,nlayermx |
---|
| 205 | do ig=1,ngridmx |
---|
[161] | 206 | if ((fm0(ig,k+1)+detr0(ig,k))*ptimestep.gt. & |
---|
| 207 | & 1.e-5*masse0(ig,k)) then |
---|
| 208 | qa(ig,k)=(fm0(ig,k)*qa(ig,k-1)+entr0(ig,k)*q(ig,k)) & |
---|
| 209 | & /(fm0(ig,k+1)+detr0(ig,k)) |
---|
| 210 | |
---|
| 211 | if ((qa(ig,k).lt.0.) .and. (charvar .ne. 'momentum')) then |
---|
[278] | 212 | if ((q(ig,k) .gt. 0.) .and. ( q(ig,k) .gt. 1.e-15 )) then |
---|
| 213 | !! only print if it is the thermal scheme which makes qa<0 |
---|
| 214 | print*,'qa<0 created by thermals!!!',charvar,ig,k,fm0(ig,k),qa(ig,k-1),entr0(ig,k),q(ig,k),fm0(ig,k+1),detr0(ig,k) |
---|
| 215 | print*,'---------> Cancelling qa' |
---|
| 216 | endif |
---|
| 217 | qa(ig,k)=q(ig,k) !! no action of thermals in this case. |
---|
[161] | 218 | endif |
---|
| 219 | else |
---|
| 220 | qa(ig,k)=q(ig,k) |
---|
| 221 | endif |
---|
| 222 | enddo |
---|
| 223 | enddo |
---|
| 224 | |
---|
| 225 | ! ========== qd : q in downdraft ================================= |
---|
| 226 | ! |
---|
[185] | 227 | ! do k=nlayermx-1,1,-1 |
---|
| 228 | ! do ig=1,ngridmx |
---|
[161] | 229 | ! if (active(ig,k)) then |
---|
| 230 | ! qd(ig,k)=(fmd(ig,k+1)*qd(ig,k+1)+entrd(ig,k)*q(ig,k)) & |
---|
| 231 | ! & /(fmd(ig,k)+detrd(ig,k)) |
---|
| 232 | ! |
---|
| 233 | ! if ((qd(ig,k).lt.0.) .and. (charvar .ne. 'momentum')) then |
---|
| 234 | ! print*,'qd<0!!!',charvar,ig,k,fmd(ig,k),qd(ig,k),entrd(ig,k),q(ig,k),fmd(ig,k+1),detrd(ig,k),lmin_down(ig),lmax_down(ig) |
---|
| 235 | ! print*, '---------> cancelling qd, no downdraft for this gridpoint' |
---|
| 236 | ! qd(ig,k)=q(ig,k) |
---|
| 237 | ! active(ig,:)=.false. |
---|
| 238 | ! endif |
---|
| 239 | ! else |
---|
| 240 | ! qd(ig,k)=q(ig,k) |
---|
| 241 | ! endif |
---|
| 242 | !! print*,'active,k,entr,detr,q,qd (down) :',active(ig,k),k,entrd(ig,k),detrd(ig,k),q(ig,k),qd(ig,k) |
---|
| 243 | ! enddo |
---|
| 244 | ! enddo |
---|
| 245 | ! |
---|
| 246 | ! |
---|
| 247 | ! ====== dq ====================================================== |
---|
| 248 | |
---|
[185] | 249 | do ig=1,ngridmx |
---|
[161] | 250 | if(active(ig,1)) then |
---|
| 251 | |
---|
| 252 | dq_therm(ig,1)=(detr0(ig,1)*qa(ig,1)+detrd(ig,1)*qd(ig,1) & |
---|
| 253 | & +fm0(ig,2)*q(ig,2) & |
---|
| 254 | & -entr0(ig,1)*q(ig,1)-entrd(ig,1)*q(ig,1) & |
---|
| 255 | & -fmd(ig,2)*q(ig,1)) & |
---|
| 256 | & *ptimestep/masse0(ig,1) |
---|
| 257 | |
---|
| 258 | else |
---|
| 259 | dq_therm(ig,1)=(detr0(ig,1)*qa(ig,1)+fm0(ig,2)*q(ig,2) & |
---|
| 260 | & -entr0(ig,1)*q(ig,1)) & |
---|
| 261 | & *ptimestep/masse0(ig,1) |
---|
| 262 | |
---|
| 263 | endif |
---|
| 264 | enddo |
---|
| 265 | |
---|
[185] | 266 | do k=2,nlayermx-1 |
---|
| 267 | do ig=1, ngridmx |
---|
[161] | 268 | |
---|
| 269 | if(active(ig,k)) then |
---|
| 270 | |
---|
| 271 | dq_therm(ig,k)=(detr0(ig,k)*qa(ig,k)+detrd(ig,k)*qd(ig,k) & |
---|
| 272 | & +fm0(ig,k+1)*q(ig,k+1)+fmd(ig,k)*q(ig,k-1) & |
---|
| 273 | & -entr0(ig,k)*q(ig,k)-entrd(ig,k)*q(ig,k) & |
---|
| 274 | & -fm0(ig,k)*q(ig,k)-fmd(ig,k+1)*q(ig,k)) & |
---|
| 275 | & *ptimestep/masse0(ig,k) |
---|
| 276 | |
---|
| 277 | |
---|
| 278 | else |
---|
| 279 | dq_therm(ig,k)=(detr0(ig,k)*qa(ig,k)+fm0(ig,k+1)*q(ig,k+1) & |
---|
| 280 | & -entr0(ig,k)*q(ig,k)-fm0(ig,k)*q(ig,k)) & |
---|
| 281 | & *ptimestep/masse0(ig,k) |
---|
| 282 | |
---|
| 283 | |
---|
| 284 | endif |
---|
| 285 | |
---|
| 286 | enddo |
---|
| 287 | enddo |
---|
| 288 | |
---|
[185] | 289 | do ig=1, ngridmx |
---|
[161] | 290 | |
---|
[185] | 291 | if(active(ig,nlayermx)) then |
---|
[161] | 292 | |
---|
[185] | 293 | dq_therm(ig,nlayermx)=(detr0(ig,nlayermx)*qa(ig,nlayermx)+detrd(ig,nlayermx)*qd(ig,nlayermx) & |
---|
| 294 | & +fmd(ig,nlayermx)*q(ig,nlayermx-1) & |
---|
| 295 | & -entr0(ig,nlayermx)*q(ig,nlayermx)-entrd(ig,nlayermx)*q(ig,nlayermx) & |
---|
| 296 | & -fm0(ig,nlayermx)*q(ig,nlayermx)) & |
---|
| 297 | & *ptimestep/masse0(ig,nlayermx) |
---|
[161] | 298 | |
---|
| 299 | else |
---|
[185] | 300 | dq_therm(ig,nlayermx)=(detr0(ig,nlayermx)*qa(ig,nlayermx) & |
---|
| 301 | & -entr0(ig,nlayermx)*q(ig,nlayermx)-fm0(ig,nlayermx)*q(ig,nlayermx)) & |
---|
| 302 | & *ptimestep/masse0(ig,nlayermx) |
---|
[161] | 303 | endif |
---|
| 304 | |
---|
| 305 | enddo |
---|
| 306 | return |
---|
| 307 | end |
---|