[4590] | 1 | MODULE lmdz_thermcell_dv2 |
---|
| 2 | CONTAINS |
---|
| 3 | |
---|
[5103] | 4 | SUBROUTINE thermcell_dv2(ngrid,nlay,ptimestep,fm,entr,masse & |
---|
[5087] | 5 | ,fraca,larga & |
---|
| 6 | ,u,v,du,dv,ua,va,lev_out) |
---|
[5112] | 7 | USE lmdz_print_control, ONLY: prt_level,lunout |
---|
[5113] | 8 | IMPLICIT NONE |
---|
[878] | 9 | |
---|
| 10 | !======================================================================= |
---|
[5099] | 11 | |
---|
[878] | 12 | ! Calcul du transport verticale dans la couche limite en presence |
---|
| 13 | ! de "thermiques" explicitement representes |
---|
| 14 | ! calcul du dq/dt une fois qu'on connait les ascendances |
---|
[5099] | 15 | |
---|
[1403] | 16 | ! Vectorisation, FH : 2010/03/08 |
---|
[5099] | 17 | |
---|
[878] | 18 | !======================================================================= |
---|
| 19 | |
---|
| 20 | |
---|
[5117] | 21 | INTEGER ngrid,nlay |
---|
[878] | 22 | |
---|
[5117] | 23 | REAL ptimestep |
---|
| 24 | REAL masse(ngrid,nlay),fm(ngrid,nlay+1) |
---|
| 25 | REAL fraca(ngrid,nlay+1) |
---|
| 26 | REAL larga(ngrid) |
---|
| 27 | REAL entr(ngrid,nlay) |
---|
| 28 | REAL u(ngrid,nlay) |
---|
| 29 | REAL ua(ngrid,nlay) |
---|
| 30 | REAL du(ngrid,nlay) |
---|
| 31 | REAL v(ngrid,nlay) |
---|
| 32 | REAL va(ngrid,nlay) |
---|
| 33 | REAL dv(ngrid,nlay) |
---|
| 34 | INTEGER lev_out ! niveau pour les print |
---|
[878] | 35 | |
---|
[5117] | 36 | REAL qa(ngrid,nlay),detr(ngrid,nlay),zf,zf2 |
---|
| 37 | REAL wvd(ngrid,nlay+1),wud(ngrid,nlay+1) |
---|
| 38 | REAL gamma0(ngrid,nlay+1),gamma(ngrid,nlay+1) |
---|
| 39 | REAL ue(ngrid,nlay),ve(ngrid,nlay) |
---|
[1403] | 40 | LOGICAL ltherm(ngrid,nlay) |
---|
[5117] | 41 | REAL dua(ngrid,nlay),dva(ngrid,nlay) |
---|
| 42 | INTEGER iter |
---|
[878] | 43 | |
---|
[5117] | 44 | INTEGER ig,k,nlarga0 |
---|
[878] | 45 | |
---|
[1403] | 46 | !------------------------------------------------------------------------- |
---|
| 47 | |
---|
[878] | 48 | ! calcul du detrainement |
---|
[1403] | 49 | !--------------------------- |
---|
[878] | 50 | |
---|
[5103] | 51 | ! PRINT*,'THERMCELL DV2 OPTIMISE 3' |
---|
[1403] | 52 | |
---|
| 53 | nlarga0=0. |
---|
| 54 | |
---|
[5158] | 55 | DO k=1,nlay |
---|
| 56 | DO ig=1,ngrid |
---|
[878] | 57 | detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k) |
---|
| 58 | enddo |
---|
| 59 | enddo |
---|
| 60 | |
---|
| 61 | ! calcul de la valeur dans les ascendances |
---|
[5158] | 62 | DO ig=1,ngrid |
---|
[878] | 63 | ua(ig,1)=u(ig,1) |
---|
| 64 | va(ig,1)=v(ig,1) |
---|
| 65 | ue(ig,1)=u(ig,1) |
---|
| 66 | ve(ig,1)=v(ig,1) |
---|
| 67 | enddo |
---|
| 68 | |
---|
[1146] | 69 | IF(prt_level>9)WRITE(lunout,*) & |
---|
[5087] | 70 | 'WARNING on initialise gamma(1:ngrid,1)=0.' |
---|
[972] | 71 | gamma(1:ngrid,1)=0. |
---|
[5158] | 72 | DO k=2,nlay |
---|
| 73 | DO ig=1,ngrid |
---|
[1403] | 74 | ltherm(ig,k)=(fm(ig,k+1)+detr(ig,k))*ptimestep > 1.e-5*masse(ig,k) |
---|
[5117] | 75 | IF(ltherm(ig,k).AND.larga(ig)>0.) THEN |
---|
[1403] | 76 | gamma0(ig,k)=masse(ig,k) & |
---|
[5087] | 77 | *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) ) & |
---|
| 78 | *0.5/larga(ig) & |
---|
| 79 | *1. |
---|
[1403] | 80 | else |
---|
| 81 | gamma0(ig,k)=0. |
---|
| 82 | endif |
---|
[5117] | 83 | IF (ltherm(ig,k).AND.larga(ig)<=0.) nlarga0=nlarga0+1 |
---|
[1403] | 84 | enddo |
---|
| 85 | enddo |
---|
| 86 | |
---|
| 87 | gamma(:,:)=0. |
---|
| 88 | |
---|
[5158] | 89 | DO k=2,nlay |
---|
[1403] | 90 | |
---|
[5158] | 91 | DO ig=1,ngrid |
---|
[5117] | 92 | IF (ltherm(ig,k)) THEN |
---|
[1403] | 93 | dua(ig,k)=ua(ig,k-1)-u(ig,k-1) |
---|
| 94 | dva(ig,k)=va(ig,k-1)-v(ig,k-1) |
---|
| 95 | else |
---|
| 96 | ua(ig,k)=u(ig,k) |
---|
| 97 | va(ig,k)=v(ig,k) |
---|
| 98 | ue(ig,k)=u(ig,k) |
---|
| 99 | ve(ig,k)=v(ig,k) |
---|
| 100 | endif |
---|
| 101 | enddo |
---|
| 102 | |
---|
| 103 | |
---|
| 104 | ! Debut des iterations |
---|
| 105 | !---------------------- |
---|
[5103] | 106 | DO iter=1,5 |
---|
[5158] | 107 | DO ig=1,ngrid |
---|
[1403] | 108 | ! Pour memoire : calcul prenant en compte la fraction reelle |
---|
| 109 | ! zf=0.5*(fraca(ig,k)+fraca(ig,k+1)) |
---|
| 110 | ! zf2=1./(1.-zf) |
---|
| 111 | ! Calcul avec fraction infiniement petite |
---|
[878] | 112 | zf=0. |
---|
[1403] | 113 | zf2=1. |
---|
| 114 | |
---|
[5093] | 115 | ! la première fois on multiplie le coefficient de freinage |
---|
[1403] | 116 | ! par le module du vent dans la couche en dessous. |
---|
| 117 | ! Mais pourquoi donc ??? |
---|
[5117] | 118 | IF (ltherm(ig,k)) THEN |
---|
[878] | 119 | ! On choisit une relaxation lineaire. |
---|
[1403] | 120 | ! gamma(ig,k)=gamma0(ig,k) |
---|
[878] | 121 | ! On choisit une relaxation quadratique. |
---|
[1403] | 122 | gamma(ig,k)=gamma0(ig,k)*sqrt(dua(ig,k)**2+dva(ig,k)**2) |
---|
[878] | 123 | ua(ig,k)=(fm(ig,k)*ua(ig,k-1) & |
---|
[5087] | 124 | +(zf2*entr(ig,k)+gamma(ig,k))*u(ig,k)) & |
---|
| 125 | /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2 & |
---|
| 126 | +gamma(ig,k)) |
---|
[878] | 127 | va(ig,k)=(fm(ig,k)*va(ig,k-1) & |
---|
[5087] | 128 | +(zf2*entr(ig,k)+gamma(ig,k))*v(ig,k)) & |
---|
| 129 | /(fm(ig,k+1)+detr(ig,k)+entr(ig,k)*zf*zf2 & |
---|
| 130 | +gamma(ig,k)) |
---|
[5103] | 131 | ! PRINT*,k,ua(ig,k),va(ig,k),u(ig,k),v(ig,k),dua(ig,k),dva(ig,k) |
---|
[1403] | 132 | dua(ig,k)=ua(ig,k)-u(ig,k) |
---|
| 133 | dva(ig,k)=va(ig,k)-v(ig,k) |
---|
[878] | 134 | ue(ig,k)=(u(ig,k)-zf*ua(ig,k))*zf2 |
---|
| 135 | ve(ig,k)=(v(ig,k)-zf*va(ig,k))*zf2 |
---|
[1403] | 136 | endif |
---|
[878] | 137 | enddo |
---|
[1403] | 138 | ! Fin des iterations |
---|
| 139 | !-------------------- |
---|
[5103] | 140 | END DO |
---|
[878] | 141 | |
---|
[1403] | 142 | enddo ! k=2,nlay |
---|
| 143 | |
---|
| 144 | |
---|
| 145 | ! Calcul du flux vertical de moment dans l'environnement. |
---|
| 146 | !--------------------------------------------------------- |
---|
[5158] | 147 | DO k=2,nlay |
---|
| 148 | DO ig=1,ngrid |
---|
[878] | 149 | wud(ig,k)=fm(ig,k)*ue(ig,k) |
---|
| 150 | wvd(ig,k)=fm(ig,k)*ve(ig,k) |
---|
| 151 | enddo |
---|
| 152 | enddo |
---|
[5158] | 153 | DO ig=1,ngrid |
---|
[878] | 154 | wud(ig,1)=0. |
---|
| 155 | wud(ig,nlay+1)=0. |
---|
| 156 | wvd(ig,1)=0. |
---|
| 157 | wvd(ig,nlay+1)=0. |
---|
| 158 | enddo |
---|
| 159 | |
---|
[1403] | 160 | ! calcul des tendances. |
---|
| 161 | !----------------------- |
---|
[5158] | 162 | DO k=1,nlay |
---|
| 163 | DO ig=1,ngrid |
---|
[878] | 164 | du(ig,k)=((detr(ig,k)+gamma(ig,k))*ua(ig,k) & |
---|
[5087] | 165 | -(entr(ig,k)+gamma(ig,k))*ue(ig,k) & |
---|
| 166 | -wud(ig,k)+wud(ig,k+1)) & |
---|
| 167 | /masse(ig,k) |
---|
[878] | 168 | dv(ig,k)=((detr(ig,k)+gamma(ig,k))*va(ig,k) & |
---|
[5087] | 169 | -(entr(ig,k)+gamma(ig,k))*ve(ig,k) & |
---|
| 170 | -wvd(ig,k)+wvd(ig,k+1)) & |
---|
| 171 | /masse(ig,k) |
---|
[878] | 172 | enddo |
---|
| 173 | enddo |
---|
| 174 | |
---|
[1403] | 175 | |
---|
| 176 | ! Sorties eventuelles. |
---|
| 177 | !---------------------- |
---|
| 178 | |
---|
[5116] | 179 | IF(prt_level>=10) THEN |
---|
[5158] | 180 | DO k=1,nlay |
---|
| 181 | DO ig=1,ngrid |
---|
[5103] | 182 | PRINT*,'th_dv2 ig k gamma entr detr ua ue va ve wud wvd masse',ig,k,gamma(ig,k), & |
---|
[5087] | 183 | entr(ig,k),detr(ig,k),ua(ig,k),ue(ig,k),va(ig,k),ve(ig,k),wud(ig,k),wvd(ig,k),wud(ig,k+1),wvd(ig,k+1), & |
---|
| 184 | masse(ig,k) |
---|
[1403] | 185 | enddo |
---|
| 186 | enddo |
---|
| 187 | endif |
---|
[5099] | 188 | |
---|
[5117] | 189 | IF (nlarga0>0) THEN |
---|
[5103] | 190 | PRINT*,'WARNING !!!!!! DANS THERMCELL_DV2 ' |
---|
| 191 | PRINT*,nlarga0,' points pour lesquels laraga=0. dans un thermique' |
---|
| 192 | PRINT*,'Il faudrait decortiquer ces points' |
---|
[1403] | 193 | endif |
---|
| 194 | |
---|
[5116] | 195 | RETURN |
---|
[5119] | 196 | END |
---|
[4590] | 197 | END MODULE lmdz_thermcell_dv2 |
---|