| 1 | ! |
|---|
| 2 | ! $Id $ |
|---|
| 3 | ! |
|---|
| 4 | SUBROUTINE cvltr_spl(pdtime, da, phi,phi2,d1a,dam, mpIN,epIN, & |
|---|
| 5 | sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm, & |
|---|
| 6 | pmflxrIN,pmflxsIN,ev,te,wdtrainA,wdtrainM, & |
|---|
| 7 | paprs,it,tr,upd,dnd,inb,icb, & |
|---|
| 8 | kk,henry,zrho, ccntrAA_spla,ccntrENV_spla,coefcoli_spla, & |
|---|
| 9 | id_prec,id_fine,id_coss, id_codu, id_scdu, & |
|---|
| 10 | dtrcv,trsptd,dtrSscav,dtrsat,dtrUscav,qDi,qPr, & |
|---|
| 11 | qPa,qMel,qTrdi,dtrcvMA,Mint, & |
|---|
| 12 | zmfd1a,zmfphi2,zmfdam) |
|---|
| 13 | USE IOIPSL |
|---|
| 14 | USE dimphy |
|---|
| 15 | USE infotrac_phy, ONLY : nbtr,tname |
|---|
| 16 | IMPLICIT NONE |
|---|
| 17 | !===================================================================== |
|---|
| 18 | ! Objet : convection des traceurs / KE |
|---|
| 19 | ! Auteurs: M-A Filiberti and J-Y Grandpeix |
|---|
| 20 | ! modifiee par R Pilon : lessivage des traceurs / KE |
|---|
| 21 | !===================================================================== |
|---|
| 22 | |
|---|
| 23 | include "YOMCST.h" |
|---|
| 24 | include "YOECUMF.h" |
|---|
| 25 | include "conema3.h" |
|---|
| 26 | include "chem.h" |
|---|
| 27 | |
|---|
| 28 | ! Entree |
|---|
| 29 | REAL,INTENT(IN) :: pdtime |
|---|
| 30 | REAL,DIMENSION(klon,klev),INTENT(IN) :: da |
|---|
| 31 | REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi |
|---|
| 32 | ! RomP |
|---|
| 33 | REAL,DIMENSION(klon,klev),INTENT(IN) :: d1a,dam ! matrices pour simplifier |
|---|
| 34 | REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2 ! l'ecriture des tendances |
|---|
| 35 | ! |
|---|
| 36 | REAL,DIMENSION(klon,klev),INTENT(IN) :: mpIN |
|---|
| 37 | REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression aux 1/2 couches (bas en haut) |
|---|
| 38 | ! REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression aux 1/2 couches (bas en haut) |
|---|
| 39 | REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr ! q de traceur (bas en haut) |
|---|
| 40 | INTEGER,INTENT(IN) :: it |
|---|
| 41 | REAL,DIMENSION(klon,klev),INTENT(IN) :: upd ! saturated updraft mass flux |
|---|
| 42 | REAL,DIMENSION(klon,klev),INTENT(IN) :: dnd ! saturated downdraft mass flux |
|---|
| 43 | ! |
|---|
| 44 | REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainA ! masses precipitantes de l'asc adiab |
|---|
| 45 | REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainM ! masses precipitantes des melanges |
|---|
| 46 | !JE REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxrIN ! vprecip: eau |
|---|
| 47 | REAL,DIMENSION(klon,klev+1),INTENT(IN) :: pmflxrIN ! vprecip: eau |
|---|
| 48 | !JE REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxsIN ! vprecip: neige |
|---|
| 49 | REAL,DIMENSION(klon,klev+1),INTENT(IN) :: pmflxsIN ! vprecip: neige |
|---|
| 50 | REAL,DIMENSION(klon,klev),INTENT(IN) :: ev ! evaporation cv30_routine |
|---|
| 51 | REAL,DIMENSION(klon,klev),INTENT(IN) :: epIN |
|---|
| 52 | REAL,DIMENSION(klon,klev),INTENT(IN) :: te |
|---|
| 53 | REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij ! fraction dair de lenv |
|---|
| 54 | REAL,DIMENSION(klon,klev),INTENT(IN) :: wght_cvfd ! weights of the layers feeding convection |
|---|
| 55 | REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij ! contenu en eau condensée spécifique/conc deau condensée massique |
|---|
| 56 | REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm ! eau condensee precipitee dans mel masse dair sat |
|---|
| 57 | REAL,DIMENSION(klon,klev),INTENT(IN) :: eplaMm ! eau condensee precipitee dans aa masse dair sat |
|---|
| 58 | |
|---|
| 59 | REAL,DIMENSION(klon,klev),INTENT(IN) :: clw ! contenu en eau condensée dans lasc adiab |
|---|
| 60 | REAL,DIMENSION(klon),INTENT(IN) :: sigd |
|---|
| 61 | INTEGER,DIMENSION(klon),INTENT(IN) :: icb,inb |
|---|
| 62 | ! Sortie |
|---|
| 63 | REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcv ! tendance totale (bas en haut) |
|---|
| 64 | REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcvMA ! M-A Filiberti |
|---|
| 65 | REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: trsptd ! tendance du transport |
|---|
| 66 | REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrSscav ! tendance du lessivage courant sat |
|---|
| 67 | REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrsat ! tendance trsp+sat scav |
|---|
| 68 | REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrUscav ! tendance du lessivage courant unsat |
|---|
| 69 | ! |
|---|
| 70 | ! Variables locales |
|---|
| 71 | INTEGER :: i,j,k |
|---|
| 72 | REAL,DIMENSION(klon,klev) :: dxpres ! difference de pression entre niveau (j+1) et (j) |
|---|
| 73 | REAL :: pdtimeRG ! pas de temps * gravite |
|---|
| 74 | REAL,DIMENSION(klon,nbtr) :: qfeed ! tracer concentration feeding convection |
|---|
| 75 | ! variables pour les courants satures |
|---|
| 76 | REAL,DIMENSION(klon,klev,klev) :: zmd |
|---|
| 77 | REAL,DIMENSION(klon,klev,klev) :: za |
|---|
| 78 | REAL,DIMENSION(klon,klev,nbtr) :: zmfd,zmfa |
|---|
| 79 | REAL,DIMENSION(klon,klev,nbtr) :: zmfp,zmfu |
|---|
| 80 | |
|---|
| 81 | REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfd1a |
|---|
| 82 | REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfdam |
|---|
| 83 | REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfphi2 |
|---|
| 84 | |
|---|
| 85 | ! RomP ! les variables sont nettoyees des valeurs aberrantes |
|---|
| 86 | REAL,DIMENSION(klon,klev) :: Pa, Pm ! pluie AA et mélanges, var temporaire |
|---|
| 87 | REAL,DIMENSION(klon,klev) :: pmflxs,pmflxr ! pmflxrIN,pmflxsIN sans valeur aberante |
|---|
| 88 | REAL,DIMENSION(klon,klev) :: mp ! flux de masse |
|---|
| 89 | REAL,DIMENSION(klon,klev) :: ep ! fraction d'eau convertie en precipitation |
|---|
| 90 | REAL,DIMENSION(klon,klev) :: evap ! evaporation : variable temporaire |
|---|
| 91 | REAL,DIMENSION(klon,klev) :: rho !environmental density |
|---|
| 92 | |
|---|
| 93 | REAL,DIMENSION(klon,klev) :: kappa ! denominateur du au calcul de la matrice |
|---|
| 94 | ! pour obtenir qd et qp |
|---|
| 95 | REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qTrdi ! traceurs descente air insature transport MA |
|---|
| 96 | REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qDi ! traceurs descente insaturees |
|---|
| 97 | REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qPr ! traceurs colonne precipitante |
|---|
| 98 | REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qPa ! traceurs dans les precip issues lasc. adiab. |
|---|
| 99 | REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qMel ! traceurs dans les precip issues des melanges |
|---|
| 100 | REAL,DIMENSION(klon,klev,nbtr) :: qMeltmp ! variable temporaire |
|---|
| 101 | REAL,DIMENSION(klon,klev,nbtr) :: qpmMint |
|---|
| 102 | REAL,DIMENSION(klon,klev),INTENT(OUT) :: Mint |
|---|
| 103 | ! tendances |
|---|
| 104 | REAL :: tdcvMA ! terme de transport de traceur (schema Marie Angele) |
|---|
| 105 | REAL :: trsptrac ! terme de transport de traceur par l'air |
|---|
| 106 | REAL :: scavtrac ! terme de lessivage courant sature |
|---|
| 107 | REAL :: uscavtrac ! terme de lessivage courant insature |
|---|
| 108 | ! impaction |
|---|
| 109 | !!! Correction apres discussion Romain P. / Olivier B. |
|---|
| 110 | !!! REAL,PARAMETER :: rdrop=2.5e-3 ! rayon des gouttes d'eau |
|---|
| 111 | REAL,PARAMETER :: rdrop=1.e-3 ! rayon des gouttes d'eau |
|---|
| 112 | !!! |
|---|
| 113 | REAL,DIMENSION(klon,klev) :: imp ! coefficient d'impaction |
|---|
| 114 | ! parametres lessivage |
|---|
| 115 | REAL :: ccntrAA_coef ! \alpha_a : fract aerosols de l'AA convertis en CCN |
|---|
| 116 | REAL :: ccntrENV_coef ! \beta_m : fract aerosols de l'env convertis en CCN |
|---|
| 117 | REAL :: coefcoli ! coefficient de collision des gouttes sur les aerosols |
|---|
| 118 | ! |
|---|
| 119 | LOGICAL,DIMENSION(klon,klev) :: NO_precip |
|---|
| 120 | ! LOGICAL :: scavON |
|---|
| 121 | ! var tmp tests |
|---|
| 122 | REAL :: conserv |
|---|
| 123 | real :: conservMA |
|---|
| 124 | ! JE SPLA adaptation |
|---|
| 125 | |
|---|
| 126 | INTEGER :: id_prec,id_fine,id_coss, id_codu, id_scdu |
|---|
| 127 | REAL,DIMENSION(nbtr) :: ccntrAA_spla,ccntrENV_spla,coefcoli_spla |
|---|
| 128 | REAL,DIMENSION(klon,klev) :: ccntrAA_coef3d |
|---|
| 129 | REAL,DIMENSION(klon,klev) :: ccntrENV_coef3d |
|---|
| 130 | REAL,DIMENSION(klon,klev) :: coefcoli3d |
|---|
| 131 | |
|---|
| 132 | REAL,DIMENSION(nbtr) :: henry !--cste de Henry mol/l/atm |
|---|
| 133 | REAL,DIMENSION(nbtr) :: kk !--coefficient de var avec T (K) |
|---|
| 134 | REAL :: henry_t !--constante de Henry a T t (mol/l/atm) |
|---|
| 135 | REAL :: f_a !--rapport de la phase aqueuse a la phase gazeuse |
|---|
| 136 | |
|---|
| 137 | REAL, PARAMETER :: ph=5. |
|---|
| 138 | REAL :: K1, K2 |
|---|
| 139 | REAL,DIMENSION(klon,klev) :: zrho |
|---|
| 140 | REAL, PARAMETER :: qliq=1.e-3 ! CONVECTIVE ONLY |
|---|
| 141 | |
|---|
| 142 | ! Je end 20140105 |
|---|
| 143 | |
|---|
| 144 | !JE20140724<< |
|---|
| 145 | !JE SPLA adapt: |
|---|
| 146 | ! |
|---|
| 147 | ! |
|---|
| 148 | !! coefficient lessivage |
|---|
| 149 | ! ccntrAA_coef = 0. |
|---|
| 150 | ! ccntrENV_coef = 0. |
|---|
| 151 | ! coefcoli = 0. |
|---|
| 152 | ! |
|---|
| 153 | !!$OMP MASTER |
|---|
| 154 | ! call getin('ccntrAA_coef',ccntrAA_coef) |
|---|
| 155 | ! call getin('ccntrENV_coef',ccntrENV_coef) |
|---|
| 156 | ! call getin('coefcoli',coefcoli) |
|---|
| 157 | !!$OMP END MASTER |
|---|
| 158 | !!$OMP BARRIER |
|---|
| 159 | !!! JE |
|---|
| 160 | ! do j=1,klev |
|---|
| 161 | ! do i=1,klon |
|---|
| 162 | ! imp(i,j)=0. |
|---|
| 163 | ! ccntrAA_coef3d(i,j)=ccntrAA_coef |
|---|
| 164 | ! ccntrENV_coef3d(i,j)=ccntrENV_coef |
|---|
| 165 | ! coefcoli3d(i,j)=coefcoli |
|---|
| 166 | ! enddo |
|---|
| 167 | ! enddo |
|---|
| 168 | ! |
|---|
| 169 | ! |
|---|
| 170 | !!! for SPLA |
|---|
| 171 | !!!JEend |
|---|
| 172 | ! print*,'cvltr coef lessivage convectif', ccntrAA_coef,ccntrENV_coef,coefcoli |
|---|
| 173 | ! |
|---|
| 174 | !JE20140724>> |
|---|
| 175 | |
|---|
| 176 | ! scavON=.TRUE. |
|---|
| 177 | ! if(scavON) then |
|---|
| 178 | ! ccntrAA_coef = 1. |
|---|
| 179 | ! ccntrENV_coef = 1. |
|---|
| 180 | ! coefcoli = 1. |
|---|
| 181 | ! else |
|---|
| 182 | ! ccntrAA_coef = 0. |
|---|
| 183 | ! ccntrENV_coef = 0. |
|---|
| 184 | ! coefcoli = 0. |
|---|
| 185 | ! endif |
|---|
| 186 | |
|---|
| 187 | ! ====================================================== |
|---|
| 188 | ! calcul de l'impaction |
|---|
| 189 | ! ====================================================== |
|---|
| 190 | !initialisation |
|---|
| 191 | do j=1,klev |
|---|
| 192 | do i=1,klon |
|---|
| 193 | imp(i,j)=0. |
|---|
| 194 | enddo |
|---|
| 195 | enddo |
|---|
| 196 | !JE init 20140103 |
|---|
| 197 | !! impaction sur la surface de la colonne de la descente insaturee |
|---|
| 198 | !! On prend la moyenne des precip entre le niveau i+1 et i |
|---|
| 199 | !! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l) |
|---|
| 200 | !! 1000kg/m3= densité de l'eau |
|---|
| 201 | !! 0.75e-3 = 3/4 /1000 |
|---|
| 202 | !! Par la suite, I est tout le temps multiplié par sig_d pour avoir l'impaction sur la surface de la maille |
|---|
| 203 | !! on le néglige ici pour simplifier le code |
|---|
| 204 | ! do j=1,klev-1 |
|---|
| 205 | ! do i=1,klon |
|---|
| 206 | ! imp(i,j) = coefcoli*0.75e-3/rdrop *& |
|---|
| 207 | ! 0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j)) |
|---|
| 208 | !! rho(i,j)=pplay(i,j)/(rd*te(i,j)) |
|---|
| 209 | ! enddo |
|---|
| 210 | ! enddo |
|---|
| 211 | !JE end 20140103 |
|---|
| 212 | |
|---|
| 213 | ! |
|---|
| 214 | ! initialisation pour flux de traceurs, td et autre |
|---|
| 215 | trsptrac = 0. |
|---|
| 216 | scavtrac = 0. |
|---|
| 217 | uscavtrac = 0. |
|---|
| 218 | qfeed(:,it) = 0. !RL |
|---|
| 219 | DO j=1,klev |
|---|
| 220 | DO i=1,klon |
|---|
| 221 | zmfd(i,j,it)=0. |
|---|
| 222 | zmfa(i,j,it)=0. |
|---|
| 223 | zmfu(i,j,it)=0. |
|---|
| 224 | zmfp(i,j,it)=0. |
|---|
| 225 | zmfphi2(i,j,it)=0. |
|---|
| 226 | zmfd1a(i,j,it)=0. |
|---|
| 227 | zmfdam(i,j,it)=0. |
|---|
| 228 | qDi(i,j,it)=0. |
|---|
| 229 | qPr(i,j,it)=0. |
|---|
| 230 | qPa(i,j,it)=0. |
|---|
| 231 | qMel(i,j,it)=0. |
|---|
| 232 | qMeltmp(i,j,it)=0. |
|---|
| 233 | qTrdi(i,j,it)=0. |
|---|
| 234 | kappa(i,j)=0. |
|---|
| 235 | trsptd(i,j,it)=0. |
|---|
| 236 | dtrsat(i,j,it)=0. |
|---|
| 237 | dtrSscav(i,j,it)=0. |
|---|
| 238 | dtrUscav(i,j,it)=0. |
|---|
| 239 | dtrcv(i,j,it)=0. |
|---|
| 240 | dtrcvMA(i,j,it)=0. |
|---|
| 241 | evap(i,j)=0. |
|---|
| 242 | dxpres(i,j)=0. |
|---|
| 243 | qpmMint(i,j,it)=0. |
|---|
| 244 | Mint(i,j)=0. |
|---|
| 245 | END DO |
|---|
| 246 | END DO |
|---|
| 247 | |
|---|
| 248 | ! suppression des valeurs très faibles (~1e-320) |
|---|
| 249 | ! multiplication de levaporation pour lavoir par unite de temps |
|---|
| 250 | ! et par unite de surface de la maille |
|---|
| 251 | ! -> cv30_unsat : evap : masse evaporee/s/(m2 de la descente) |
|---|
| 252 | DO j=1,klev |
|---|
| 253 | DO i=1,klon |
|---|
| 254 | if(ev(i,j).lt.1.e-16) then |
|---|
| 255 | evap(i,j)=0. |
|---|
| 256 | else |
|---|
| 257 | evap(i,j)=ev(i,j)*sigd(i) |
|---|
| 258 | endif |
|---|
| 259 | END DO |
|---|
| 260 | END DO |
|---|
| 261 | |
|---|
| 262 | DO j=1,klev |
|---|
| 263 | DO i=1,klon |
|---|
| 264 | if(j.lt.klev) then |
|---|
| 265 | if(epIN(i,j).lt.1.e-32) then |
|---|
| 266 | ep(i,j)=0. |
|---|
| 267 | else |
|---|
| 268 | ep(i,j)=epIN(i,j) |
|---|
| 269 | endif |
|---|
| 270 | else |
|---|
| 271 | ep(i,j)=epmax |
|---|
| 272 | endif |
|---|
| 273 | if(mpIN(i,j).lt.1.e-32) then |
|---|
| 274 | mp(i,j)=0. |
|---|
| 275 | else |
|---|
| 276 | mp(i,j)=mpIN(i,j) |
|---|
| 277 | endif |
|---|
| 278 | if(pmflxsIN(i,j).lt.1.e-32) then |
|---|
| 279 | pmflxs(i,j)=0. |
|---|
| 280 | else |
|---|
| 281 | pmflxs(i,j)=pmflxsIN(i,j) |
|---|
| 282 | endif |
|---|
| 283 | if(pmflxrIN(i,j).lt.1.e-32) then |
|---|
| 284 | pmflxr(i,j)=0. |
|---|
| 285 | else |
|---|
| 286 | pmflxr(i,j)=pmflxrIN(i,j) |
|---|
| 287 | endif |
|---|
| 288 | if(wdtrainA(i,j).lt.1.e-32) then |
|---|
| 289 | Pa(i,j)=0. |
|---|
| 290 | else |
|---|
| 291 | Pa(i,j)=wdtrainA(i,j) |
|---|
| 292 | endif |
|---|
| 293 | if(wdtrainM(i,j).lt.1.e-32) then |
|---|
| 294 | Pm(i,j)=0. |
|---|
| 295 | else |
|---|
| 296 | Pm(i,j)=wdtrainM(i,j) |
|---|
| 297 | endif |
|---|
| 298 | END DO |
|---|
| 299 | END DO |
|---|
| 300 | |
|---|
| 301 | !========================================== |
|---|
| 302 | DO j = klev-1,1,-1 |
|---|
| 303 | DO i = 1,klon |
|---|
| 304 | NO_precip(i,j) = (pmflxr(i,j+1)+pmflxs(i,j+1)).lt.1.e-10& |
|---|
| 305 | .and.Pa(i,j).lt.1.e-10.and.Pm(i,j).lt.1.e-10 |
|---|
| 306 | END DO |
|---|
| 307 | END DO |
|---|
| 308 | !============================================================================== |
|---|
| 309 | ! JE SPLA: Calc of ccntrAA_coef,ccntrENV_coef, coefcoli for SPLA aerosols and |
|---|
| 310 | ! precursors. From SPLA code inscav_spl.F |
|---|
| 311 | !print *,'Using SPLA values for cvltr_spl ccntr and coefcoli params' |
|---|
| 312 | ! |
|---|
| 313 | ! |
|---|
| 314 | !IF (it.EQ.2) THEN !--fine aerosol |
|---|
| 315 | ! DO j=1,klev |
|---|
| 316 | ! DO i=1,klon |
|---|
| 317 | ! ccntrAA_coef3d(i,j)=0.7 |
|---|
| 318 | ! ccntrENV_coef3d(i,j)=0.7 |
|---|
| 319 | ! coefcoli3d(i,j)=0.001 |
|---|
| 320 | ! ENDDO |
|---|
| 321 | ! ENDDO |
|---|
| 322 | !ELSEIF (it.EQ.3) THEN !-- Coarse Sea salt aerosol |
|---|
| 323 | ! DO j=1,klev |
|---|
| 324 | ! DO i=1,klon |
|---|
| 325 | ! ccntrAA_coef3d(i,j)=1.0 |
|---|
| 326 | ! ccntrENV_coef3d(i,j)=1.0 |
|---|
| 327 | ! coefcoli3d(i,j)=0.001 |
|---|
| 328 | ! ENDDO |
|---|
| 329 | ! ENDDO |
|---|
| 330 | ! |
|---|
| 331 | !ELSEIF (it.EQ.4) THEN !--Coarse Dust aerosol |
|---|
| 332 | ! DO j=1,klev |
|---|
| 333 | ! DO i=1,klon |
|---|
| 334 | ! ccntrAA_coef3d(i,j)=0.7 |
|---|
| 335 | ! ccntrENV_coef3d(i,j)=0.7 |
|---|
| 336 | ! coefcoli3d(i,j)=0.001 |
|---|
| 337 | ! |
|---|
| 338 | ! ENDDO |
|---|
| 339 | ! ENDDO |
|---|
| 340 | ! Gas precursor: Henry's law |
|---|
| 341 | |
|---|
| 342 | IF (it .EQ. id_prec) THEN |
|---|
| 343 | DO k=1, klev |
|---|
| 344 | DO i=1, klon |
|---|
| 345 | henry_t=henry(it)*exp(-kk(it)*(1./298.-1./te(i,k))) !--mol/l/atm |
|---|
| 346 | K1=1.2e-2*exp(-2010*(1/298.-1/te(i,k))) |
|---|
| 347 | K2=6.6e-8*exp(-1510*(1/298.-1/te(i,k))) |
|---|
| 348 | henry_t=henry_t*(1. + K1/10.**(-ph) + K1*K2/(10.**(-ph))**2) |
|---|
| 349 | f_a=henry_t/101.325*R*te(i,k)*qliq*zrho(i,k)/rho_water |
|---|
| 350 | ! scav(i,k)=f_a/(1.+f_a) |
|---|
| 351 | ccntrAA_coef3d(i,k)= f_a/(1.+f_a) |
|---|
| 352 | ccntrENV_coef3d(i,k)= f_a/(1.+f_a) |
|---|
| 353 | coefcoli3d(i,k)=0.0 |
|---|
| 354 | ENDDO |
|---|
| 355 | ENDDO |
|---|
| 356 | ! CALL minmaxqfi2(clw,1.e33,-1.e33,'minmax clw') |
|---|
| 357 | ELSE |
|---|
| 358 | DO j=1,klev |
|---|
| 359 | DO i=1,klon |
|---|
| 360 | ccntrAA_coef3d(i,j)=ccntrAA_spla(it) |
|---|
| 361 | ccntrENV_coef3d(i,j)=ccntrENV_spla(it) |
|---|
| 362 | coefcoli3d(i,j)=coefcoli_spla(it) |
|---|
| 363 | ENDDO |
|---|
| 364 | ENDDO |
|---|
| 365 | |
|---|
| 366 | |
|---|
| 367 | ENDIF |
|---|
| 368 | |
|---|
| 369 | ! JE end SPLA modifs in ccn parameters |
|---|
| 370 | !============================================================================== |
|---|
| 371 | |
|---|
| 372 | |
|---|
| 373 | |
|---|
| 374 | |
|---|
| 375 | !JE init 20140103 |
|---|
| 376 | ! impaction sur la surface de la colonne de la descente insaturee |
|---|
| 377 | ! On prend la moyenne des precip entre le niveau i+1 et i |
|---|
| 378 | ! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l) |
|---|
| 379 | ! 1000kg/m3= densité de l'eau |
|---|
| 380 | ! 0.75e-3 = 3/4 /1000 |
|---|
| 381 | ! Par la suite, I est tout le temps multiplié par sig_d pour avoir l'impaction sur la surface de la maille |
|---|
| 382 | ! on le néglige ici pour simplifier le code |
|---|
| 383 | do j=1,klev-1 |
|---|
| 384 | do i=1,klon |
|---|
| 385 | !JE imp(i,j) = coefcoli*0.75e-3/rdrop *& |
|---|
| 386 | !JE 0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j)) |
|---|
| 387 | imp(i,j) = coefcoli3d(i,j)*0.75e-3/rdrop *& |
|---|
| 388 | 0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j)) |
|---|
| 389 | ! rho(i,j)=pplay(i,j)/(rd*te(i,j)) |
|---|
| 390 | enddo |
|---|
| 391 | enddo |
|---|
| 392 | !JE end 20140103 |
|---|
| 393 | |
|---|
| 394 | |
|---|
| 395 | ! ========================================= |
|---|
| 396 | ! calcul des tendances liees au downdraft |
|---|
| 397 | ! ========================================= |
|---|
| 398 | !cdir collapse |
|---|
| 399 | DO k=1,klev |
|---|
| 400 | DO j=1,klev |
|---|
| 401 | DO i=1,klon |
|---|
| 402 | zmd(i,j,k)=0. |
|---|
| 403 | za (i,j,k)=0. |
|---|
| 404 | END DO |
|---|
| 405 | END DO |
|---|
| 406 | END DO |
|---|
| 407 | ! calcul de la matrice d echange |
|---|
| 408 | ! matrice de distribution de la masse entrainee en k |
|---|
| 409 | ! commmentaire RomP : mp > 0 |
|---|
| 410 | DO k=1,klev-1 |
|---|
| 411 | DO i=1,klon |
|---|
| 412 | zmd(i,k,k)=max(0.,mp(i,k)-mp(i,k+1)) ! ~ mk(k) |
|---|
| 413 | END DO |
|---|
| 414 | END DO |
|---|
| 415 | DO k=2,klev |
|---|
| 416 | DO j=k-1,1,-1 |
|---|
| 417 | DO i=1,klon |
|---|
| 418 | if(mp(i,j+1).gt.1.e-10) then |
|---|
| 419 | zmd(i,j,k)=zmd(i,j+1,k)*min(1.,mp(i,j)/mp(i,j+1)) !det ~ mk(j)=mk(j+1)*mp(i,j)/mp(i,j+1) |
|---|
| 420 | ENDif |
|---|
| 421 | END DO |
|---|
| 422 | END DO |
|---|
| 423 | END DO |
|---|
| 424 | DO k=1,klev |
|---|
| 425 | DO j=1,klev-1 |
|---|
| 426 | DO i=1,klon |
|---|
| 427 | za(i,j,k)=max(0.,zmd(i,j+1,k)-zmd(i,j,k)) |
|---|
| 428 | END DO |
|---|
| 429 | END DO |
|---|
| 430 | END DO |
|---|
| 431 | !!!!! quantite de traceur dans la descente d'air insaturee : 4 juin 2012 |
|---|
| 432 | DO k=1,klev |
|---|
| 433 | DO j=1,klev-1 |
|---|
| 434 | DO i=1,klon |
|---|
| 435 | if(mp(i,j+1).gt.1.e-10) then |
|---|
| 436 | qTrdi(i,j+1,it)=qTrdi(i,j+1,it)+(zmd(i,j+1,k)/mp(i,j+1))*tr(i,k,it) |
|---|
| 437 | else |
|---|
| 438 | qTrdi(i,j,it)=0.!tr(i,j,it) |
|---|
| 439 | endif |
|---|
| 440 | ENDDO |
|---|
| 441 | ENDDO |
|---|
| 442 | ENDDO |
|---|
| 443 | !!!!! |
|---|
| 444 | ! |
|---|
| 445 | ! rajout du terme lie a l ascendance induite |
|---|
| 446 | ! |
|---|
| 447 | DO j=2,klev |
|---|
| 448 | DO i=1,klon |
|---|
| 449 | za(i,j,j-1)=za(i,j,j-1)+mp(i,j) |
|---|
| 450 | END DO |
|---|
| 451 | END DO |
|---|
| 452 | ! |
|---|
| 453 | ! tendance courants insatures ! sans lessivage ancien schema |
|---|
| 454 | ! |
|---|
| 455 | DO k=1,klev |
|---|
| 456 | DO j=1,klev |
|---|
| 457 | DO i=1,klon |
|---|
| 458 | zmfd(i,j,it)=zmfd(i,j,it)+za(i,j,k)*(tr(i,k,it)-tr(i,j,it)) |
|---|
| 459 | END DO |
|---|
| 460 | END DO |
|---|
| 461 | END DO |
|---|
| 462 | ! |
|---|
| 463 | ! ========================================= |
|---|
| 464 | ! calcul des tendances liees aux courants satures j <-> z ; k <-> z' |
|---|
| 465 | ! ========================================= |
|---|
| 466 | !RL |
|---|
| 467 | ! Feeding concentrations |
|---|
| 468 | DO j=1,klev |
|---|
| 469 | DO i=1,klon |
|---|
| 470 | qfeed(i,it)=qfeed(i,it)+wght_cvfd(i,j)*tr(i,j,it) |
|---|
| 471 | END DO |
|---|
| 472 | END DO |
|---|
| 473 | !RL |
|---|
| 474 | ! |
|---|
| 475 | DO j=1,klev |
|---|
| 476 | DO i=1,klon |
|---|
| 477 | !RL |
|---|
| 478 | !! zmfa(i,j,it)=da(i,j)*(tr(i,1,it)-tr(i,j,it)) ! da |
|---|
| 479 | zmfa(i,j,it)=da(i,j)*(qfeed(i,it)-tr(i,j,it)) ! da |
|---|
| 480 | !RL |
|---|
| 481 | END DO |
|---|
| 482 | END DO |
|---|
| 483 | ! |
|---|
| 484 | DO k=1,klev |
|---|
| 485 | DO j=1,klev |
|---|
| 486 | DO i=1,klon |
|---|
| 487 | zmfp(i,j,it)=zmfp(i,j,it)+phi(i,j,k)*(tr(i,k,it)-tr(i,j,it)) ! phi |
|---|
| 488 | END DO |
|---|
| 489 | END DO |
|---|
| 490 | END DO |
|---|
| 491 | ! RomP ajout des matrices liees au lessivage |
|---|
| 492 | DO j=1,klev |
|---|
| 493 | DO i=1,klon |
|---|
| 494 | zmfd1a(i,j,it)=d1a(i,j)*tr(i,1,it) ! da1 |
|---|
| 495 | zmfdam(i,j,it)=dam(i,j)*tr(i,1,it) ! dam |
|---|
| 496 | END DO |
|---|
| 497 | END DO |
|---|
| 498 | DO k=1,klev |
|---|
| 499 | DO j=1,klev |
|---|
| 500 | DO i=1,klon |
|---|
| 501 | zmfphi2(i,j,it)=zmfphi2(i,j,it)+phi2(i,j,k)*tr(i,k,it) ! psi |
|---|
| 502 | END DO |
|---|
| 503 | END DO |
|---|
| 504 | END DO |
|---|
| 505 | DO j=1,klev-1 |
|---|
| 506 | DO i=1,klon |
|---|
| 507 | zmfu(i,j,it)=max(0.,upd(i,j+1)+dnd(i,j+1))*(tr(i,j+1,it)-tr(i,j,it)) |
|---|
| 508 | END DO |
|---|
| 509 | END DO |
|---|
| 510 | DO j=2,klev |
|---|
| 511 | DO i=1,klon |
|---|
| 512 | zmfu(i,j,it)=zmfu(i,j,it)+min(0.,upd(i,j)+dnd(i,j))*(tr(i,j,it)-tr(i,j-1,it)) |
|---|
| 513 | END DO |
|---|
| 514 | END DO |
|---|
| 515 | ! =================================================== |
|---|
| 516 | ! calcul des tendances liees aux courants insatures |
|---|
| 517 | ! =================================================== |
|---|
| 518 | ! pression |
|---|
| 519 | DO k=1, klev |
|---|
| 520 | DO i=1, klon |
|---|
| 521 | dxpres(i,k)=paprs(i,k)-paprs(i,k+1) |
|---|
| 522 | ENDDO |
|---|
| 523 | ENDDO |
|---|
| 524 | pdtimeRG=pdtime*RG |
|---|
| 525 | |
|---|
| 526 | ! q_pa et q_pm traceurs issues des courants satures se retrouvant dans les precipitations |
|---|
| 527 | DO j=1,klev |
|---|
| 528 | DO i=1,klon |
|---|
| 529 | if(j.ge.icb(i).and.j.le.inb(i)) then |
|---|
| 530 | if(clw(i,j).gt.1.e-16) then |
|---|
| 531 | ! qPa(i,j,it)=ccntrAA_coef*tr(i,1,it)/clw(i,j) |
|---|
| 532 | qPa(i,j,it)=ccntrAA_coef3d(i,j)*tr(i,1,it)/clw(i,j) |
|---|
| 533 | else |
|---|
| 534 | qPa(i,j,it)=0. |
|---|
| 535 | endif |
|---|
| 536 | endif |
|---|
| 537 | END DO |
|---|
| 538 | END DO |
|---|
| 539 | |
|---|
| 540 | ! calcul de q_pm en 2 parties : |
|---|
| 541 | ! 1) calcul de sa valeur pour un niveau z' donne |
|---|
| 542 | ! 2) integration sur la verticale sur z' |
|---|
| 543 | DO j=1,klev |
|---|
| 544 | DO k=1,j-1 |
|---|
| 545 | DO i=1,klon |
|---|
| 546 | if(k.ge.icb(i).and.k.le.inb(i).and.& |
|---|
| 547 | j.le.inb(i)) then |
|---|
| 548 | if(elij(i,k,j).gt.1.e-16) then |
|---|
| 549 | !JE qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_coef*tr(i,1,it)& |
|---|
| 550 | !JE *(1.-sij(i,k,j)) +ccntrENV_coef& |
|---|
| 551 | !JE *tr(i,k,it)*sij(i,k,j)) / elij(i,k,j) |
|---|
| 552 | qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_coef3d(i,k)*tr(i,1,it)& |
|---|
| 553 | *(1.-sij(i,k,j)) +ccntrENV_coef3d(i,k)& |
|---|
| 554 | *tr(i,k,it)*sij(i,k,j)) / elij(i,k,j) |
|---|
| 555 | else |
|---|
| 556 | qMeltmp(i,j,it)=0. |
|---|
| 557 | endif |
|---|
| 558 | qpmMint(i,j,it)=qpmMint(i,j,it) + qMeltmp(i,j,it)*epmlmMm(i,j,k) |
|---|
| 559 | Mint(i,j)=Mint(i,j) + epmlmMm(i,j,k) |
|---|
| 560 | endif ! end if dans nuage |
|---|
| 561 | END DO |
|---|
| 562 | END DO |
|---|
| 563 | END DO |
|---|
| 564 | |
|---|
| 565 | DO j=1,klev |
|---|
| 566 | DO i=1,klon |
|---|
| 567 | if(Mint(i,j).gt.1.e-16) then |
|---|
| 568 | qMel(i,j,it)=qpmMint(i,j,it)/Mint(i,j) |
|---|
| 569 | else |
|---|
| 570 | qMel(i,j,it)=0. |
|---|
| 571 | endif |
|---|
| 572 | END DO |
|---|
| 573 | END DO |
|---|
| 574 | |
|---|
| 575 | ! calcul de q_d et q_p traceurs de la descente precipitante |
|---|
| 576 | DO j=klev-1,1,-1 |
|---|
| 577 | DO i=1,klon |
|---|
| 578 | if(mp(i,j+1).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then ! detrainement |
|---|
| 579 | kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*& |
|---|
| 580 | (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))& |
|---|
| 581 | + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j))) |
|---|
| 582 | |
|---|
| 583 | elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then! entrainement |
|---|
| 584 | if(j.eq.1) then |
|---|
| 585 | kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*& |
|---|
| 586 | (-mp(i,2)-imp(i,j)/RG*dxpres(i,j))& |
|---|
| 587 | + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j))) |
|---|
| 588 | else |
|---|
| 589 | kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*& |
|---|
| 590 | (-mp(i,j)-imp(i,j)/RG*dxpres(i,j))& |
|---|
| 591 | + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j))) |
|---|
| 592 | endif |
|---|
| 593 | else |
|---|
| 594 | kappa(i,j)=1. |
|---|
| 595 | endif |
|---|
| 596 | ENDDO |
|---|
| 597 | ENDDO |
|---|
| 598 | |
|---|
| 599 | DO j=klev-1,1,-1 |
|---|
| 600 | DO i=1,klon |
|---|
| 601 | if (abs(kappa(i,j)).lt.1.e-25) then !si denominateur nul (il peut y avoir des mp!=0) |
|---|
| 602 | kappa(i,j)=1. |
|---|
| 603 | if(j.eq.1) then |
|---|
| 604 | qDi(i,j,it)=qDi(i,j+1,it) !orig tr(i,j,it) ! mp(1)=0 donc tout vient de la couche supérieure |
|---|
| 605 | elseif(mp(i,j+1).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then |
|---|
| 606 | qDi(i,j,it)=qDi(i,j+1,it) |
|---|
| 607 | elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then! entrainement |
|---|
| 608 | qDi(i,j,it)=(-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))/(-mp(i,j)) |
|---|
| 609 | else ! si mp (i)=0 et mp(j+1)=0 |
|---|
| 610 | qDi(i,j,it)=tr(i,j,it) ! orig 0. |
|---|
| 611 | endif |
|---|
| 612 | |
|---|
| 613 | if(NO_precip(i,j)) then |
|---|
| 614 | qPr(i,j,it)=0. |
|---|
| 615 | else |
|---|
| 616 | qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
|---|
| 617 | Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)& |
|---|
| 618 | +imp(i,j)/RG*dxpres(i,j)*qDi(i,j,it))/& |
|---|
| 619 | (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j)) |
|---|
| 620 | endif |
|---|
| 621 | else ! denominateur non nul |
|---|
| 622 | kappa(i,j)=1./kappa(i,j) |
|---|
| 623 | ! calcul de qd et qp |
|---|
| 624 | !!jyg (20130119) correction pour le sommet du nuage |
|---|
| 625 | !! if(j.ge.inb(i)) then !au-dessus du nuage, sommet inclu |
|---|
| 626 | if(j.gt.inb(i)) then !au-dessus du nuage |
|---|
| 627 | qDi(i,j,it)=tr(i,j,it) ! pas de descente => environnement = descente insaturee |
|---|
| 628 | qPr(i,j,it)=0. |
|---|
| 629 | |
|---|
| 630 | ! vvv premiere couche du modele ou mp(1)=0 ! det tout le temps vvv |
|---|
| 631 | elseif(j.eq.1) then |
|---|
| 632 | if(mp(i,2).gt.1.e-10) then !mp(2) non nul -> detrainement (car mp(1) = 0) !ent pas possible |
|---|
| 633 | if(NO_precip(i,j)) then !pas de precip en (i) |
|---|
| 634 | qDi(i,j,it)=qDi(i,j+1,it) |
|---|
| 635 | qPr(i,j,it)=0. |
|---|
| 636 | else |
|---|
| 637 | qDi(i,j,it)=kappa(i,j)*(& |
|---|
| 638 | (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
|---|
| 639 | Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +& |
|---|
| 640 | (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*& |
|---|
| 641 | (-mp(i,j+1)*qDi(i,j+1,it))) |
|---|
| 642 | |
|---|
| 643 | qPr(i,j,it)=kappa(i,j)*(& |
|---|
| 644 | (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))*& |
|---|
| 645 | ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
|---|
| 646 | Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))& |
|---|
| 647 | +(-mp(i,j+1)*qDi(i,j+1,it)) * (imp(i,j)/RG*dxpres(i,j))) |
|---|
| 648 | endif |
|---|
| 649 | |
|---|
| 650 | else !mp(2) nul -> plus de descente insaturee -> pluie agit sur environnement |
|---|
| 651 | qDi(i,j,it)=tr(i,j,it) ! orig 0. |
|---|
| 652 | if(NO_precip(i,j)) then |
|---|
| 653 | qPr(i,j,it)=0. |
|---|
| 654 | else |
|---|
| 655 | qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
|---|
| 656 | Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)& |
|---|
| 657 | +imp(i,j)/RG*dxpres(i,j)*tr(i,j,it))/& |
|---|
| 658 | (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j)) |
|---|
| 659 | endif |
|---|
| 660 | |
|---|
| 661 | endif !mp(2) nul ou non |
|---|
| 662 | |
|---|
| 663 | ! vvv (j!=1.and.j.lt.inb(i)) en-dessous du sommet nuage vvv |
|---|
| 664 | else |
|---|
| 665 | !------------------------------------------------------------- detrainement |
|---|
| 666 | if(mp(i,j+1).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then !mp(i,j).gt.1.e-10) then |
|---|
| 667 | if(NO_precip(i,j)) then |
|---|
| 668 | qDi(i,j,it)=qDi(i,j+1,it) |
|---|
| 669 | qPr(i,j,it)=0. |
|---|
| 670 | else |
|---|
| 671 | qDi(i,j,it)=kappa(i,j)*(& |
|---|
| 672 | (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
|---|
| 673 | Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +& |
|---|
| 674 | (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*& |
|---|
| 675 | (-mp(i,j+1)*qDi(i,j+1,it))) |
|---|
| 676 | ! |
|---|
| 677 | qPr(i,j,it)=kappa(i,j)*(& |
|---|
| 678 | (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))*& |
|---|
| 679 | ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
|---|
| 680 | Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))& |
|---|
| 681 | +(-mp(i,j+1)*qDi(i,j+1,it)) * (imp(i,j)/RG*dxpres(i,j))) |
|---|
| 682 | endif !precip |
|---|
| 683 | !------------------------------------------------------------- entrainement |
|---|
| 684 | elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then |
|---|
| 685 | if(NO_precip(i,j)) then |
|---|
| 686 | qDi(i,j,it)=(-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))/(-mp(i,j)) |
|---|
| 687 | qPr(i,j,it)=0. |
|---|
| 688 | else |
|---|
| 689 | qDi(i,j,it)=kappa(i,j)*(& |
|---|
| 690 | (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
|---|
| 691 | Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +& |
|---|
| 692 | (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*& |
|---|
| 693 | (-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))) |
|---|
| 694 | ! |
|---|
| 695 | qPr(i,j,it)=kappa(i,j)*(& |
|---|
| 696 | (-mp(i,j)-imp(i,j)/RG*dxpres(i,j))*& |
|---|
| 697 | ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
|---|
| 698 | Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))& |
|---|
| 699 | +(-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))*& |
|---|
| 700 | (imp(i,j)/RG*dxpres(i,j))) |
|---|
| 701 | endif !precip |
|---|
| 702 | !------------------------------------------------------------- endif ! ent/det |
|---|
| 703 | else !mp nul |
|---|
| 704 | qDi(i,j,it)=tr(i,j,it) ! orig 0. |
|---|
| 705 | if(NO_precip(i,j)) then |
|---|
| 706 | qPr(i,j,it)=0. |
|---|
| 707 | else |
|---|
| 708 | qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
|---|
| 709 | Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)& |
|---|
| 710 | +imp(i,j)/RG*dxpres(i,j)*tr(i,j,it))/& |
|---|
| 711 | (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j)) |
|---|
| 712 | endif |
|---|
| 713 | endif ! mp nul ou non |
|---|
| 714 | endif ! condition sur j |
|---|
| 715 | endif ! kappa |
|---|
| 716 | ENDDO |
|---|
| 717 | ENDDO |
|---|
| 718 | |
|---|
| 719 | !! print test descente insaturee |
|---|
| 720 | ! DO j=klev,1,-1 |
|---|
| 721 | ! DO i=1,klon |
|---|
| 722 | ! if(it.eq.3) then |
|---|
| 723 | ! write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') j,& |
|---|
| 724 | !! 'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),& |
|---|
| 725 | ! 'zmfdam+zmfpsi',zmfdam(i,j,it)+zmfphi2(i,j,it),'qpmMint',qpmMint(i,j,it),& |
|---|
| 726 | ! 'Pm',Pm(i,j),'Mint',Mint(i,j),& |
|---|
| 727 | !! 'zmfa',zmfa(i,j,it),'zmfp',zmfp(i,j,it),& |
|---|
| 728 | ! 'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),'zmfd1a',zmfd1a(i,j,it) |
|---|
| 729 | !! 'Pa',Pa(i,j),'eplaMm',eplaMm(i,j) |
|---|
| 730 | !! 'zmfd1a=da1*qa',zmfd1a(i,j,it),'Pa*qPa',wdtrainA(i,j)*qPa(i,j,it),'da1',d1a(i,j) |
|---|
| 731 | ! endif |
|---|
| 732 | ! ENDDO |
|---|
| 733 | ! ENDDO |
|---|
| 734 | |
|---|
| 735 | |
|---|
| 736 | ! =================================================== |
|---|
| 737 | ! calcul final des tendances |
|---|
| 738 | ! =================================================== |
|---|
| 739 | |
|---|
| 740 | DO k=klev-1,1,-1 |
|---|
| 741 | DO i=1, klon |
|---|
| 742 | ! transport |
|---|
| 743 | tdcvMA=zmfd(i,k,it)+zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it) ! double comptage des downdraft insatures |
|---|
| 744 | trsptrac=zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it) |
|---|
| 745 | ! lessivage courants satures |
|---|
| 746 | !JE scavtrac=-ccntrAA_coef*zmfd1a(i,k,it)& |
|---|
| 747 | !JE -zmfphi2(i,k,it)*ccntrENV_coef& |
|---|
| 748 | !JE -zmfdam(i,k,it)*ccntrAA_coef |
|---|
| 749 | scavtrac=-ccntrAA_coef3d(i,k)*zmfd1a(i,k,it)& |
|---|
| 750 | -zmfphi2(i,k,it)*ccntrENV_coef3d(i,k)& |
|---|
| 751 | -zmfdam(i,k,it)*ccntrAA_coef3d(i,k) |
|---|
| 752 | ! lessivage courants insatures |
|---|
| 753 | if(k.le.inb(i).and.k.gt.1) then ! tendances dans le nuage |
|---|
| 754 | !------------------------------------------------------------- detrainement |
|---|
| 755 | if(mp(i,k+1).gt.mp(i,k).and.mp(i,k+1).gt.1.e-10) then |
|---|
| 756 | uscavtrac= (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))& |
|---|
| 757 | + mp(i,k)*(tr(i,k-1,it)-tr(i,k,it)) |
|---|
| 758 | ! |
|---|
| 759 | ! if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' det incloud',& |
|---|
| 760 | ! (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+& |
|---|
| 761 | ! mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k),& |
|---|
| 762 | ! 'mp',mp(i,k) |
|---|
| 763 | !------------------------------------------------------------- entrainement |
|---|
| 764 | elseif(mp(i,k).gt.mp(i,k+1).and.mp(i,k).gt.1.e-10) then |
|---|
| 765 | uscavtrac= mp(i,k)*(tr(i,k-1,it)-tr(i,k,it)) |
|---|
| 766 | ! |
|---|
| 767 | ! if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' ent incloud',uscavtrac*pdtimeRG/dxpres(i,k), 'mp',mp(i,k) |
|---|
| 768 | !=!------------------------------------------------------------- end ent/det |
|---|
| 769 | else ! mp(i,k+1)=0. et mp(i,k)=0. pluie directement sur l environnement |
|---|
| 770 | |
|---|
| 771 | if(NO_precip(i,k)) then |
|---|
| 772 | uscavtrac=0. |
|---|
| 773 | ! if(it.eq.3) write(*,'(I2,1X,a,e20.12,82X,a,e20.12)')k,' no P ent incloud',uscavtrac*pdtimeRG/dxpres(i,k), 'mp',mp(i,k) |
|---|
| 774 | else |
|---|
| 775 | uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG |
|---|
| 776 | ! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac*pdtimeRG/dxpres(i,k), 'mp',mp(i,k) |
|---|
| 777 | !!JE adds |
|---|
| 778 | ! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'imp',imp(i,k) |
|---|
| 779 | ! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'tr',tr(i,k,it) |
|---|
| 780 | ! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'evap',evap(i,k) |
|---|
| 781 | ! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'qPr',qPr(i,k,it) |
|---|
| 782 | ! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'dxpres',dxpres(i,k) |
|---|
| 783 | !!Je end |
|---|
| 784 | |
|---|
| 785 | endif |
|---|
| 786 | endif ! mp/det/ent |
|---|
| 787 | !------------------------------------------------------------- premiere couche |
|---|
| 788 | elseif(k.eq.1) then ! mp(1)=0. |
|---|
| 789 | if(mp(i,2).gt.1.e-10) then !detrainement |
|---|
| 790 | uscavtrac= (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it)) !& |
|---|
| 791 | ! + mp(i,2)*(0.-tr(i,k,it)) |
|---|
| 792 | ! |
|---|
| 793 | ! if(it.eq.3) write(*,'(I2,1X,a,e20.12,84X,a,e20.12)')k,' 1 det',& |
|---|
| 794 | ! (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+& |
|---|
| 795 | ! mp(i,2)*(0.-tr(i,k,it))*pdtimeRG/dxpres(i,k),& |
|---|
| 796 | ! 'mp',mp(i,k) |
|---|
| 797 | else ! mp(2) = 0 = mp(1) pas de descente insaturee, rien ne se passe s'il ne pleut pas, sinon pluie->env |
|---|
| 798 | if(NO_precip(i,1)) then |
|---|
| 799 | uscavtrac=0. |
|---|
| 800 | else |
|---|
| 801 | uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG |
|---|
| 802 | endif |
|---|
| 803 | ! if(it.eq.3) write(*,'(I2,1X,a,2X,e20.12,82X,a,e20.12)')k,'1 P env incloud',uscavtrac*pdtimeRG/dxpres(i,k), 'mp',mp(i,k) |
|---|
| 804 | endif |
|---|
| 805 | |
|---|
| 806 | else ! k > INB au-dessus du nuage |
|---|
| 807 | uscavtrac=0. |
|---|
| 808 | endif |
|---|
| 809 | |
|---|
| 810 | ! ===== tendances finales ====== |
|---|
| 811 | trsptd(i,k,it)=trsptrac*pdtimeRG/dxpres(i,k) ! td transport sans eau dans courants satures |
|---|
| 812 | dtrSscav(i,k,it)=scavtrac*pdtimeRG/dxpres(i,k) ! td du lessivage dans courants satures |
|---|
| 813 | dtrUscav(i,k,it)=uscavtrac*pdtimeRG/dxpres(i,k) ! td courant insat |
|---|
| 814 | dtrsat(i,k,it)=(trsptrac+scavtrac)*pdtimeRG/dxpres(i,k) ! td courant sat |
|---|
| 815 | dtrcv(i,k,it)=(trsptrac+scavtrac+uscavtrac)*pdtimeRG/dxpres(i,k)!dtrsat(i,k,it)+dtrUscav(i,k,it) td conv |
|---|
| 816 | !!!!!! |
|---|
| 817 | dtrcvMA(i,k,it)=tdcvMA*pdtimeRG/dxpres(i,k) ! MA tendance convection |
|---|
| 818 | ENDDO |
|---|
| 819 | ENDDO |
|---|
| 820 | |
|---|
| 821 | ! test de conservation du traceur |
|---|
| 822 | !print*,"_____________________________________________________________" |
|---|
| 823 | !print*," " |
|---|
| 824 | ! conserv=0. |
|---|
| 825 | ! conservMA=0. |
|---|
| 826 | ! DO k= klev-1,1,-1 |
|---|
| 827 | ! DO i=1, klon |
|---|
| 828 | ! conserv=conserv+dtrcv(i,k,it)* & |
|---|
| 829 | ! (paprs(i,k)-paprs(i,k+1))/RG |
|---|
| 830 | ! conservMA=conservMA+dtrcvMA(i,k,it)* & |
|---|
| 831 | ! (paprs(i,k)-paprs(i,k+1))/RG |
|---|
| 832 | ! |
|---|
| 833 | ! if(it.eq.3) write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') k,& |
|---|
| 834 | ! 'MA td ',dtrcvMA(i,k,it)*dxpres(i,k)/RG,& |
|---|
| 835 | ! ' td',dtrcv(i,k,it)*dxpres(i,k)/RG,' conservMA ',conservMA,'conserv ',conserv |
|---|
| 836 | !! |
|---|
| 837 | ! ENDDO |
|---|
| 838 | ! ENDDO |
|---|
| 839 | ! if(it.eq.3) print *,'it',it,'conserv ',conserv,'conservMA ',conservMA |
|---|
| 840 | |
|---|
| 841 | END SUBROUTINE cvltr_spl |
|---|
