[1403] | 1 | ! |
---|
| 2 | ! $Id: wake.F 1403 2010-07-01 09:02:53Z lfalletti $ |
---|
| 3 | ! |
---|
| 4 | Subroutine WAKE (p,ph,pi,dtime,sigd_con |
---|
[953] | 5 | : ,te0,qe0,omgb |
---|
[879] | 6 | : ,dtdwn,dqdwn,amdwn,amup,dta,dqa |
---|
| 7 | : ,wdtPBL,wdqPBL,udtPBL,udqPBL |
---|
| 8 | o ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl |
---|
| 9 | o ,dtls,dqls |
---|
| 10 | o ,ktopw,omgbdth,dp_omgb,wdens |
---|
| 11 | o ,tu,qu |
---|
| 12 | o ,dtKE,dqKE |
---|
| 13 | o ,dtPBL,dqPBL |
---|
| 14 | o ,omg,dp_deltomg,spread |
---|
| 15 | o ,Cstar,d_deltat_gw |
---|
| 16 | o ,d_deltatw2,d_deltaqw2) |
---|
| 17 | |
---|
[1146] | 18 | |
---|
[879] | 19 | *************************************************************** |
---|
| 20 | * * |
---|
| 21 | * WAKE * |
---|
| 22 | * retour a un Pupper fixe * |
---|
| 23 | * * |
---|
| 24 | * written by : GRANDPEIX Jean-Yves 09/03/2000 * |
---|
| 25 | * modified by : ROEHRIG Romain 01/29/2007 * |
---|
| 26 | *************************************************************** |
---|
| 27 | c |
---|
[974] | 28 | use dimphy |
---|
| 29 | IMPLICIT none |
---|
| 30 | c============================================================================ |
---|
| 31 | C |
---|
| 32 | C |
---|
| 33 | C But : Decrire le comportement des poches froides apparaissant dans les |
---|
| 34 | C grands systemes convectifs, et fournir l'energie disponible pour |
---|
| 35 | C le declenchement de nouvelles colonnes convectives. |
---|
| 36 | C |
---|
| 37 | C Variables d'etat : deltatw : ecart de temperature wake-undisturbed area |
---|
| 38 | C deltaqw : ecart d'humidite wake-undisturbed area |
---|
| 39 | C sigmaw : fraction d'aire occupee par la poche. |
---|
| 40 | C |
---|
| 41 | C Variable de sortie : |
---|
| 42 | c |
---|
| 43 | c wape : WAke Potential Energy |
---|
| 44 | c fip : Front Incident Power (W/m2) - ALP |
---|
| 45 | c gfl : Gust Front Length per unit area (m-1) |
---|
| 46 | C dtls : large scale temperature tendency due to wake |
---|
| 47 | C dqls : large scale humidity tendency due to wake |
---|
| 48 | C hw : hauteur de la poche |
---|
| 49 | C dp_omgb : vertical gradient of large scale omega |
---|
[1403] | 50 | C wdens : densite de poches |
---|
[974] | 51 | C omgbdth: flux of Delta_Theta transported by LS omega |
---|
| 52 | C dtKE : differential heating (wake - unpertubed) |
---|
| 53 | C dqKE : differential moistening (wake - unpertubed) |
---|
| 54 | C omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s) |
---|
| 55 | C dp_deltomg : vertical gradient of omg (s-1) |
---|
| 56 | C spread : spreading term in dt_wake and dq_wake |
---|
| 57 | C deltatw : updated temperature difference (T_w-T_u). |
---|
| 58 | C deltaqw : updated humidity difference (q_w-q_u). |
---|
| 59 | C sigmaw : updated wake fractional area. |
---|
| 60 | C d_deltat_gw : delta T tendency due to GW |
---|
| 61 | c |
---|
| 62 | C Variables d'entree : |
---|
| 63 | c |
---|
| 64 | c aire : aire de la maille |
---|
| 65 | c te0 : temperature dans l'environnement (K) |
---|
| 66 | C qe0 : humidite dans l'environnement (kg/kg) |
---|
| 67 | C omgb : vitesse verticale moyenne sur la maille (Pa/s) |
---|
| 68 | C dtdwn: source de chaleur due aux descentes (K/s) |
---|
| 69 | C dqdwn: source d'humidite due aux descentes (kg/kg/s) |
---|
| 70 | C dta : source de chaleur due courants satures et detrain (K/s) |
---|
| 71 | C dqa : source d'humidite due aux courants satures et detra (kg/kg/s) |
---|
| 72 | C amdwn: flux de masse total des descentes, par unite de |
---|
| 73 | C surface de la maille (kg/m2/s) |
---|
| 74 | C amup : flux de masse total des ascendances, par unite de |
---|
| 75 | C surface de la maille (kg/m2/s) |
---|
| 76 | C p : pressions aux milieux des couches (Pa) |
---|
| 77 | C ph : pressions aux interfaces (Pa) |
---|
[1403] | 78 | C pi : (p/p_0)**kapa (adim) |
---|
[974] | 79 | C dtime: increment temporel (s) |
---|
| 80 | c |
---|
| 81 | C Variables internes : |
---|
| 82 | c |
---|
| 83 | c rhow : masse volumique de la poche froide |
---|
| 84 | C rho : environment density at P levels |
---|
| 85 | C rhoh : environment density at Ph levels |
---|
| 86 | C te : environment temperature | may change within |
---|
| 87 | C qe : environment humidity | sub-time-stepping |
---|
| 88 | C the : environment potential temperature |
---|
| 89 | C thu : potential temperature in undisturbed area |
---|
| 90 | C tu : temperature in undisturbed area |
---|
| 91 | C qu : humidity in undisturbed area |
---|
| 92 | C dp_omgb: vertical gradient og LS omega |
---|
| 93 | C omgbw : wake average vertical omega |
---|
| 94 | C dp_omgbw: vertical gradient of omgbw |
---|
| 95 | C omgbdq : flux of Delta_q transported by LS omega |
---|
| 96 | C dth : potential temperature diff. wake-undist. |
---|
| 97 | C th1 : first pot. temp. for vertical advection (=thu) |
---|
| 98 | C th2 : second pot. temp. for vertical advection (=thw) |
---|
| 99 | C q1 : first humidity for vertical advection |
---|
| 100 | C q2 : second humidity for vertical advection |
---|
| 101 | C d_deltatw : terme de redistribution pour deltatw |
---|
| 102 | C d_deltaqw : terme de redistribution pour deltaqw |
---|
| 103 | C deltatw0 : deltatw initial |
---|
| 104 | C deltaqw0 : deltaqw initial |
---|
| 105 | C hw0 : hw initial |
---|
| 106 | C sigmaw0: sigmaw initial |
---|
| 107 | C amflux : horizontal mass flux through wake boundary |
---|
[1403] | 108 | C wdens_ref: initial number of wakes per unit area (3D) or per |
---|
| 109 | C unit length (2D), at the beginning of each time step |
---|
[974] | 110 | C Tgw : 1 sur la période de onde de gravité |
---|
| 111 | c Cgw : vitesse de propagation de onde de gravité |
---|
| 112 | c LL : distance entre 2 poches |
---|
| 113 | |
---|
| 114 | c------------------------------------------------------------------------- |
---|
| 115 | c Déclaration de variables |
---|
| 116 | c------------------------------------------------------------------------- |
---|
| 117 | |
---|
| 118 | #include "dimensions.h" |
---|
| 119 | #include "YOMCST.h" |
---|
| 120 | #include "cvthermo.h" |
---|
| 121 | #include "iniprint.h" |
---|
| 122 | |
---|
| 123 | c Arguments en entree |
---|
| 124 | c-------------------- |
---|
| 125 | |
---|
[1403] | 126 | REAL, dimension(klon,klev) :: p, pi |
---|
[974] | 127 | REAL, dimension(klon,klev+1) :: ph, omgb |
---|
| 128 | REAL dtime |
---|
| 129 | REAL, dimension(klon,klev) :: te0,qe0 |
---|
| 130 | REAL, dimension(klon,klev) :: dtdwn, dqdwn |
---|
| 131 | REAL, dimension(klon,klev) :: wdtPBL,wdqPBL |
---|
| 132 | REAL, dimension(klon,klev) :: udtPBL,udqPBL |
---|
| 133 | REAL, dimension(klon,klev) :: amdwn, amup |
---|
| 134 | REAL, dimension(klon,klev) :: dta, dqa |
---|
| 135 | REAL, dimension(klon) :: sigd_con |
---|
| 136 | |
---|
| 137 | c Sorties |
---|
| 138 | c-------- |
---|
| 139 | |
---|
| 140 | REAL, dimension(klon,klev) :: deltatw, deltaqw, dth |
---|
| 141 | REAL, dimension(klon,klev) :: tu, qu |
---|
| 142 | REAL, dimension(klon,klev) :: dtls, dqls |
---|
| 143 | REAL, dimension(klon,klev) :: dtKE, dqKE |
---|
| 144 | REAL, dimension(klon,klev) :: dtPBL, dqPBL |
---|
| 145 | REAL, dimension(klon,klev) :: spread |
---|
| 146 | REAL, dimension(klon,klev) :: d_deltatgw |
---|
| 147 | REAL, dimension(klon,klev) :: d_deltatw2, d_deltaqw2 |
---|
| 148 | REAL, dimension(klon,klev+1) :: omgbdth, omg |
---|
| 149 | REAL, dimension(klon,klev) :: dp_omgb, dp_deltomg |
---|
| 150 | REAL, dimension(klon,klev) :: d_deltat_gw |
---|
| 151 | REAL, dimension(klon) :: hw, sigmaw, wape, fip, gfl, Cstar |
---|
[1403] | 152 | REAL, dimension(klon) :: wdens |
---|
[974] | 153 | INTEGER, dimension(klon) :: ktopw |
---|
| 154 | |
---|
| 155 | c Variables internes |
---|
| 156 | c------------------- |
---|
| 157 | |
---|
| 158 | c Variables à fixer |
---|
| 159 | REAL ALON |
---|
| 160 | REAL coefgw |
---|
[1403] | 161 | REAL :: wdens_ref |
---|
[974] | 162 | REAL stark |
---|
| 163 | REAL alpk |
---|
| 164 | REAL delta_t_min |
---|
| 165 | INTEGER nsub |
---|
| 166 | REAL dtimesub |
---|
[1403] | 167 | REAL sigmad, hwmin,wapecut |
---|
[1146] | 168 | REAL :: sigmaw_max |
---|
[1403] | 169 | REAL :: dens_rate |
---|
| 170 | REAL wdens0 |
---|
[974] | 171 | cIM 080208 |
---|
| 172 | LOGICAL, dimension(klon) :: gwake |
---|
| 173 | |
---|
| 174 | c Variables de sauvegarde |
---|
| 175 | REAL, dimension(klon,klev) :: deltatw0 |
---|
| 176 | REAL, dimension(klon,klev) :: deltaqw0 |
---|
| 177 | REAL, dimension(klon,klev) :: te, qe |
---|
| 178 | REAL, dimension(klon) :: sigmaw0, sigmaw1 |
---|
| 179 | |
---|
| 180 | c Variables pour les GW |
---|
| 181 | REAL, DIMENSION(klon) :: LL |
---|
| 182 | REAL, dimension(klon,klev) :: N2 |
---|
| 183 | REAL, dimension(klon,klev) :: Cgw |
---|
| 184 | REAL, dimension(klon,klev) :: Tgw |
---|
| 185 | |
---|
| 186 | c Variables liées au calcul de hw |
---|
| 187 | REAL, DIMENSION(klon) :: ptop_provis, ptop, ptop_new |
---|
| 188 | REAL, DIMENSION(klon) :: sum_dth |
---|
| 189 | REAL, DIMENSION(klon) :: dthmin |
---|
| 190 | REAL, DIMENSION(klon) :: z, dz, hw0 |
---|
| 191 | INTEGER, DIMENSION(klon) :: ktop, kupper |
---|
| 192 | |
---|
[1146] | 193 | c Sub-timestep tendencies and related variables |
---|
| 194 | REAL d_deltatw(klon,klev),d_deltaqw(klon,klev) |
---|
| 195 | REAL d_te(klon,klev),d_qe(klon,klev) |
---|
| 196 | REAL d_sigmaw(klon),alpha(klon) |
---|
| 197 | REAL q0_min(klon),q1_min(klon) |
---|
| 198 | LOGICAL wk_adv(klon), OK_qx_qw(klon) |
---|
[1403] | 199 | REAL epsilon |
---|
| 200 | DATA epsilon/1.e-15/ |
---|
[1146] | 201 | |
---|
[974] | 202 | c Autres variables internes |
---|
| 203 | INTEGER isubstep, k, i |
---|
| 204 | |
---|
| 205 | REAL, DIMENSION(klon) :: sum_thu, sum_tu, sum_qu,sum_thvu |
---|
| 206 | REAL, DIMENSION(klon) :: sum_dq, sum_rho |
---|
| 207 | REAL, DIMENSION(klon) :: sum_dtdwn, sum_dqdwn |
---|
| 208 | REAL, DIMENSION(klon) :: av_thu, av_tu, av_qu, av_thvu |
---|
| 209 | REAL, DIMENSION(klon) :: av_dth, av_dq, av_rho |
---|
| 210 | REAL, DIMENSION(klon) :: av_dtdwn, av_dqdwn |
---|
| 211 | |
---|
| 212 | REAL, DIMENSION(klon,klev) :: rho, rhow |
---|
| 213 | REAL, DIMENSION(klon,klev+1) :: rhoh |
---|
| 214 | REAL, DIMENSION(klon,klev) :: rhow_moyen |
---|
| 215 | REAL, DIMENSION(klon,klev) :: zh |
---|
| 216 | REAL, DIMENSION(klon,klev+1) :: zhh |
---|
| 217 | REAL, DIMENSION(klon,klev) :: epaisseur1, epaisseur2 |
---|
| 218 | |
---|
| 219 | REAL, DIMENSION(klon,klev) :: the, thu |
---|
| 220 | |
---|
[1146] | 221 | ! REAL, DIMENSION(klon,klev) :: d_deltatw, d_deltaqw |
---|
[974] | 222 | |
---|
[1403] | 223 | REAL, DIMENSION(klon,klev+1) :: omgbw |
---|
[1146] | 224 | REAL, DIMENSION(klon) :: pupper |
---|
[974] | 225 | REAL, DIMENSION(klon) :: omgtop |
---|
| 226 | REAL, DIMENSION(klon,klev) :: dp_omgbw |
---|
| 227 | REAL, DIMENSION(klon) :: ztop, dztop |
---|
| 228 | REAL, DIMENSION(klon,klev) :: alpha_up |
---|
| 229 | |
---|
| 230 | REAL, dimension(klon) :: RRe1, RRe2 |
---|
| 231 | REAL :: RRd1, RRd2 |
---|
[1403] | 232 | REAL, DIMENSION(klon,klev) :: Th1, Th2, q1, q2 |
---|
[974] | 233 | REAL, DIMENSION(klon,klev) :: D_Th1, D_Th2, D_dth |
---|
| 234 | REAL, DIMENSION(klon,klev) :: D_q1, D_q2, D_dq |
---|
| 235 | REAL, DIMENSION(klon,klev) :: omgbdq |
---|
| 236 | |
---|
| 237 | REAL, dimension(klon) :: ff, gg |
---|
| 238 | REAL, dimension(klon) :: wape2, Cstar2, heff |
---|
| 239 | |
---|
| 240 | REAL, DIMENSION(klon,klev) :: Crep |
---|
| 241 | REAL Crep_upper, Crep_sol |
---|
| 242 | |
---|
[1403] | 243 | REAL, DIMENSION(klon,klev) :: ppi |
---|
| 244 | |
---|
| 245 | ccc nrlmd |
---|
| 246 | real, dimension(klon) :: death_rate,nat_rate |
---|
| 247 | real, dimension(klon,klev) :: entr |
---|
| 248 | real, dimension(klon,klev) :: detr |
---|
| 249 | |
---|
[974] | 250 | C------------------------------------------------------------------------- |
---|
| 251 | c Initialisations |
---|
| 252 | c------------------------------------------------------------------------- |
---|
| 253 | |
---|
| 254 | c print*, 'wake initialisations' |
---|
| 255 | |
---|
| 256 | c Essais d'initialisation avec sigmaw = 0.02 et hw = 10. |
---|
| 257 | c------------------------------------------------------------------------- |
---|
| 258 | |
---|
[1403] | 259 | DATA wapecut,sigmad, hwmin /5.,.02,10./ |
---|
| 260 | ccc nrlmd |
---|
| 261 | DATA sigmaw_max /0.4/ |
---|
| 262 | DATA dens_rate /0.1/ |
---|
| 263 | ccc |
---|
[974] | 264 | C Longueur de maille (en m) |
---|
| 265 | c------------------------------------------------------------------------- |
---|
| 266 | |
---|
| 267 | c ALON = 3.e5 |
---|
| 268 | ALON = 1.e6 |
---|
| 269 | |
---|
| 270 | |
---|
| 271 | C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei) |
---|
| 272 | c |
---|
| 273 | c coefgw : Coefficient pour les ondes de gravité |
---|
| 274 | c stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) |
---|
| 275 | c wdens : Densité de poche froide par maille |
---|
| 276 | c------------------------------------------------------------------------- |
---|
| 277 | |
---|
[1403] | 278 | ccc nrlmd coefgw=10 |
---|
[974] | 279 | c coefgw=1 |
---|
[1403] | 280 | c wdens0 = 1.0/(alon**2) |
---|
| 281 | ccc nrlmd wdens = 1.0/(alon**2) |
---|
| 282 | ccc nrlmd stark = 0.50 |
---|
[974] | 283 | cCRtest |
---|
[1403] | 284 | ccc nrlmd alpk=0.1 |
---|
| 285 | c alpk = 1.0 |
---|
[974] | 286 | c alpk = 0.5 |
---|
| 287 | c alpk = 0.05 |
---|
[1403] | 288 | c |
---|
| 289 | stark = 0.33 |
---|
| 290 | Alpk = 0.25 |
---|
| 291 | wdens_ref = 8.e-12 |
---|
| 292 | coefgw = 4. |
---|
[974] | 293 | Crep_upper=0.9 |
---|
| 294 | Crep_sol=1.0 |
---|
| 295 | |
---|
[1403] | 296 | ccc nrlmd Lecture du fichier wake_param.data |
---|
| 297 | OPEN(99,file='wake_param.data',status='old', |
---|
| 298 | $ form='formatted',err=9999) |
---|
| 299 | READ(99,*,end=9998) stark |
---|
| 300 | READ(99,*,end=9998) Alpk |
---|
| 301 | READ(99,*,end=9998) wdens_ref |
---|
| 302 | READ(99,*,end=9998) coefgw |
---|
| 303 | 9998 Continue |
---|
| 304 | CLOSE(99) |
---|
| 305 | 9999 Continue |
---|
| 306 | c |
---|
| 307 | c Initialisation de toutes des densites a wdens_ref. |
---|
| 308 | c Les densites peuvent evoluer si les poches debordent |
---|
| 309 | c (voir au tout debut de la boucle sur les substeps) |
---|
| 310 | wdens = wdens_ref |
---|
| 311 | c |
---|
| 312 | c print*,'stark',stark |
---|
| 313 | c print*,'alpk',alpk |
---|
| 314 | c print*,'wdens',wdens |
---|
| 315 | c print*,'coefgw',coefgw |
---|
| 316 | ccc |
---|
[974] | 317 | C Minimum value for |T_wake - T_undist|. Used for wake top definition |
---|
| 318 | c------------------------------------------------------------------------- |
---|
| 319 | |
---|
| 320 | delta_t_min = 0.2 |
---|
| 321 | |
---|
| 322 | C 1. - Save initial values and initialize tendencies |
---|
| 323 | C -------------------------------------------------- |
---|
| 324 | |
---|
| 325 | DO k=1,klev |
---|
| 326 | DO i=1, klon |
---|
[1403] | 327 | ppi(i,k)=pi(i,k) |
---|
[974] | 328 | deltatw0(i,k) = deltatw(i,k) |
---|
| 329 | deltaqw0(i,k)= deltaqw(i,k) |
---|
| 330 | te(i,k) = te0(i,k) |
---|
| 331 | qe(i,k) = qe0(i,k) |
---|
| 332 | dtls(i,k) = 0. |
---|
| 333 | dqls(i,k) = 0. |
---|
| 334 | d_deltat_gw(i,k)=0. |
---|
[1146] | 335 | d_te(i,k) = 0. |
---|
| 336 | d_qe(i,k) = 0. |
---|
| 337 | d_deltatw(i,k) = 0. |
---|
| 338 | d_deltaqw(i,k) = 0. |
---|
[974] | 339 | !IM 060508 beg |
---|
| 340 | d_deltatw2(i,k)=0. |
---|
| 341 | d_deltaqw2(i,k)=0. |
---|
| 342 | !IM 060508 end |
---|
| 343 | ENDDO |
---|
| 344 | ENDDO |
---|
| 345 | c sigmaw1=sigmaw |
---|
| 346 | c IF (sigd_con.GT.sigmaw1) THEN |
---|
| 347 | c print*, 'sigmaw,sigd_con', sigmaw, sigd_con |
---|
| 348 | c ENDIF |
---|
| 349 | DO i=1, klon |
---|
| 350 | cc sigmaw(i) = amax1(sigmaw(i),sigd_con(i)) |
---|
| 351 | sigmaw(i) = amax1(sigmaw(i),sigmad) |
---|
| 352 | sigmaw(i) = amin1(sigmaw(i),0.99) |
---|
| 353 | sigmaw0(i) = sigmaw(i) |
---|
[1146] | 354 | wape(i) = 0. |
---|
| 355 | wape2(i) = 0. |
---|
| 356 | d_sigmaw(i) = 0. |
---|
| 357 | ktopw(i) = 0 |
---|
[974] | 358 | ENDDO |
---|
| 359 | C |
---|
| 360 | C |
---|
| 361 | C 2. - Prognostic part |
---|
| 362 | C -------------------- |
---|
| 363 | C |
---|
| 364 | C |
---|
| 365 | C 2.1 - Undisturbed area and Wake integrals |
---|
| 366 | C --------------------------------------------------------- |
---|
| 367 | |
---|
| 368 | DO i=1, klon |
---|
| 369 | z(i) = 0. |
---|
| 370 | ktop(i)=0 |
---|
| 371 | kupper(i) = 0 |
---|
| 372 | sum_thu(i) = 0. |
---|
| 373 | sum_tu(i) = 0. |
---|
| 374 | sum_qu(i) = 0. |
---|
| 375 | sum_thvu(i) = 0. |
---|
| 376 | sum_dth(i) = 0. |
---|
| 377 | sum_dq(i) = 0. |
---|
| 378 | sum_rho(i) = 0. |
---|
| 379 | sum_dtdwn(i) = 0. |
---|
| 380 | sum_dqdwn(i) = 0. |
---|
| 381 | |
---|
| 382 | av_thu(i) = 0. |
---|
| 383 | av_tu(i) =0. |
---|
| 384 | av_qu(i) =0. |
---|
| 385 | av_thvu(i) = 0. |
---|
| 386 | av_dth(i) = 0. |
---|
| 387 | av_dq(i) = 0. |
---|
| 388 | av_rho(i) =0. |
---|
| 389 | av_dtdwn(i) =0. |
---|
| 390 | av_dqdwn(i) = 0. |
---|
| 391 | ENDDO |
---|
| 392 | c |
---|
| 393 | c Distance between wakes |
---|
| 394 | DO i = 1,klon |
---|
[1403] | 395 | LL(i) = (1-sqrt(sigmaw(i)))/sqrt(wdens(i)) |
---|
[974] | 396 | ENDDO |
---|
| 397 | C Potential temperatures and humidity |
---|
| 398 | c---------------------------------------------------------- |
---|
| 399 | DO k =1,klev |
---|
| 400 | DO i=1, klon |
---|
[1403] | 401 | ! write(*,*)'wake 1',i,k,rd,te(i,k) |
---|
[974] | 402 | rho(i,k) = p(i,k)/(rd*te(i,k)) |
---|
[1403] | 403 | ! write(*,*)'wake 2',rho(i,k) |
---|
[974] | 404 | IF(k .eq. 1) THEN |
---|
[1403] | 405 | ! write(*,*)'wake 3',i,k,rd,te(i,k) |
---|
[974] | 406 | rhoh(i,k) = ph(i,k)/(rd*te(i,k)) |
---|
[1403] | 407 | ! write(*,*)'wake 4',i,k,rd,te(i,k) |
---|
[974] | 408 | zhh(i,k)=0 |
---|
| 409 | ELSE |
---|
[1403] | 410 | ! write(*,*)'wake 5',rd,(te(i,k)+te(i,k-1)) |
---|
[974] | 411 | rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1))) |
---|
[1403] | 412 | ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1) |
---|
[974] | 413 | zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1) |
---|
| 414 | ENDIF |
---|
[1403] | 415 | ! write(*,*)'wake 7',ppi(i,k) |
---|
[974] | 416 | the(i,k) = te(i,k)/ppi(i,k) |
---|
| 417 | thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k) |
---|
| 418 | tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i) |
---|
| 419 | qu(i,k) = qe(i,k) - deltaqw(i,k)*sigmaw(i) |
---|
[1403] | 420 | ! write(*,*)'wake 8',(rd*(te(i,k)+deltatw(i,k))) |
---|
| 421 | rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k))) |
---|
[974] | 422 | dth(i,k) = deltatw(i,k)/ppi(i,k) |
---|
| 423 | ENDDO |
---|
| 424 | ENDDO |
---|
| 425 | |
---|
| 426 | DO k = 1, klev-1 |
---|
| 427 | DO i=1, klon |
---|
| 428 | IF(k.eq.1) THEN |
---|
| 429 | N2(i,k)=0 |
---|
| 430 | ELSE |
---|
| 431 | N2(i,k)=amax1(0.,-RG**2/the(i,k)*rho(i,k)*(the(i,k+1)- |
---|
| 432 | $ the(i,k-1))/(p(i,k+1)-p(i,k-1))) |
---|
| 433 | ENDIF |
---|
| 434 | ZH(i,k)=(zhh(i,k)+zhh(i,k+1))/2 |
---|
| 435 | |
---|
| 436 | Cgw(i,k)=sqrt(N2(i,k))*ZH(i,k) |
---|
| 437 | Tgw(i,k)=coefgw*Cgw(i,k)/LL(i) |
---|
| 438 | ENDDO |
---|
| 439 | ENDDO |
---|
| 440 | |
---|
| 441 | DO i=1, klon |
---|
| 442 | N2(i,klev)=0 |
---|
| 443 | ZH(i,klev)=0 |
---|
| 444 | Cgw(i,klev)=0 |
---|
| 445 | Tgw(i,klev)=0 |
---|
| 446 | ENDDO |
---|
| 447 | |
---|
| 448 | c Calcul de la masse volumique moyenne de la colonne (bdlmd) |
---|
| 449 | c----------------------------------------------------------------- |
---|
| 450 | |
---|
| 451 | DO k=1,klev |
---|
| 452 | DO i=1, klon |
---|
| 453 | epaisseur1(i,k)=0. |
---|
| 454 | epaisseur2(i,k)=0. |
---|
| 455 | ENDDO |
---|
| 456 | ENDDO |
---|
| 457 | |
---|
| 458 | DO i=1, klon |
---|
| 459 | epaisseur1(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1. |
---|
| 460 | epaisseur2(i,1)= -(ph(i,2)-ph(i,1))/(rho(i,1)*rg)+1. |
---|
| 461 | rhow_moyen(i,1) = rhow(i,1) |
---|
| 462 | ENDDO |
---|
| 463 | |
---|
| 464 | DO k = 2, klev |
---|
| 465 | DO i=1, klon |
---|
| 466 | epaisseur1(i,k)= -(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg) +1. |
---|
| 467 | epaisseur2(i,k)=epaisseur2(i,k-1)+epaisseur1(i,k) |
---|
| 468 | rhow_moyen(i,k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+ |
---|
| 469 | $ rhow(i,k)*epaisseur1(i,k))/epaisseur2(i,k) |
---|
| 470 | ENDDO |
---|
| 471 | ENDDO |
---|
| 472 | |
---|
| 473 | C |
---|
| 474 | C Choose an integration bound well above wake top |
---|
| 475 | c----------------------------------------------------------------- |
---|
| 476 | c |
---|
| 477 | C Pupper = 50000. ! melting level |
---|
[1146] | 478 | c Pupper = 60000. |
---|
[974] | 479 | c Pupper = 80000. ! essais pour case_e |
---|
[1146] | 480 | DO i = 1,klon |
---|
[1403] | 481 | Pupper(i) = 0.6*ph(i,1) |
---|
| 482 | Pupper(i) = max(Pupper(i), 45000.) |
---|
| 483 | ccc Pupper(i) = 60000. |
---|
[1146] | 484 | ENDDO |
---|
| 485 | |
---|
[974] | 486 | C |
---|
| 487 | C Determine Wake top pressure (Ptop) from buoyancy integral |
---|
| 488 | C -------------------------------------------------------- |
---|
| 489 | c |
---|
| 490 | c-1/ Pressure of the level where dth becomes less than delta_t_min. |
---|
| 491 | |
---|
| 492 | DO i=1,klon |
---|
| 493 | ptop_provis(i)=ph(i,1) |
---|
| 494 | ENDDO |
---|
| 495 | DO k= 2,klev |
---|
| 496 | DO i=1,klon |
---|
| 497 | c |
---|
| 498 | cIM v3JYG; ptop_provis(i).LT. ph(i,1) |
---|
| 499 | c |
---|
| 500 | IF (dth(i,k) .GT. -delta_t_min .and. |
---|
| 501 | $ dth(i,k-1).LT. -delta_t_min .and. |
---|
| 502 | $ ptop_provis(i).EQ. ph(i,1)) THEN |
---|
| 503 | ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) |
---|
| 504 | $ - (dth(i,k-1)+delta_t_min)*p(i,k)) / |
---|
| 505 | $ (dth(i,k) - dth(i,k-1)) |
---|
| 506 | ENDIF |
---|
| 507 | ENDDO |
---|
| 508 | ENDDO |
---|
| 509 | |
---|
| 510 | c-2/ dth integral |
---|
| 511 | |
---|
| 512 | DO i=1,klon |
---|
| 513 | sum_dth(i) = 0. |
---|
| 514 | dthmin(i) = -delta_t_min |
---|
| 515 | z(i) = 0. |
---|
| 516 | ENDDO |
---|
| 517 | |
---|
| 518 | DO k = 1,klev |
---|
| 519 | DO i=1,klon |
---|
| 520 | dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg) |
---|
| 521 | IF (dz(i) .gt. 0) THEN |
---|
| 522 | z(i) = z(i)+dz(i) |
---|
| 523 | sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i) |
---|
| 524 | dthmin(i) = amin1(dthmin(i),dth(i,k)) |
---|
| 525 | ENDIF |
---|
| 526 | ENDDO |
---|
| 527 | ENDDO |
---|
| 528 | |
---|
| 529 | c-3/ height of triangle with area= sum_dth and base = dthmin |
---|
| 530 | |
---|
| 531 | DO i=1,klon |
---|
| 532 | hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5) |
---|
| 533 | hw0(i) = amax1(hwmin,hw0(i)) |
---|
| 534 | ENDDO |
---|
| 535 | |
---|
| 536 | c-4/ now, get Ptop |
---|
| 537 | |
---|
| 538 | DO i=1,klon |
---|
| 539 | z(i) = 0. |
---|
| 540 | ptop(i) = ph(i,1) |
---|
| 541 | ENDDO |
---|
| 542 | |
---|
| 543 | DO k = 1,klev |
---|
| 544 | DO i=1,klon |
---|
| 545 | dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*rg),hw0(i)-z(i)) |
---|
| 546 | IF (dz(i) .gt. 0) THEN |
---|
| 547 | z(i) = z(i)+dz(i) |
---|
| 548 | ptop(i) = ph(i,k)-rho(i,k)*rg*dz(i) |
---|
| 549 | ENDIF |
---|
| 550 | ENDDO |
---|
| 551 | ENDDO |
---|
| 552 | |
---|
| 553 | |
---|
| 554 | C-5/ Determination de ktop et kupper |
---|
| 555 | |
---|
| 556 | DO k=klev,1,-1 |
---|
| 557 | DO i=1,klon |
---|
| 558 | IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k |
---|
[1146] | 559 | IF (ph(i,k+1) .lt. pupper(i)) kupper(i)=k |
---|
[974] | 560 | ENDDO |
---|
| 561 | ENDDO |
---|
| 562 | |
---|
[1403] | 563 | c On evite kupper = 1 |
---|
| 564 | DO i=1,klon |
---|
| 565 | kupper(i) = max(kupper(i),2) |
---|
| 566 | ENDDO |
---|
| 567 | |
---|
| 568 | |
---|
[974] | 569 | c-6/ Correct ktop and ptop |
---|
| 570 | |
---|
| 571 | DO i = 1,klon |
---|
| 572 | ptop_new(i)=ptop(i) |
---|
| 573 | ENDDO |
---|
| 574 | DO k= klev,2,-1 |
---|
| 575 | DO i=1,klon |
---|
| 576 | IF (k .LE. ktop(i) .and. |
---|
| 577 | $ ptop_new(i) .EQ. ptop(i) .and. |
---|
| 578 | $ dth(i,k) .GT. -delta_t_min .and. |
---|
| 579 | $ dth(i,k-1).LT. -delta_t_min) THEN |
---|
| 580 | ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) |
---|
| 581 | $ - (dth(i,k-1)+delta_t_min)*p(i,k)) / |
---|
| 582 | $ (dth(i,k) - dth(i,k-1)) |
---|
| 583 | ENDIF |
---|
| 584 | ENDDO |
---|
| 585 | ENDDO |
---|
| 586 | |
---|
| 587 | DO i=1,klon |
---|
| 588 | ptop(i) = ptop_new(i) |
---|
| 589 | ENDDO |
---|
| 590 | |
---|
| 591 | DO k=klev,1,-1 |
---|
| 592 | DO i=1,klon |
---|
| 593 | IF (ph(i,k+1) .lt. ptop(i)) ktop(i)=k |
---|
| 594 | ENDDO |
---|
| 595 | ENDDO |
---|
| 596 | c |
---|
| 597 | c-5/ Set deltatw & deltaqw to 0 above kupper |
---|
| 598 | c |
---|
| 599 | DO k = 1,klev |
---|
| 600 | DO i=1,klon |
---|
| 601 | IF (k.GE. kupper(i)) THEN |
---|
| 602 | deltatw(i,k) = 0. |
---|
| 603 | deltaqw(i,k) = 0. |
---|
| 604 | ENDIF |
---|
| 605 | ENDDO |
---|
| 606 | ENDDO |
---|
| 607 | c |
---|
| 608 | C |
---|
| 609 | C Vertical gradient of LS omega |
---|
| 610 | C |
---|
| 611 | DO k = 1,klev |
---|
| 612 | DO i=1,klon |
---|
| 613 | IF (k.LE. kupper(i)) THEN |
---|
| 614 | dp_omgb(i,k) = (omgb(i,k+1) - omgb(i,k))/(ph(i,k+1)-ph(i,k)) |
---|
| 615 | ENDIF |
---|
| 616 | ENDDO |
---|
| 617 | ENDDO |
---|
| 618 | C |
---|
| 619 | C Integrals (and wake top level number) |
---|
| 620 | C -------------------------------------- |
---|
| 621 | C |
---|
| 622 | C Initialize sum_thvu to 1st level virt. pot. temp. |
---|
| 623 | |
---|
| 624 | DO i=1,klon |
---|
| 625 | z(i) = 1. |
---|
| 626 | dz(i) = 1. |
---|
| 627 | sum_thvu(i) = thu(i,1)*(1.+eps*qu(i,1))*dz(i) |
---|
| 628 | sum_dth(i) = 0. |
---|
| 629 | ENDDO |
---|
| 630 | |
---|
| 631 | DO k = 1,klev |
---|
| 632 | DO i=1,klon |
---|
| 633 | dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) |
---|
| 634 | IF (dz(i) .GT. 0) THEN |
---|
| 635 | z(i) = z(i)+dz(i) |
---|
| 636 | sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i) |
---|
| 637 | sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i) |
---|
| 638 | sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i) |
---|
| 639 | sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i) |
---|
| 640 | sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i) |
---|
| 641 | sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i) |
---|
| 642 | sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i) |
---|
| 643 | sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i) |
---|
| 644 | sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i) |
---|
| 645 | ENDIF |
---|
| 646 | ENDDO |
---|
| 647 | ENDDO |
---|
| 648 | c |
---|
| 649 | DO i=1,klon |
---|
| 650 | hw0(i) = z(i) |
---|
| 651 | ENDDO |
---|
| 652 | c |
---|
| 653 | C |
---|
| 654 | C 2.1 - WAPE and mean forcing computation |
---|
| 655 | C --------------------------------------- |
---|
| 656 | C |
---|
| 657 | C --------------------------------------- |
---|
| 658 | C |
---|
| 659 | C Means |
---|
| 660 | |
---|
| 661 | DO i=1,klon |
---|
| 662 | av_thu(i) = sum_thu(i)/hw0(i) |
---|
| 663 | av_tu(i) = sum_tu(i)/hw0(i) |
---|
| 664 | av_qu(i) = sum_qu(i)/hw0(i) |
---|
| 665 | av_thvu(i) = sum_thvu(i)/hw0(i) |
---|
| 666 | c av_thve = sum_thve/hw0 |
---|
| 667 | av_dth(i) = sum_dth(i)/hw0(i) |
---|
| 668 | av_dq(i) = sum_dq(i)/hw0(i) |
---|
| 669 | av_rho(i) = sum_rho(i)/hw0(i) |
---|
| 670 | av_dtdwn(i) = sum_dtdwn(i)/hw0(i) |
---|
| 671 | av_dqdwn(i) = sum_dqdwn(i)/hw0(i) |
---|
| 672 | |
---|
| 673 | wape(i) = - rg*hw0(i)*(av_dth(i) |
---|
| 674 | $ + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)* |
---|
| 675 | $ av_dq(i) ))/av_thvu(i) |
---|
| 676 | ENDDO |
---|
| 677 | C |
---|
| 678 | C 2.2 Prognostic variable update |
---|
| 679 | C ------------------------------ |
---|
| 680 | C |
---|
| 681 | C Filter out bad wakes |
---|
| 682 | |
---|
| 683 | DO k = 1,klev |
---|
| 684 | DO i=1,klon |
---|
| 685 | IF ( wape(i) .LT. 0.) THEN |
---|
| 686 | deltatw(i,k) = 0. |
---|
| 687 | deltaqw(i,k) = 0. |
---|
| 688 | dth(i,k) = 0. |
---|
| 689 | ENDIF |
---|
| 690 | ENDDO |
---|
| 691 | ENDDO |
---|
| 692 | c |
---|
| 693 | DO i=1,klon |
---|
| 694 | IF ( wape(i) .LT. 0.) THEN |
---|
| 695 | wape(i) = 0. |
---|
| 696 | Cstar(i) = 0. |
---|
| 697 | hw(i) = hwmin |
---|
| 698 | sigmaw(i) = amax1(sigmad,sigd_con(i)) |
---|
| 699 | fip(i) = 0. |
---|
| 700 | gwake(i) = .FALSE. |
---|
| 701 | ELSE |
---|
| 702 | Cstar(i) = stark*sqrt(2.*wape(i)) |
---|
| 703 | gwake(i) = .TRUE. |
---|
| 704 | ENDIF |
---|
| 705 | ENDDO |
---|
[1146] | 706 | |
---|
[974] | 707 | c |
---|
[1146] | 708 | c Check qx and qw positivity |
---|
| 709 | c -------------------------- |
---|
| 710 | DO i = 1,klon |
---|
| 711 | q0_min(i)=min( (qe(i,1)-sigmaw(i)*deltaqw(i,1)), |
---|
| 712 | $ (qe(i,1)+(1.-sigmaw(i))*deltaqw(i,1)) ) |
---|
| 713 | ENDDO |
---|
| 714 | DO k = 2,klev |
---|
| 715 | DO i = 1,klon |
---|
| 716 | q1_min(i)=min( (qe(i,k)-sigmaw(i)*deltaqw(i,k)), |
---|
| 717 | $ (qe(i,k)+(1.-sigmaw(i))*deltaqw(i,k)) ) |
---|
| 718 | IF (q1_min(i).le.q0_min(i)) THEN |
---|
| 719 | q0_min(i)=q1_min(i) |
---|
| 720 | ENDIF |
---|
| 721 | ENDDO |
---|
| 722 | ENDDO |
---|
| 723 | c |
---|
| 724 | DO i = 1,klon |
---|
| 725 | OK_qx_qw(i) = q0_min(i) .GE. 0. |
---|
| 726 | alpha(i) = 1. |
---|
| 727 | ENDDO |
---|
| 728 | c |
---|
[974] | 729 | CC ----------------------------------------------------------------- |
---|
| 730 | C Sub-time-stepping |
---|
| 731 | C ----------------- |
---|
| 732 | C |
---|
| 733 | nsub=10 |
---|
| 734 | dtimesub=dtime/nsub |
---|
| 735 | c |
---|
| 736 | c------------------------------------------------------------ |
---|
| 737 | DO isubstep = 1,nsub |
---|
| 738 | c------------------------------------------------------------ |
---|
[1146] | 739 | c |
---|
| 740 | c wk_adv is the logical flag enabling wake evolution in the time advance loop |
---|
| 741 | DO i = 1,klon |
---|
| 742 | wk_adv(i) = OK_qx_qw(i) .AND. alpha(i) .GE. 1. |
---|
| 743 | ENDDO |
---|
| 744 | c |
---|
[1403] | 745 | ccc nrlmd Ajout d'un recalcul de wdens dans le cas d'un entrainement négatif de ktop à kupper -------- |
---|
| 746 | ccc On calcule pour cela une densité wdens0 pour laquelle on aurait un entrainement nul --- |
---|
[974] | 747 | DO i=1,klon |
---|
[1403] | 748 | cc print *,' isubstep,wk_adv(i),cstar(i),wape(i) ', |
---|
| 749 | cc $ isubstep,wk_adv(i),cstar(i),wape(i) |
---|
| 750 | IF (wk_adv(i) .AND. cstar(i).GT.0.01) THEN |
---|
| 751 | omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i) |
---|
| 752 | $ + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i)) |
---|
| 753 | wdens0 = ( sigmaw(i) / (4.*3.14) ) * |
---|
| 754 | $ ( (1.-sigmaw(i)) * omg(i,kupper(i)+1) / |
---|
| 755 | $ ( (ph(i,1)-pupper(i)) * cstar(i) ) ) **(2) |
---|
| 756 | IF ( wdens(i) .LE. wdens0*1.1 ) THEN |
---|
| 757 | wdens(i) = wdens0 |
---|
| 758 | ENDIF |
---|
| 759 | cc print*,'omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i) |
---|
| 760 | cc $ ,ph(i,1)-pupper(i)', |
---|
| 761 | cc $ omg(i,kupper(i)+1),wdens0,wdens(i),cstar(i) |
---|
| 762 | cc $ ,ph(i,1)-pupper(i) |
---|
[1146] | 763 | ENDIF |
---|
[974] | 764 | ENDDO |
---|
[1403] | 765 | |
---|
| 766 | ccc nrlmd |
---|
| 767 | |
---|
[974] | 768 | DO i=1,klon |
---|
[1403] | 769 | IF (wk_adv(i)) THEN |
---|
| 770 | gfl(i) = 2.*sqrt(3.14*wdens(i)*sigmaw(i)) |
---|
| 771 | sigmaw(i)=amin1(sigmaw(i),sigmaw_max) |
---|
| 772 | ENDIF |
---|
| 773 | ENDDO |
---|
| 774 | DO i=1,klon |
---|
[1146] | 775 | IF (wk_adv(i)) THEN |
---|
[1403] | 776 | ccc nrlmd Introduction du taux de mortalité des poches et test sur sigmaw_max=0.4 |
---|
| 777 | ccc d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub |
---|
| 778 | IF (sigmaw(i).ge.sigmaw_max) THEN |
---|
| 779 | death_rate(i)=gfl(i)*Cstar(i)/sigmaw(i) |
---|
| 780 | ELSE |
---|
| 781 | death_rate(i)=0. |
---|
| 782 | END IF |
---|
| 783 | d_sigmaw(i) = gfl(i)*Cstar(i)*dtimesub |
---|
| 784 | $ - death_rate(i)*sigmaw(i)*dtimesub |
---|
| 785 | c $ - nat_rate(i)*sigmaw(i)*dtimesub |
---|
| 786 | cc print*, 'd_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i), |
---|
| 787 | cc $ death_rate(i),ktop(i),kupper(i)', |
---|
| 788 | cc $ d_sigmaw(i),sigmaw(i),gfl(i),Cstar(i),wape(i), |
---|
| 789 | cc $ death_rate(i),ktop(i),kupper(i) |
---|
| 790 | |
---|
[1146] | 791 | c sigmaw(i) =sigmaw(i) + gfl(i)*Cstar(i)*dtimesub |
---|
| 792 | c sigmaw(i) =min(sigmaw(i),0.99) !!!!!!!! |
---|
[974] | 793 | c wdens = wdens0/(10.*sigmaw) |
---|
| 794 | c sigmaw =max(sigmaw,sigd_con) |
---|
| 795 | c sigmaw =max(sigmaw,sigmad) |
---|
[1146] | 796 | ENDIF |
---|
[974] | 797 | ENDDO |
---|
| 798 | C |
---|
| 799 | C |
---|
| 800 | c calcul de la difference de vitesse verticale poche - zone non perturbee |
---|
| 801 | cIM 060208 differences par rapport au code initial; init. a 0 dp_deltomg |
---|
| 802 | cIM 060208 et omg sur les niveaux de 1 a klev+1, alors que avant l'on definit |
---|
| 803 | cIM 060208 au niveau k=1..? |
---|
[1146] | 804 | DO k= 1,klev |
---|
| 805 | DO i = 1,klon |
---|
[1403] | 806 | if (wk_adv(i)) THEN !!! nrlmd |
---|
[1146] | 807 | dp_deltomg(i,k)=0. |
---|
[1403] | 808 | end if |
---|
[1146] | 809 | ENDDO |
---|
| 810 | ENDDO |
---|
[974] | 811 | DO k= 1,klev+1 |
---|
| 812 | DO i = 1,klon |
---|
[1403] | 813 | if (wk_adv(i)) THEN !!! nrlmd |
---|
[974] | 814 | omg(i,k)=0. |
---|
[1403] | 815 | end if |
---|
[974] | 816 | ENDDO |
---|
| 817 | ENDDO |
---|
| 818 | c |
---|
| 819 | DO i=1,klon |
---|
[1146] | 820 | IF (wk_adv(i)) THEN |
---|
[974] | 821 | z(i)= 0. |
---|
| 822 | omg(i,1) = 0. |
---|
| 823 | dp_deltomg(i,1) = -(gfl(i)*Cstar(i))/(sigmaw(i) * (1-sigmaw(i))) |
---|
[1146] | 824 | ENDIF |
---|
[974] | 825 | ENDDO |
---|
| 826 | c |
---|
| 827 | DO k= 2,klev |
---|
| 828 | DO i = 1,klon |
---|
[1146] | 829 | IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN |
---|
[974] | 830 | dz(i) = -(ph(i,k)-ph(i,k-1))/(rho(i,k-1)*rg) |
---|
| 831 | z(i) = z(i)+dz(i) |
---|
| 832 | dp_deltomg(i,k)= dp_deltomg(i,1) |
---|
| 833 | omg(i,k)= dp_deltomg(i,1)*z(i) |
---|
| 834 | ENDIF |
---|
| 835 | ENDDO |
---|
| 836 | ENDDO |
---|
| 837 | c |
---|
| 838 | DO i = 1,klon |
---|
[1146] | 839 | IF (wk_adv(i)) THEN |
---|
[974] | 840 | dztop(i)=-(ptop(i)-ph(i,ktop(i)))/(rho(i,ktop(i))*rg) |
---|
| 841 | ztop(i) = z(i)+dztop(i) |
---|
| 842 | omgtop(i)=dp_deltomg(i,1)*ztop(i) |
---|
[1146] | 843 | ENDIF |
---|
[974] | 844 | ENDDO |
---|
| 845 | c |
---|
| 846 | c ----------------- |
---|
| 847 | c From m/s to Pa/s |
---|
| 848 | c ----------------- |
---|
| 849 | c |
---|
| 850 | DO i=1,klon |
---|
[1146] | 851 | IF (wk_adv(i)) THEN |
---|
[974] | 852 | omgtop(i) = -rho(i,ktop(i))*rg*omgtop(i) |
---|
| 853 | dp_deltomg(i,1) = omgtop(i)/(ptop(i)-ph(i,1)) |
---|
[1146] | 854 | ENDIF |
---|
[974] | 855 | ENDDO |
---|
| 856 | c |
---|
| 857 | DO k= 1,klev |
---|
| 858 | DO i = 1,klon |
---|
[1146] | 859 | IF( wk_adv(i) .AND. k .LE. ktop(i)) THEN |
---|
[974] | 860 | omg(i,k) = - rho(i,k)*rg*omg(i,k) |
---|
| 861 | dp_deltomg(i,k) = dp_deltomg(i,1) |
---|
| 862 | ENDIF |
---|
| 863 | ENDDO |
---|
| 864 | ENDDO |
---|
| 865 | c |
---|
| 866 | c raccordement lineaire de omg de ptop a pupper |
---|
| 867 | |
---|
| 868 | DO i=1,klon |
---|
[1146] | 869 | IF ( wk_adv(i) .AND. kupper(i) .GT. ktop(i)) THEN |
---|
[974] | 870 | omg(i,kupper(i)+1) = - Rg*amdwn(i,kupper(i)+1)/sigmaw(i) |
---|
| 871 | $ + Rg*amup(i,kupper(i)+1)/(1.-sigmaw(i)) |
---|
| 872 | dp_deltomg(i,kupper(i)) = (omgtop(i)-omg(i,kupper(i)+1))/ |
---|
[1146] | 873 | $ (ptop(i)-pupper(i)) |
---|
[974] | 874 | ENDIF |
---|
| 875 | ENDDO |
---|
| 876 | c |
---|
[1403] | 877 | cc DO i=1,klon |
---|
| 878 | cc print*,'Pente entre 0 et kupper (référence)' |
---|
| 879 | cc $ ,omg(i,kupper(i)+1)/(pupper(i)-ph(i,1)) |
---|
| 880 | cc print*,'Pente entre ktop et kupper' |
---|
| 881 | cc $ ,(omg(i,kupper(i)+1)-omgtop(i))/(pupper(i)-ptop(i)) |
---|
| 882 | cc ENDDO |
---|
| 883 | cc |
---|
[974] | 884 | DO k= 1,klev |
---|
| 885 | DO i = 1,klon |
---|
[1146] | 886 | IF( wk_adv(i) .AND. k .GT. ktop(i) .AND. k .LE. kupper(i)) THEN |
---|
[974] | 887 | dp_deltomg(i,k) = dp_deltomg(i,kupper(i)) |
---|
| 888 | omg(i,k) = omgtop(i)+(ph(i,k)-ptop(i))*dp_deltomg(i,kupper(i)) |
---|
| 889 | ENDIF |
---|
| 890 | ENDDO |
---|
| 891 | ENDDO |
---|
[1403] | 892 | ccc nrlmd |
---|
| 893 | cc DO i=1,klon |
---|
| 894 | cc print*,'deltaw_ktop,deltaw_conv',omgtop(i),omg(i,kupper(i)+1) |
---|
| 895 | cc END DO |
---|
| 896 | ccc |
---|
[974] | 897 | c |
---|
[1146] | 898 | c |
---|
[974] | 899 | c-- Compute wake average vertical velocity omgbw |
---|
| 900 | c |
---|
| 901 | c |
---|
| 902 | DO k = 1,klev+1 |
---|
| 903 | DO i=1,klon |
---|
[1146] | 904 | IF ( wk_adv(i)) THEN |
---|
[974] | 905 | omgbw(i,k) = omgb(i,k)+(1.-sigmaw(i))*omg(i,k) |
---|
[1146] | 906 | ENDIF |
---|
[974] | 907 | ENDDO |
---|
| 908 | ENDDO |
---|
| 909 | c-- and its vertical gradient dp_omgbw |
---|
| 910 | c |
---|
| 911 | DO k = 1,klev |
---|
| 912 | DO i=1,klon |
---|
[1146] | 913 | IF ( wk_adv(i)) THEN |
---|
[974] | 914 | dp_omgbw(i,k) = (omgbw(i,k+1)-omgbw(i,k))/(ph(i,k+1)-ph(i,k)) |
---|
[1146] | 915 | ENDIF |
---|
[974] | 916 | ENDDO |
---|
| 917 | ENDDO |
---|
| 918 | C |
---|
| 919 | c-- Upstream coefficients for omgb velocity |
---|
| 920 | c-- (alpha_up(k) is the coefficient of the value at level k) |
---|
| 921 | c-- (1-alpha_up(k) is the coefficient of the value at level k-1) |
---|
| 922 | DO k = 1,klev |
---|
| 923 | DO i=1,klon |
---|
[1146] | 924 | IF ( wk_adv(i)) THEN |
---|
| 925 | alpha_up(i,k) = 0. |
---|
| 926 | IF (omgb(i,k) .GT. 0.) alpha_up(i,k) = 1. |
---|
| 927 | ENDIF |
---|
[974] | 928 | ENDDO |
---|
| 929 | ENDDO |
---|
| 930 | |
---|
| 931 | c Matrix expressing [The,deltatw] from [Th1,Th2] |
---|
| 932 | |
---|
| 933 | DO i=1,klon |
---|
[1146] | 934 | IF ( wk_adv(i)) THEN |
---|
| 935 | RRe1(i) = 1.-sigmaw(i) |
---|
| 936 | RRe2(i) = sigmaw(i) |
---|
| 937 | ENDIF |
---|
[974] | 938 | ENDDO |
---|
| 939 | RRd1 = -1. |
---|
| 940 | RRd2 = 1. |
---|
| 941 | c |
---|
| 942 | c-- Get [Th1,Th2], dth and [q1,q2] |
---|
| 943 | c |
---|
| 944 | DO k= 1,klev |
---|
| 945 | DO i = 1,klon |
---|
[1146] | 946 | IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN |
---|
[974] | 947 | dth(i,k) = deltatw(i,k)/ppi(i,k) |
---|
| 948 | Th1(i,k) = the(i,k) - sigmaw(i) *dth(i,k) ! undisturbed area |
---|
| 949 | Th2(i,k) = the(i,k) + (1.-sigmaw(i))*dth(i,k) ! wake |
---|
| 950 | q1(i,k) = qe(i,k) - sigmaw(i) *deltaqw(i,k) ! undisturbed area |
---|
| 951 | q2(i,k) = qe(i,k) + (1.-sigmaw(i))*deltaqw(i,k) ! wake |
---|
| 952 | ENDIF |
---|
| 953 | ENDDO |
---|
| 954 | ENDDO |
---|
| 955 | |
---|
| 956 | DO i=1,klon |
---|
[1403] | 957 | if (wk_adv(i)) then !!! nrlmd |
---|
| 958 | D_Th1(i,1) = 0. |
---|
[974] | 959 | D_Th2(i,1) = 0. |
---|
| 960 | D_dth(i,1) = 0. |
---|
| 961 | D_q1(i,1) = 0. |
---|
| 962 | D_q2(i,1) = 0. |
---|
| 963 | D_dq(i,1) = 0. |
---|
[1403] | 964 | end if |
---|
[974] | 965 | ENDDO |
---|
| 966 | |
---|
| 967 | DO k= 2,klev |
---|
| 968 | DO i = 1,klon |
---|
[1146] | 969 | IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN |
---|
[974] | 970 | D_Th1(i,k) = Th1(i,k-1)-Th1(i,k) |
---|
| 971 | D_Th2(i,k) = Th2(i,k-1)-Th2(i,k) |
---|
| 972 | D_dth(i,k) = dth(i,k-1)-dth(i,k) |
---|
| 973 | D_q1(i,k) = q1(i,k-1)-q1(i,k) |
---|
| 974 | D_q2(i,k) = q2(i,k-1)-q2(i,k) |
---|
| 975 | D_dq(i,k) = deltaqw(i,k-1)-deltaqw(i,k) |
---|
| 976 | ENDIF |
---|
| 977 | ENDDO |
---|
| 978 | ENDDO |
---|
| 979 | |
---|
| 980 | DO i=1,klon |
---|
[1146] | 981 | IF( wk_adv(i)) THEN |
---|
| 982 | omgbdth(i,1) = 0. |
---|
| 983 | omgbdq(i,1) = 0. |
---|
| 984 | ENDIF |
---|
[974] | 985 | ENDDO |
---|
| 986 | |
---|
| 987 | DO k= 2,klev |
---|
| 988 | DO i = 1,klon |
---|
[1146] | 989 | IF( wk_adv(i) .AND. k .LE. kupper(i)+1) THEN ! loop on interfaces |
---|
[974] | 990 | omgbdth(i,k) = omgb(i,k)*( dth(i,k-1) - dth(i,k)) |
---|
| 991 | omgbdq(i,k) = omgb(i,k)*(deltaqw(i,k-1) - deltaqw(i,k)) |
---|
| 992 | ENDIF |
---|
| 993 | ENDDO |
---|
| 994 | ENDDO |
---|
| 995 | c |
---|
| 996 | c----------------------------------------------------------------- |
---|
| 997 | DO k= 1,klev |
---|
| 998 | DO i = 1,klon |
---|
[1146] | 999 | IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN |
---|
[974] | 1000 | c----------------------------------------------------------------- |
---|
| 1001 | c |
---|
| 1002 | c Compute redistribution (advective) term |
---|
| 1003 | c |
---|
| 1004 | d_deltatw(i,k) = |
---|
| 1005 | $ dtimesub/(Ph(i,k)-Ph(i,k+1))*( |
---|
| 1006 | $ RRd1*omg(i,k )*sigmaw(i) *D_Th1(i,k) |
---|
| 1007 | $ -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1) |
---|
| 1008 | $ -(1.-alpha_up(i,k))*omgbdth(i,k) - alpha_up(i,k+1)* |
---|
| 1009 | $ omgbdth(i,k+1))*ppi(i,k) |
---|
| 1010 | c print*,'d_deltatw=',d_deltatw(i,k) |
---|
| 1011 | c |
---|
| 1012 | d_deltaqw(i,k) = |
---|
| 1013 | $ dtimesub/(Ph(i,k)-Ph(i,k+1))*( |
---|
| 1014 | $ RRd1*omg(i,k )*sigmaw(i) *D_q1(i,k) |
---|
| 1015 | $ -RRd2*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1) |
---|
| 1016 | $ -(1.-alpha_up(i,k))*omgbdq(i,k) - alpha_up(i,k+1)* |
---|
| 1017 | $ omgbdq(i,k+1)) |
---|
| 1018 | c print*,'d_deltaqw=',d_deltaqw(i,k) |
---|
| 1019 | c |
---|
| 1020 | c and increment large scale tendencies |
---|
| 1021 | c |
---|
[1146] | 1022 | |
---|
| 1023 | c |
---|
| 1024 | C |
---|
| 1025 | CC ----------------------------------------------------------------- |
---|
| 1026 | d_te(i,k) = dtimesub*( |
---|
[974] | 1027 | $ ( RRe1(i)*omg(i,k )*sigmaw(i) *D_Th1(i,k) |
---|
| 1028 | $ -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_Th2(i,k+1) ) |
---|
| 1029 | $ /(Ph(i,k)-Ph(i,k+1)) |
---|
[1403] | 1030 | ccc nrlmd $ -sigmaw(i)*(1.-sigmaw(i))*dth(i,k)*dp_deltomg(i,k) |
---|
| 1031 | $ -sigmaw(i)*(1.-sigmaw(i))*dth(i,k) |
---|
| 1032 | $ *(omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1)) |
---|
| 1033 | ccc |
---|
[974] | 1034 | $ )*ppi(i,k) |
---|
| 1035 | c |
---|
[1146] | 1036 | d_qe(i,k) = dtimesub*( |
---|
[974] | 1037 | $ ( RRe1(i)*omg(i,k )*sigmaw(i) *D_q1(i,k) |
---|
| 1038 | $ -RRe2(i)*omg(i,k+1)*(1.-sigmaw(i))*D_q2(i,k+1) ) |
---|
| 1039 | $ /(Ph(i,k)-Ph(i,k+1)) |
---|
[1403] | 1040 | ccc nrlmd $ -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k)*dp_deltomg(i,k) |
---|
| 1041 | $ -sigmaw(i)*(1.-sigmaw(i))*deltaqw(i,k) |
---|
| 1042 | $ *(omg(i,k)-omg(i,k+1))/(Ph(i,k)-Ph(i,k+1)) |
---|
| 1043 | ccc |
---|
[974] | 1044 | $ ) |
---|
[1403] | 1045 | ccc nrlmd |
---|
| 1046 | ELSE IF(wk_adv(i) .AND. k .EQ. kupper(i)) THEN |
---|
[1277] | 1047 | d_te(i,k) = dtimesub*( |
---|
[1403] | 1048 | $ ( RRe1(i)*omg(i,k )*sigmaw(i) *D_Th1(i,k) |
---|
| 1049 | $ /(Ph(i,k)-Ph(i,k+1))) |
---|
| 1050 | $ )*ppi(i,k) |
---|
[1277] | 1051 | |
---|
| 1052 | d_qe(i,k) = dtimesub*( |
---|
[1403] | 1053 | $ ( RRe1(i)*omg(i,k )*sigmaw(i) *D_q1(i,k) |
---|
| 1054 | $ /(Ph(i,k)-Ph(i,k+1))) |
---|
| 1055 | $ ) |
---|
| 1056 | |
---|
[974] | 1057 | ENDIF |
---|
[1403] | 1058 | ccc |
---|
[974] | 1059 | ENDDO |
---|
| 1060 | ENDDO |
---|
| 1061 | c------------------------------------------------------------------ |
---|
| 1062 | C |
---|
| 1063 | C Increment state variables |
---|
| 1064 | |
---|
| 1065 | DO k= 1,klev |
---|
| 1066 | DO i = 1,klon |
---|
[1403] | 1067 | ccc nrlmd IF( wk_adv(i) .AND. k .LE. kupper(i)-1) THEN |
---|
| 1068 | IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN |
---|
| 1069 | ccc |
---|
| 1070 | |
---|
| 1071 | |
---|
[974] | 1072 | c |
---|
| 1073 | c Coefficient de répartition |
---|
| 1074 | |
---|
| 1075 | Crep(i,k)=Crep_sol*(ph(i,kupper(i))-ph(i,k))/(ph(i,kupper(i)) |
---|
| 1076 | $ -ph(i,1)) |
---|
| 1077 | Crep(i,k)=Crep(i,k)+Crep_upper*(ph(i,1)-ph(i,k))/(p(i,1)- |
---|
| 1078 | $ ph(i,kupper(i))) |
---|
| 1079 | |
---|
| 1080 | |
---|
| 1081 | c Reintroduce compensating subsidence term. |
---|
| 1082 | |
---|
| 1083 | c dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw |
---|
| 1084 | c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)) |
---|
| 1085 | c . /(1-sigmaw) |
---|
| 1086 | c dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw |
---|
| 1087 | c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)) |
---|
| 1088 | c . /(1-sigmaw) |
---|
| 1089 | c |
---|
| 1090 | c dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw |
---|
| 1091 | c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k)) |
---|
| 1092 | c . /(1-sigmaw) |
---|
| 1093 | c dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw |
---|
| 1094 | c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k)) |
---|
| 1095 | c . /(1-sigmaw) |
---|
| 1096 | |
---|
| 1097 | dtKE(i,k)=(dtdwn(i,k)/sigmaw(i) - dta(i,k)/(1.-sigmaw(i))) |
---|
| 1098 | dqKE(i,k)=(dqdwn(i,k)/sigmaw(i) - dqa(i,k)/(1.-sigmaw(i))) |
---|
[1403] | 1099 | c print*,'dtKE= ',dtKE(i,k),' dqKE= ',dqKE(i,k) |
---|
[974] | 1100 | c |
---|
| 1101 | dtPBL(i,k)=(wdtPBL(i,k)/sigmaw(i) - udtPBL(i,k)/(1.-sigmaw(i))) |
---|
| 1102 | dqPBL(i,k)=(wdqPBL(i,k)/sigmaw(i) - udqPBL(i,k)/(1.-sigmaw(i))) |
---|
[1403] | 1103 | c print*,'dtPBL= ',dtPBL(i,k),' dqPBL= ',dqPBL(i,k) |
---|
[974] | 1104 | c |
---|
[1403] | 1105 | ccc nrlmd Prise en compte du taux de mortalité |
---|
| 1106 | ccc Définitions de entr, detr |
---|
| 1107 | detr(i,k)=0. |
---|
[974] | 1108 | |
---|
[1403] | 1109 | entr(i,k)=detr(i,k)+gfl(i)*cstar(i)+ |
---|
| 1110 | $ sigmaw(i)*(1.-sigmaw(i))*dp_deltomg(i,k) |
---|
[974] | 1111 | |
---|
[1403] | 1112 | spread(i,k) = (entr(i,k)-detr(i,k))/sigmaw(i) |
---|
| 1113 | ccc spread(i,k) = (1.-sigmaw(i))*dp_deltomg(i,k)+gfl(i)*Cstar(i)/ |
---|
| 1114 | ccc $ sigmaw(i) |
---|
| 1115 | |
---|
| 1116 | |
---|
[974] | 1117 | c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei |
---|
| 1118 | |
---|
[1403] | 1119 | ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltat_gw(i,k), |
---|
| 1120 | ! & Tgw(i,k),deltatw(i,k) |
---|
| 1121 | d_deltat_gw(i,k)=d_deltat_gw(i,k)-Tgw(i,k)*deltatw(i,k)* |
---|
[974] | 1122 | $ dtimesub |
---|
[1403] | 1123 | ! write(lunout,*)'wake.F ',i,k, dtimesub,d_deltatw(i,k) |
---|
[974] | 1124 | ff(i)=d_deltatw(i,k)/dtimesub |
---|
| 1125 | |
---|
| 1126 | c Sans GW |
---|
| 1127 | c |
---|
[1403] | 1128 | c deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k)) |
---|
[974] | 1129 | c |
---|
| 1130 | c GW formule 1 |
---|
| 1131 | c |
---|
| 1132 | c deltatw(k) = deltatw(k)+dtimesub* |
---|
| 1133 | c $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k)) |
---|
| 1134 | c |
---|
| 1135 | c GW formule 2 |
---|
| 1136 | |
---|
| 1137 | IF (dtimesub*Tgw(i,k).lt.1.e-10) THEN |
---|
[1146] | 1138 | d_deltatw(i,k) = dtimesub* |
---|
[1403] | 1139 | $ (ff(i)+dtKE(i,k)+dtPBL(i,k) |
---|
| 1140 | ccc $ -spread(i,k)*deltatw(i,k) |
---|
| 1141 | $ - entr(i,k)*deltatw(i,k)/sigmaw(i) |
---|
| 1142 | $ - (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k) |
---|
| 1143 | $ / (1.-sigmaw(i)) |
---|
| 1144 | ccc |
---|
| 1145 | $ -Tgw(i,k)*deltatw(i,k)) |
---|
[974] | 1146 | ELSE |
---|
[1146] | 1147 | d_deltatw(i,k) = 1/Tgw(i,k)*(1-exp(-dtimesub* |
---|
[1403] | 1148 | $ Tgw(i,k)))* |
---|
| 1149 | $ (ff(i)+dtKE(i,k)+dtPBL(i,k) |
---|
| 1150 | ccc $ -spread(i,k)*deltatw(i,k) |
---|
| 1151 | $ - entr(i,k)*deltatw(i,k)/sigmaw(i) |
---|
| 1152 | $ - (death_rate(i)*sigmaw(i)+detr(i,k))*deltatw(i,k) |
---|
| 1153 | $ / (1.-sigmaw(i)) |
---|
| 1154 | ccc |
---|
| 1155 | $ -Tgw(i,k)*deltatw(i,k)) |
---|
[974] | 1156 | ENDIF |
---|
[1146] | 1157 | |
---|
[974] | 1158 | dth(i,k) = deltatw(i,k)/ppi(i,k) |
---|
| 1159 | |
---|
| 1160 | gg(i)=d_deltaqw(i,k)/dtimesub |
---|
| 1161 | |
---|
[1146] | 1162 | d_deltaqw(i,k) = dtimesub*(gg(i)+ dqKE(i,k)+dqPBL(i,k) |
---|
[1403] | 1163 | ccc $ -spread(i,k)*deltaqw(i,k)) |
---|
| 1164 | $ - entr(i,k)*deltaqw(i,k)/sigmaw(i) |
---|
| 1165 | $ - (death_rate(i)*sigmaw(i)+detr(i,k))*deltaqw(i,k) |
---|
| 1166 | $ /(1.-sigmaw(i))) |
---|
| 1167 | ccc |
---|
[974] | 1168 | |
---|
[1403] | 1169 | ccc nrlmd |
---|
| 1170 | ccc d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k) |
---|
| 1171 | ccc d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k) |
---|
| 1172 | ccc |
---|
[974] | 1173 | ENDIF |
---|
| 1174 | ENDDO |
---|
| 1175 | ENDDO |
---|
| 1176 | |
---|
[1146] | 1177 | C |
---|
| 1178 | C Scale tendencies so that water vapour remains positive in w and x. |
---|
| 1179 | C |
---|
[1403] | 1180 | call wake_vec_modulation(klon,klev,wk_adv,epsilon,qe,d_qe,deltaqw, |
---|
[1146] | 1181 | $ d_deltaqw,sigmaw,d_sigmaw,alpha) |
---|
| 1182 | c |
---|
[1403] | 1183 | ccc nrlmd |
---|
| 1184 | cc print*,'alpha' |
---|
| 1185 | cc do i=1,klon |
---|
| 1186 | cc print*,alpha(i) |
---|
| 1187 | cc end do |
---|
| 1188 | ccc |
---|
[1146] | 1189 | DO k = 1,klev |
---|
| 1190 | DO i = 1,klon |
---|
| 1191 | IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN |
---|
| 1192 | d_te(i,k)=alpha(i)*d_te(i,k) |
---|
| 1193 | d_qe(i,k)=alpha(i)*d_qe(i,k) |
---|
| 1194 | d_deltatw(i,k)=alpha(i)*d_deltatw(i,k) |
---|
| 1195 | d_deltaqw(i,k)=alpha(i)*d_deltaqw(i,k) |
---|
| 1196 | d_deltat_gw(i,k)=alpha(i)*d_deltat_gw(i,k) |
---|
| 1197 | ENDIF |
---|
| 1198 | ENDDO |
---|
| 1199 | ENDDO |
---|
| 1200 | DO i = 1,klon |
---|
| 1201 | IF( wk_adv(i)) THEN |
---|
| 1202 | d_sigmaw(i)=alpha(i)*d_sigmaw(i) |
---|
| 1203 | ENDIF |
---|
| 1204 | ENDDO |
---|
| 1205 | |
---|
| 1206 | C Update large scale variables and wake variables |
---|
[974] | 1207 | cIM 060208 manque DO i + remplace DO k=1,kupper(i) |
---|
| 1208 | cIM 060208 DO k = 1,kupper(i) |
---|
| 1209 | DO k= 1,klev |
---|
| 1210 | DO i = 1,klon |
---|
[1146] | 1211 | IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN |
---|
| 1212 | dtls(i,k)=dtls(i,k)+d_te(i,k) |
---|
| 1213 | dqls(i,k)=dqls(i,k)+d_qe(i,k) |
---|
[1403] | 1214 | ccc nrlmd |
---|
| 1215 | d_deltatw2(i,k)=d_deltatw2(i,k)+d_deltatw(i,k) |
---|
| 1216 | d_deltaqw2(i,k)=d_deltaqw2(i,k)+d_deltaqw(i,k) |
---|
| 1217 | ccc |
---|
[1146] | 1218 | ENDIF |
---|
| 1219 | ENDDO |
---|
| 1220 | ENDDO |
---|
| 1221 | DO k= 1,klev |
---|
| 1222 | DO i = 1,klon |
---|
| 1223 | IF( wk_adv(i) .AND. k .LE. kupper(i)) THEN |
---|
[974] | 1224 | te(i,k) = te0(i,k) + dtls(i,k) |
---|
| 1225 | qe(i,k) = qe0(i,k) + dqls(i,k) |
---|
| 1226 | the(i,k) = te(i,k)/ppi(i,k) |
---|
[1146] | 1227 | deltatw(i,k) = deltatw(i,k)+d_deltatw(i,k) |
---|
| 1228 | deltaqw(i,k) = deltaqw(i,k)+d_deltaqw(i,k) |
---|
| 1229 | dth(i,k) = deltatw(i,k)/ppi(i,k) |
---|
[1403] | 1230 | cc print*,'k,qx,qw',k,qe(i,k)-sigmaw(i)*deltaqw(i,k) |
---|
| 1231 | cc $ ,qe(i,k)+(1-sigmaw(i))*deltaqw(i,k) |
---|
[974] | 1232 | ENDIF |
---|
| 1233 | ENDDO |
---|
| 1234 | ENDDO |
---|
[1146] | 1235 | DO i = 1,klon |
---|
| 1236 | IF( wk_adv(i)) THEN |
---|
| 1237 | sigmaw(i) = sigmaw(i)+d_sigmaw(i) |
---|
| 1238 | ENDIF |
---|
| 1239 | ENDDO |
---|
[974] | 1240 | c |
---|
| 1241 | C |
---|
| 1242 | c Determine Ptop from buoyancy integral |
---|
| 1243 | c --------------------------------------- |
---|
| 1244 | c |
---|
| 1245 | c- 1/ Pressure of the level where dth changes sign. |
---|
| 1246 | c |
---|
| 1247 | DO i=1,klon |
---|
[1146] | 1248 | IF ( wk_adv(i)) THEN |
---|
| 1249 | Ptop_provis(i)=ph(i,1) |
---|
| 1250 | ENDIF |
---|
[974] | 1251 | ENDDO |
---|
| 1252 | c |
---|
| 1253 | DO k= 2,klev |
---|
| 1254 | DO i=1,klon |
---|
[1146] | 1255 | IF ( wk_adv(i) .AND. |
---|
| 1256 | $ Ptop_provis(i) .EQ. ph(i,1) .AND. |
---|
[974] | 1257 | $ dth(i,k) .GT. -delta_t_min .and. |
---|
| 1258 | $ dth(i,k-1).LT. -delta_t_min) THEN |
---|
| 1259 | Ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) |
---|
| 1260 | $ - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k) |
---|
| 1261 | $ - dth(i,k-1)) |
---|
| 1262 | ENDIF |
---|
| 1263 | ENDDO |
---|
| 1264 | ENDDO |
---|
| 1265 | c |
---|
| 1266 | c- 2/ dth integral |
---|
| 1267 | c |
---|
| 1268 | DO i=1,klon |
---|
[1403] | 1269 | if (wk_adv(i)) then !!! nrlmd |
---|
[974] | 1270 | sum_dth(i) = 0. |
---|
| 1271 | dthmin(i) = -delta_t_min |
---|
| 1272 | z(i) = 0. |
---|
[1403] | 1273 | end if |
---|
[974] | 1274 | ENDDO |
---|
| 1275 | |
---|
| 1276 | DO k = 1,klev |
---|
| 1277 | DO i=1,klon |
---|
[1146] | 1278 | IF ( wk_adv(i)) THEN |
---|
[974] | 1279 | dz(i) = -(amax1(ph(i,k+1),Ptop_provis(i))-Ph(i,k))/(rho(i,k)*rg) |
---|
| 1280 | IF (dz(i) .gt. 0) THEN |
---|
| 1281 | z(i) = z(i)+dz(i) |
---|
| 1282 | sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i) |
---|
| 1283 | dthmin(i) = amin1(dthmin(i),dth(i,k)) |
---|
| 1284 | ENDIF |
---|
[1146] | 1285 | ENDIF |
---|
[974] | 1286 | ENDDO |
---|
| 1287 | ENDDO |
---|
| 1288 | c |
---|
| 1289 | c- 3/ height of triangle with area= sum_dth and base = dthmin |
---|
| 1290 | |
---|
| 1291 | DO i=1,klon |
---|
[1146] | 1292 | IF ( wk_adv(i)) THEN |
---|
| 1293 | hw(i) = 2.*sum_dth(i)/amin1(dthmin(i),-0.5) |
---|
| 1294 | hw(i) = amax1(hwmin,hw(i)) |
---|
| 1295 | ENDIF |
---|
[974] | 1296 | ENDDO |
---|
| 1297 | c |
---|
| 1298 | c- 4/ now, get Ptop |
---|
| 1299 | c |
---|
| 1300 | DO i=1,klon |
---|
[1403] | 1301 | if (wk_adv(i)) then !!! nrlmd |
---|
[974] | 1302 | ktop(i) = 0 |
---|
| 1303 | z(i)=0. |
---|
[1403] | 1304 | end if |
---|
[974] | 1305 | ENDDO |
---|
| 1306 | c |
---|
| 1307 | DO k = 1,klev |
---|
| 1308 | DO i=1,klon |
---|
[1146] | 1309 | IF ( wk_adv(i)) THEN |
---|
[974] | 1310 | dz(i) = amin1(-(ph(i,k+1)-Ph(i,k))/(rho(i,k)*rg),hw(i)-z(i)) |
---|
| 1311 | IF (dz(i) .gt. 0) THEN |
---|
| 1312 | z(i) = z(i)+dz(i) |
---|
| 1313 | Ptop(i) = Ph(i,k)-rho(i,k)*rg*dz(i) |
---|
| 1314 | ktop(i) = k |
---|
| 1315 | ENDIF |
---|
[1146] | 1316 | ENDIF |
---|
[974] | 1317 | ENDDO |
---|
| 1318 | ENDDO |
---|
| 1319 | c |
---|
| 1320 | c 4.5/Correct ktop and ptop |
---|
| 1321 | c |
---|
| 1322 | DO i=1,klon |
---|
[1146] | 1323 | IF ( wk_adv(i)) THEN |
---|
[974] | 1324 | Ptop_new(i)=ptop(i) |
---|
[1146] | 1325 | ENDIF |
---|
[974] | 1326 | ENDDO |
---|
| 1327 | c |
---|
| 1328 | DO k= klev,2,-1 |
---|
| 1329 | DO i=1,klon |
---|
| 1330 | cIM v3JYG; IF (k .GE. ktop(i) |
---|
[1146] | 1331 | IF ( wk_adv(i) .AND. |
---|
| 1332 | $ k .LE. ktop(i) .AND. |
---|
[974] | 1333 | $ ptop_new(i) .EQ. ptop(i) .AND. |
---|
| 1334 | $ dth(i,k) .GT. -delta_t_min .and. |
---|
| 1335 | $ dth(i,k-1).LT. -delta_t_min) THEN |
---|
| 1336 | Ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) |
---|
| 1337 | $ - (dth(i,k-1)+delta_t_min)*p(i,k)) /(dth(i,k) |
---|
| 1338 | $ - dth(i,k-1)) |
---|
| 1339 | ENDIF |
---|
| 1340 | ENDDO |
---|
| 1341 | ENDDO |
---|
| 1342 | c |
---|
| 1343 | c |
---|
| 1344 | DO i=1,klon |
---|
[1146] | 1345 | IF ( wk_adv(i)) THEN |
---|
| 1346 | ptop(i) = ptop_new(i) |
---|
| 1347 | ENDIF |
---|
[974] | 1348 | ENDDO |
---|
| 1349 | |
---|
| 1350 | DO k=klev,1,-1 |
---|
| 1351 | DO i=1,klon |
---|
[1403] | 1352 | if (wk_adv(i)) then !!! nrlmd |
---|
[974] | 1353 | IF (ph(i,k+1) .LT. ptop(i)) ktop(i)=k |
---|
[1403] | 1354 | end if |
---|
[974] | 1355 | ENDDO |
---|
| 1356 | ENDDO |
---|
| 1357 | c |
---|
| 1358 | c 5/ Set deltatw & deltaqw to 0 above kupper |
---|
| 1359 | c |
---|
| 1360 | DO k = 1,klev |
---|
| 1361 | DO i=1,klon |
---|
[1146] | 1362 | IF ( wk_adv(i) .AND. k .GE. kupper(i)) THEN |
---|
[974] | 1363 | deltatw(i,k) = 0. |
---|
| 1364 | deltaqw(i,k) = 0. |
---|
| 1365 | ENDIF |
---|
| 1366 | ENDDO |
---|
| 1367 | ENDDO |
---|
| 1368 | c |
---|
| 1369 | C |
---|
[1146] | 1370 | c-------------Cstar computation--------------------------------- |
---|
| 1371 | DO i=1, klon |
---|
[1403] | 1372 | if (wk_adv(i)) then !!! nrlmd |
---|
[1146] | 1373 | sum_thu(i) = 0. |
---|
| 1374 | sum_tu(i) = 0. |
---|
| 1375 | sum_qu(i) = 0. |
---|
| 1376 | sum_thvu(i) = 0. |
---|
| 1377 | sum_dth(i) = 0. |
---|
| 1378 | sum_dq(i) = 0. |
---|
| 1379 | sum_rho(i) = 0. |
---|
| 1380 | sum_dtdwn(i) = 0. |
---|
| 1381 | sum_dqdwn(i) = 0. |
---|
| 1382 | |
---|
| 1383 | av_thu(i) = 0. |
---|
| 1384 | av_tu(i) =0. |
---|
| 1385 | av_qu(i) =0. |
---|
| 1386 | av_thvu(i) = 0. |
---|
| 1387 | av_dth(i) = 0. |
---|
| 1388 | av_dq(i) = 0. |
---|
| 1389 | av_rho(i) =0. |
---|
| 1390 | av_dtdwn(i) =0. |
---|
| 1391 | av_dqdwn(i) = 0. |
---|
[1403] | 1392 | end if |
---|
[1146] | 1393 | ENDDO |
---|
| 1394 | C |
---|
| 1395 | C Integrals (and wake top level number) |
---|
| 1396 | C -------------------------------------- |
---|
| 1397 | C |
---|
| 1398 | C Initialize sum_thvu to 1st level virt. pot. temp. |
---|
| 1399 | |
---|
| 1400 | DO i=1,klon |
---|
[1403] | 1401 | if (wk_adv(i)) then !!! nrlmd |
---|
[1146] | 1402 | z(i) = 1. |
---|
| 1403 | dz(i) = 1. |
---|
| 1404 | sum_thvu(i) = thu(i,1)*(1.+eps*qu(i,1))*dz(i) |
---|
| 1405 | sum_dth(i) = 0. |
---|
[1403] | 1406 | end if |
---|
[1146] | 1407 | ENDDO |
---|
| 1408 | |
---|
| 1409 | DO k = 1,klev |
---|
| 1410 | DO i=1,klon |
---|
[1403] | 1411 | if (wk_adv(i)) then !!! nrlmd |
---|
[1146] | 1412 | dz(i) = -(max(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) |
---|
| 1413 | IF (dz(i) .GT. 0) THEN |
---|
| 1414 | z(i) = z(i)+dz(i) |
---|
[1403] | 1415 | sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i) |
---|
| 1416 | sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i) |
---|
| 1417 | sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i) |
---|
| 1418 | sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i) |
---|
[1146] | 1419 | sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i) |
---|
| 1420 | sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i) |
---|
| 1421 | sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i) |
---|
| 1422 | sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i) |
---|
| 1423 | sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i) |
---|
| 1424 | ENDIF |
---|
[1403] | 1425 | end if |
---|
[1146] | 1426 | ENDDO |
---|
| 1427 | ENDDO |
---|
| 1428 | c |
---|
| 1429 | DO i=1,klon |
---|
[1403] | 1430 | if (wk_adv(i)) then !!! nrlmd |
---|
[1146] | 1431 | hw0(i) = z(i) |
---|
[1403] | 1432 | end if |
---|
[1146] | 1433 | ENDDO |
---|
| 1434 | c |
---|
| 1435 | C |
---|
| 1436 | C - WAPE and mean forcing computation |
---|
| 1437 | C --------------------------------------- |
---|
| 1438 | C |
---|
| 1439 | C --------------------------------------- |
---|
| 1440 | C |
---|
| 1441 | C Means |
---|
| 1442 | |
---|
| 1443 | DO i=1,klon |
---|
[1403] | 1444 | if (wk_adv(i)) then !!! nrlmd |
---|
[1146] | 1445 | av_thu(i) = sum_thu(i)/hw0(i) |
---|
| 1446 | av_tu(i) = sum_tu(i)/hw0(i) |
---|
| 1447 | av_qu(i) = sum_qu(i)/hw0(i) |
---|
| 1448 | av_thvu(i) = sum_thvu(i)/hw0(i) |
---|
| 1449 | av_dth(i) = sum_dth(i)/hw0(i) |
---|
| 1450 | av_dq(i) = sum_dq(i)/hw0(i) |
---|
| 1451 | av_rho(i) = sum_rho(i)/hw0(i) |
---|
| 1452 | av_dtdwn(i) = sum_dtdwn(i)/hw0(i) |
---|
| 1453 | av_dqdwn(i) = sum_dqdwn(i)/hw0(i) |
---|
| 1454 | c |
---|
| 1455 | wape(i) = - rg*hw0(i)*(av_dth(i) |
---|
| 1456 | $ + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)* |
---|
| 1457 | $ av_dq(i) ))/av_thvu(i) |
---|
[1403] | 1458 | end if |
---|
[1146] | 1459 | ENDDO |
---|
| 1460 | C |
---|
| 1461 | C Filter out bad wakes |
---|
| 1462 | |
---|
| 1463 | DO k = 1,klev |
---|
| 1464 | DO i=1,klon |
---|
[1403] | 1465 | if (wk_adv(i)) then !!! nrlmd |
---|
[1146] | 1466 | IF ( wape(i) .LT. 0.) THEN |
---|
| 1467 | deltatw(i,k) = 0. |
---|
| 1468 | deltaqw(i,k) = 0. |
---|
| 1469 | dth(i,k) = 0. |
---|
| 1470 | ENDIF |
---|
[1403] | 1471 | end if |
---|
[1146] | 1472 | ENDDO |
---|
| 1473 | ENDDO |
---|
| 1474 | c |
---|
| 1475 | DO i=1,klon |
---|
[1403] | 1476 | if (wk_adv(i)) then !!! nrlmd |
---|
[1146] | 1477 | IF ( wape(i) .LT. 0.) THEN |
---|
| 1478 | wape(i) = 0. |
---|
| 1479 | Cstar(i) = 0. |
---|
| 1480 | hw(i) = hwmin |
---|
| 1481 | sigmaw(i) = max(sigmad,sigd_con(i)) |
---|
| 1482 | fip(i) = 0. |
---|
| 1483 | gwake(i) = .FALSE. |
---|
| 1484 | ELSE |
---|
| 1485 | Cstar(i) = stark*sqrt(2.*wape(i)) |
---|
| 1486 | gwake(i) = .TRUE. |
---|
| 1487 | ENDIF |
---|
[1403] | 1488 | end if |
---|
[1146] | 1489 | ENDDO |
---|
| 1490 | |
---|
[974] | 1491 | ENDDO ! end sub-timestep loop |
---|
| 1492 | C |
---|
| 1493 | C ----------------------------------------------------------------- |
---|
| 1494 | c Get back to tendencies per second |
---|
| 1495 | c |
---|
| 1496 | DO k = 1,klev |
---|
| 1497 | DO i=1,klon |
---|
[1403] | 1498 | |
---|
| 1499 | ccc nrlmd IF ( wk_adv(i) .AND. k .LE. kupper(i)) THEN |
---|
| 1500 | IF ( OK_qx_qw(i) .AND. k .LE. kupper(i)) THEN |
---|
| 1501 | ccc |
---|
[974] | 1502 | dtls(i,k) = dtls(i,k)/dtime |
---|
| 1503 | dqls(i,k) = dqls(i,k)/dtime |
---|
| 1504 | d_deltatw2(i,k)=d_deltatw2(i,k)/dtime |
---|
| 1505 | d_deltaqw2(i,k)=d_deltaqw2(i,k)/dtime |
---|
| 1506 | d_deltat_gw(i,k) = d_deltat_gw(i,k)/dtime |
---|
[1403] | 1507 | cc print*,'k,dqls,omg,entr,detr',k,dqls(i,k),omg(i,k),entr(i,k) |
---|
| 1508 | cc $ ,death_rate(i)*sigmaw(i) |
---|
[974] | 1509 | ENDIF |
---|
| 1510 | ENDDO |
---|
| 1511 | ENDDO |
---|
[1403] | 1512 | |
---|
[974] | 1513 | c |
---|
| 1514 | c---------------------------------------------------------- |
---|
| 1515 | c Determine wake final state; recompute wape, cstar, ktop; |
---|
| 1516 | c filter out bad wakes. |
---|
| 1517 | c---------------------------------------------------------- |
---|
| 1518 | c |
---|
| 1519 | C 2.1 - Undisturbed area and Wake integrals |
---|
| 1520 | C --------------------------------------------------------- |
---|
| 1521 | |
---|
| 1522 | DO i=1,klon |
---|
[1403] | 1523 | ccc nrlmd if (wk_adv(i)) then !!! nrlmd |
---|
| 1524 | if (OK_qx_qw(i)) then |
---|
| 1525 | ccc |
---|
[974] | 1526 | z(i) = 0. |
---|
| 1527 | sum_thu(i) = 0. |
---|
| 1528 | sum_tu(i) = 0. |
---|
| 1529 | sum_qu(i) = 0. |
---|
| 1530 | sum_thvu(i) = 0. |
---|
| 1531 | sum_dth(i) = 0. |
---|
| 1532 | sum_dq(i) = 0. |
---|
| 1533 | sum_rho(i) = 0. |
---|
| 1534 | sum_dtdwn(i) = 0. |
---|
| 1535 | sum_dqdwn(i) = 0. |
---|
| 1536 | |
---|
| 1537 | av_thu(i) = 0. |
---|
| 1538 | av_tu(i) =0. |
---|
| 1539 | av_qu(i) =0. |
---|
| 1540 | av_thvu(i) = 0. |
---|
| 1541 | av_dth(i) = 0. |
---|
| 1542 | av_dq(i) = 0. |
---|
| 1543 | av_rho(i) =0. |
---|
| 1544 | av_dtdwn(i) =0. |
---|
| 1545 | av_dqdwn(i) = 0. |
---|
[1403] | 1546 | end if |
---|
[974] | 1547 | ENDDO |
---|
| 1548 | C Potential temperatures and humidity |
---|
| 1549 | c---------------------------------------------------------- |
---|
| 1550 | |
---|
| 1551 | DO k =1,klev |
---|
| 1552 | DO i=1,klon |
---|
[1403] | 1553 | ccc nrlmd IF ( wk_adv(i)) THEN |
---|
| 1554 | if (OK_qx_qw(i)) then |
---|
| 1555 | ccc |
---|
[974] | 1556 | rho(i,k) = p(i,k)/(rd*te(i,k)) |
---|
| 1557 | IF(k .eq. 1) THEN |
---|
| 1558 | rhoh(i,k) = ph(i,k)/(rd*te(i,k)) |
---|
| 1559 | zhh(i,k)=0 |
---|
| 1560 | ELSE |
---|
| 1561 | rhoh(i,k) = ph(i,k)*2./(rd*(te(i,k)+te(i,k-1))) |
---|
| 1562 | zhh(i,k)=(ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG)+zhh(i,k-1) |
---|
| 1563 | ENDIF |
---|
| 1564 | the(i,k) = te(i,k)/ppi(i,k) |
---|
| 1565 | thu(i,k) = (te(i,k) - deltatw(i,k)*sigmaw(i))/ppi(i,k) |
---|
| 1566 | tu(i,k) = te(i,k) - deltatw(i,k)*sigmaw(i) |
---|
| 1567 | qu(i,k) = qe(i,k) - deltaqw(i,k)*sigmaw(i) |
---|
[1403] | 1568 | rhow(i,k) = p(i,k)/(rd*(te(i,k)+deltatw(i,k))) |
---|
[974] | 1569 | dth(i,k) = deltatw(i,k)/ppi(i,k) |
---|
[1146] | 1570 | ENDIF |
---|
[974] | 1571 | ENDDO |
---|
| 1572 | ENDDO |
---|
| 1573 | |
---|
| 1574 | C Integrals (and wake top level number) |
---|
| 1575 | C ----------------------------------------------------------- |
---|
| 1576 | |
---|
| 1577 | C Initialize sum_thvu to 1st level virt. pot. temp. |
---|
| 1578 | |
---|
| 1579 | DO i=1,klon |
---|
[1403] | 1580 | ccc nrlmd IF ( wk_adv(i)) THEN |
---|
| 1581 | if (OK_qx_qw(i)) then |
---|
| 1582 | ccc |
---|
[974] | 1583 | z(i) = 1. |
---|
| 1584 | dz(i) = 1. |
---|
| 1585 | sum_thvu(i) = thu(i,1)*(1.+eps*qu(i,1))*dz(i) |
---|
| 1586 | sum_dth(i) = 0. |
---|
[1146] | 1587 | ENDIF |
---|
[974] | 1588 | ENDDO |
---|
| 1589 | |
---|
| 1590 | DO k = 1,klev |
---|
| 1591 | DO i=1,klon |
---|
[1403] | 1592 | ccc nrlmd IF ( wk_adv(i)) THEN |
---|
| 1593 | if (OK_qx_qw(i)) then |
---|
| 1594 | ccc |
---|
[974] | 1595 | dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*rg) |
---|
| 1596 | IF (dz(i) .GT. 0) THEN |
---|
| 1597 | z(i) = z(i)+dz(i) |
---|
| 1598 | sum_thu(i) = sum_thu(i) + thu(i,k)*dz(i) |
---|
| 1599 | sum_tu(i) = sum_tu(i) + tu(i,k)*dz(i) |
---|
| 1600 | sum_qu(i) = sum_qu(i) + qu(i,k)*dz(i) |
---|
| 1601 | sum_thvu(i) = sum_thvu(i) + thu(i,k)*(1.+eps*qu(i,k))*dz(i) |
---|
| 1602 | sum_dth(i) = sum_dth(i) + dth(i,k)*dz(i) |
---|
| 1603 | sum_dq(i) = sum_dq(i) + deltaqw(i,k)*dz(i) |
---|
| 1604 | sum_rho(i) = sum_rho(i) + rhow(i,k)*dz(i) |
---|
| 1605 | sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i,k)*dz(i) |
---|
| 1606 | sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i,k)*dz(i) |
---|
| 1607 | ENDIF |
---|
[1146] | 1608 | ENDIF |
---|
[974] | 1609 | ENDDO |
---|
| 1610 | ENDDO |
---|
| 1611 | c |
---|
| 1612 | DO i=1,klon |
---|
[1403] | 1613 | ccc nrlmd IF ( wk_adv(i)) THEN |
---|
| 1614 | if (OK_qx_qw(i)) then |
---|
| 1615 | ccc |
---|
[974] | 1616 | hw0(i) = z(i) |
---|
[1146] | 1617 | ENDIF |
---|
[974] | 1618 | ENDDO |
---|
| 1619 | c |
---|
[1146] | 1620 | C - WAPE and mean forcing computation |
---|
[974] | 1621 | C------------------------------------------------------------- |
---|
| 1622 | |
---|
| 1623 | C Means |
---|
| 1624 | |
---|
| 1625 | DO i=1, klon |
---|
[1403] | 1626 | ccc nrlmd IF ( wk_adv(i)) THEN |
---|
| 1627 | if (OK_qx_qw(i)) then |
---|
| 1628 | ccc |
---|
[974] | 1629 | av_thu(i) = sum_thu(i)/hw0(i) |
---|
| 1630 | av_tu(i) = sum_tu(i)/hw0(i) |
---|
| 1631 | av_qu(i) = sum_qu(i)/hw0(i) |
---|
| 1632 | av_thvu(i) = sum_thvu(i)/hw0(i) |
---|
| 1633 | av_dth(i) = sum_dth(i)/hw0(i) |
---|
| 1634 | av_dq(i) = sum_dq(i)/hw0(i) |
---|
| 1635 | av_rho(i) = sum_rho(i)/hw0(i) |
---|
| 1636 | av_dtdwn(i) = sum_dtdwn(i)/hw0(i) |
---|
| 1637 | av_dqdwn(i) = sum_dqdwn(i)/hw0(i) |
---|
| 1638 | |
---|
| 1639 | wape2(i) = - rg*hw0(i)*(av_dth(i) |
---|
[1403] | 1640 | $ + eps*(av_thu(i)*av_dq(i)+av_dth(i)*av_qu(i) |
---|
| 1641 | $ + av_dth(i)*av_dq(i) ))/av_thvu(i) |
---|
[1146] | 1642 | ENDIF |
---|
[974] | 1643 | ENDDO |
---|
| 1644 | |
---|
[1146] | 1645 | C Prognostic variable update |
---|
[974] | 1646 | C ------------------------------------------------------------ |
---|
| 1647 | |
---|
| 1648 | C Filter out bad wakes |
---|
| 1649 | c |
---|
| 1650 | DO k = 1,klev |
---|
| 1651 | DO i=1,klon |
---|
[1403] | 1652 | ccc nrlmd IF ( wk_adv(i) .AND. wape2(i) .LT. 0.) THEN |
---|
| 1653 | if (OK_qx_qw(i) .AND. wape2(i) .LT. 0.) then |
---|
| 1654 | ccc |
---|
[974] | 1655 | deltatw(i,k) = 0. |
---|
| 1656 | deltaqw(i,k) = 0. |
---|
| 1657 | dth(i,k) = 0. |
---|
| 1658 | ENDIF |
---|
| 1659 | ENDDO |
---|
| 1660 | ENDDO |
---|
| 1661 | c |
---|
| 1662 | |
---|
| 1663 | DO i=1, klon |
---|
[1403] | 1664 | ccc nrlmd IF ( wk_adv(i)) THEN |
---|
| 1665 | if (OK_qx_qw(i)) then |
---|
| 1666 | ccc |
---|
[1146] | 1667 | IF ( wape2(i) .LT. 0.) THEN |
---|
[974] | 1668 | wape2(i) = 0. |
---|
| 1669 | Cstar2(i) = 0. |
---|
| 1670 | hw(i) = hwmin |
---|
| 1671 | sigmaw(i) = amax1(sigmad,sigd_con(i)) |
---|
| 1672 | fip(i) = 0. |
---|
| 1673 | gwake(i) = .FALSE. |
---|
| 1674 | ELSE |
---|
| 1675 | if(prt_level.ge.10) print*,'wape2>0' |
---|
| 1676 | Cstar2(i) = stark*sqrt(2.*wape2(i)) |
---|
| 1677 | gwake(i) = .TRUE. |
---|
| 1678 | ENDIF |
---|
[1146] | 1679 | ENDIF |
---|
[974] | 1680 | ENDDO |
---|
| 1681 | c |
---|
| 1682 | DO i=1, klon |
---|
[1403] | 1683 | ccc nrlmd IF ( wk_adv(i)) THEN |
---|
| 1684 | if (OK_qx_qw(i)) then |
---|
| 1685 | ccc |
---|
[974] | 1686 | ktopw(i) = ktop(i) |
---|
[1146] | 1687 | ENDIF |
---|
[974] | 1688 | ENDDO |
---|
| 1689 | c |
---|
| 1690 | DO i=1, klon |
---|
[1403] | 1691 | ccc nrlmd IF ( wk_adv(i)) THEN |
---|
| 1692 | if (OK_qx_qw(i)) then |
---|
| 1693 | ccc |
---|
[1146] | 1694 | IF (ktopw(i) .gt. 0 .and. gwake(i)) then |
---|
[974] | 1695 | |
---|
| 1696 | Cjyg1 Utilisation d'un h_efficace constant ( ~ feeding layer) |
---|
| 1697 | ccc heff = 600. |
---|
| 1698 | C Utilisation de la hauteur hw |
---|
| 1699 | cc heff = 0.7*hw |
---|
| 1700 | heff(i) = hw(i) |
---|
| 1701 | |
---|
| 1702 | FIP(i) = 0.5*rho(i,ktopw(i))*Cstar2(i)**3*heff(i)*2* |
---|
[1403] | 1703 | $ sqrt(sigmaw(i)*wdens(i)*3.14) |
---|
[974] | 1704 | FIP(i) = alpk * FIP(i) |
---|
| 1705 | Cjyg2 |
---|
| 1706 | ELSE |
---|
| 1707 | FIP(i) = 0. |
---|
| 1708 | ENDIF |
---|
[1146] | 1709 | ENDIF |
---|
[974] | 1710 | ENDDO |
---|
| 1711 | c |
---|
| 1712 | C Limitation de sigmaw |
---|
[1403] | 1713 | |
---|
| 1714 | ccc nrlmd |
---|
| 1715 | c DO i=1,klon |
---|
| 1716 | c IF (OK_qx_qw(i)) THEN |
---|
| 1717 | c IF (sigmaw(i).GE.sigmaw_max) sigmaw(i)=sigmaw_max |
---|
| 1718 | c ENDIF |
---|
| 1719 | c ENDDO |
---|
| 1720 | ccc |
---|
[974] | 1721 | DO k = 1,klev |
---|
| 1722 | DO i=1, klon |
---|
[1403] | 1723 | |
---|
| 1724 | ccc nrlmd On maintient désormais constant sigmaw en régime permanent |
---|
| 1725 | ccc IF ((sigmaw(i).GT.sigmaw_max).or. |
---|
| 1726 | IF ( ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or. |
---|
| 1727 | $ (ktopw(i).le.2) .OR. |
---|
| 1728 | $ .not. OK_qx_qw(i) ) THEN |
---|
| 1729 | ccc |
---|
[974] | 1730 | dtls(i,k) = 0. |
---|
| 1731 | dqls(i,k) = 0. |
---|
| 1732 | deltatw(i,k) = 0. |
---|
| 1733 | deltaqw(i,k) = 0. |
---|
| 1734 | ENDIF |
---|
| 1735 | ENDDO |
---|
| 1736 | ENDDO |
---|
| 1737 | c |
---|
[1403] | 1738 | ccc nrlmd On maintient désormais constant sigmaw en régime permanent |
---|
[974] | 1739 | DO i=1, klon |
---|
[1403] | 1740 | IF ( ((wape(i).ge.wape2(i)).and.(wape2(i).le.1.0)).or. |
---|
| 1741 | $ (ktopw(i).le.2) .OR. |
---|
| 1742 | $ .not. OK_qx_qw(i) ) THEN |
---|
[974] | 1743 | wape(i) = 0. |
---|
[1403] | 1744 | cstar(i)=0. |
---|
[974] | 1745 | hw(i) = hwmin |
---|
| 1746 | sigmaw(i) = sigmad |
---|
| 1747 | fip(i) = 0. |
---|
| 1748 | ELSE |
---|
| 1749 | wape(i) = wape2(i) |
---|
[1403] | 1750 | cstar(i)=cstar2(i) |
---|
[974] | 1751 | ENDIF |
---|
[1403] | 1752 | cc print*,'wape wape2 ktopw OK_qx_qw =', |
---|
| 1753 | cc $ wape(i),wape2(i),ktopw(i),OK_qx_qw(i) |
---|
[974] | 1754 | ENDDO |
---|
| 1755 | c |
---|
| 1756 | c |
---|
| 1757 | RETURN |
---|
| 1758 | END |
---|
[1146] | 1759 | |
---|
[1403] | 1760 | SUBROUTINE wake_vec_modulation(nlon,nl,wk_adv,epsilon,qe,d_qe, |
---|
[1146] | 1761 | $ deltaqw,d_deltaqw,sigmaw,d_sigmaw,alpha) |
---|
| 1762 | c------------------------------------------------------ |
---|
| 1763 | cDtermination du coefficient alpha tel que les tendances |
---|
| 1764 | c corriges alpha*d_G, pour toutes les grandeurs G, correspondent |
---|
| 1765 | c a une humidite positive dans la zone (x) et dans la zone (w). |
---|
| 1766 | c------------------------------------------------------ |
---|
| 1767 | c |
---|
| 1768 | |
---|
| 1769 | c Input |
---|
| 1770 | REAL qe(nlon,nl),d_qe(nlon,nl) |
---|
| 1771 | REAL deltaqw(nlon,nl),d_deltaqw(nlon,nl) |
---|
| 1772 | REAL sigmaw(nlon),d_sigmaw(nlon) |
---|
| 1773 | LOGICAL wk_adv(nlon) |
---|
| 1774 | INTEGER nl,nlon |
---|
| 1775 | c Output |
---|
| 1776 | REAL alpha(nlon) |
---|
| 1777 | c Internal variables |
---|
[1403] | 1778 | REAL zeta(nlon,nl) |
---|
[1146] | 1779 | REAL alpha1(nlon) |
---|
[1403] | 1780 | REAL x,a,b,c,discrim |
---|
[1146] | 1781 | REAL epsilon |
---|
[1403] | 1782 | ! DATA epsilon/1.e-15/ |
---|
[1146] | 1783 | c |
---|
| 1784 | DO k=1,nl |
---|
| 1785 | DO i = 1,nlon |
---|
| 1786 | IF (wk_adv(i)) THEN |
---|
| 1787 | IF ((deltaqw(i,k)+d_deltaqw(i,k)).ge.0.) then |
---|
[1403] | 1788 | zeta(i,k)=0. |
---|
[1146] | 1789 | ELSE |
---|
[1403] | 1790 | zeta(i,k)=1. |
---|
[1146] | 1791 | END IF |
---|
| 1792 | ENDIF |
---|
| 1793 | ENDDO |
---|
| 1794 | DO i = 1,nlon |
---|
| 1795 | IF (wk_adv(i)) THEN |
---|
[1403] | 1796 | x = qe(i,k)+(zeta(i,k)-sigmaw(i))*deltaqw(i,k) |
---|
| 1797 | $ + d_qe(i,k)+(zeta(i,k)-sigmaw(i))*d_deltaqw(i,k) |
---|
| 1798 | $ - d_sigmaw(i)*(deltaqw(i,k)+d_deltaqw(i,k)) |
---|
| 1799 | a = -d_sigmaw(i)*d_deltaqw(i,k) |
---|
| 1800 | b = d_qe(i,k)+(zeta(i,k)-sigmaw(i))*d_deltaqw(i,k) |
---|
| 1801 | $ - deltaqw(i,k)*d_sigmaw(i) |
---|
| 1802 | c = qe(i,k)+(zeta(i,k)-sigmaw(i))*deltaqw(i,k)+epsilon |
---|
| 1803 | discrim = b*b-4.*a*c |
---|
| 1804 | c print*, 'x, a, b, c, discrim', x, a, b, c, discrim |
---|
| 1805 | IF (a+b .GE. 0.) THEN !! Condition suffisante pour la positivité de ovap |
---|
[1146] | 1806 | alpha1(i)=1. |
---|
| 1807 | ELSE |
---|
| 1808 | IF (x .GE. 0.) THEN |
---|
| 1809 | alpha1(i)=1. |
---|
| 1810 | ELSE |
---|
| 1811 | IF (a .GT. 0.) THEN |
---|
| 1812 | alpha1(i)=0.9*min( (2.*c)/(-b+sqrt(discrim)), |
---|
| 1813 | $ (-b+sqrt(discrim))/(2.*a) ) |
---|
[1403] | 1814 | ELSE IF (a .eq. 0.) then |
---|
[1146] | 1815 | alpha1(i)=0.9*(-c/b) |
---|
| 1816 | ELSE |
---|
[1403] | 1817 | c print*,'a,b,c discrim',a,b,c discrim |
---|
[1146] | 1818 | alpha1(i)=0.9*max( (2.*c)/(-b+sqrt(discrim)), |
---|
| 1819 | $ (-b+sqrt(discrim))/(2.*a) ) |
---|
| 1820 | ENDIF |
---|
| 1821 | ENDIF |
---|
| 1822 | ENDIF |
---|
[1403] | 1823 | alpha(i) = min(alpha(i),alpha1(i)) |
---|
[1146] | 1824 | ENDIF |
---|
| 1825 | ENDDO |
---|
| 1826 | ENDDO |
---|
[1403] | 1827 | ! |
---|
[1146] | 1828 | return |
---|
| 1829 | end |
---|
| 1830 | |
---|
[974] | 1831 | Subroutine WAKE_scal (p,ph,ppi,dtime,sigd_con |
---|
| 1832 | : ,te0,qe0,omgb |
---|
| 1833 | : ,dtdwn,dqdwn,amdwn,amup,dta,dqa |
---|
| 1834 | : ,wdtPBL,wdqPBL,udtPBL,udqPBL |
---|
| 1835 | o ,deltatw,deltaqw,dth,hw,sigmaw,wape,fip,gfl |
---|
| 1836 | o ,dtls,dqls |
---|
| 1837 | o ,ktopw,omgbdth,dp_omgb,wdens |
---|
| 1838 | o ,tu,qu |
---|
| 1839 | o ,dtKE,dqKE |
---|
| 1840 | o ,dtPBL,dqPBL |
---|
| 1841 | o ,omg,dp_deltomg,spread |
---|
| 1842 | o ,Cstar,d_deltat_gw |
---|
| 1843 | o ,d_deltatw2,d_deltaqw2) |
---|
| 1844 | |
---|
| 1845 | *************************************************************** |
---|
| 1846 | * * |
---|
| 1847 | * WAKE * |
---|
| 1848 | * retour a un Pupper fixe * |
---|
| 1849 | * * |
---|
| 1850 | * written by : GRANDPEIX Jean-Yves 09/03/2000 * |
---|
| 1851 | * modified by : ROEHRIG Romain 01/29/2007 * |
---|
| 1852 | *************************************************************** |
---|
| 1853 | c |
---|
[940] | 1854 | USE dimphy |
---|
[879] | 1855 | IMPLICIT none |
---|
| 1856 | c============================================================================ |
---|
| 1857 | C |
---|
| 1858 | C |
---|
| 1859 | C But : Decrire le comportement des poches froides apparaissant dans les |
---|
| 1860 | C grands systemes convectifs, et fournir l'energie disponible pour |
---|
| 1861 | C le declenchement de nouvelles colonnes convectives. |
---|
| 1862 | C |
---|
| 1863 | C Variables d'etat : deltatw : ecart de temperature wake-undisturbed area |
---|
| 1864 | C deltaqw : ecart d'humidite wake-undisturbed area |
---|
| 1865 | C sigmaw : fraction d'aire occupee par la poche. |
---|
| 1866 | C |
---|
| 1867 | C Variable de sortie : |
---|
| 1868 | c |
---|
| 1869 | c wape : WAke Potential Energy |
---|
| 1870 | c fip : Front Incident Power (W/m2) - ALP |
---|
| 1871 | c gfl : Gust Front Length per unit area (m-1) |
---|
| 1872 | C dtls : large scale temperature tendency due to wake |
---|
| 1873 | C dqls : large scale humidity tendency due to wake |
---|
| 1874 | C hw : hauteur de la poche |
---|
| 1875 | C dp_omgb : vertical gradient of large scale omega |
---|
| 1876 | C omgbdth: flux of Delta_Theta transported by LS omega |
---|
| 1877 | C dtKE : differential heating (wake - unpertubed) |
---|
| 1878 | C dqKE : differential moistening (wake - unpertubed) |
---|
| 1879 | C omg : Delta_omg =vertical velocity diff. wake-undist. (Pa/s) |
---|
| 1880 | C dp_deltomg : vertical gradient of omg (s-1) |
---|
| 1881 | C spread : spreading term in dt_wake and dq_wake |
---|
| 1882 | C deltatw : updated temperature difference (T_w-T_u). |
---|
| 1883 | C deltaqw : updated humidity difference (q_w-q_u). |
---|
| 1884 | C sigmaw : updated wake fractional area. |
---|
| 1885 | C d_deltat_gw : delta T tendency due to GW |
---|
| 1886 | c |
---|
| 1887 | C Variables d'entree : |
---|
| 1888 | c |
---|
| 1889 | c aire : aire de la maille |
---|
| 1890 | c te0 : temperature dans l'environnement (K) |
---|
| 1891 | C qe0 : humidite dans l'environnement (kg/kg) |
---|
| 1892 | C omgb : vitesse verticale moyenne sur la maille (Pa/s) |
---|
| 1893 | C dtdwn: source de chaleur due aux descentes (K/s) |
---|
| 1894 | C dqdwn: source d'humidite due aux descentes (kg/kg/s) |
---|
| 1895 | C dta : source de chaleur due courants satures et detrain (K/s) |
---|
| 1896 | C dqa : source d'humidite due aux courants satures et detra (kg/kg/s) |
---|
| 1897 | C amdwn: flux de masse total des descentes, par unite de |
---|
| 1898 | C surface de la maille (kg/m2/s) |
---|
| 1899 | C amup : flux de masse total des ascendances, par unite de |
---|
| 1900 | C surface de la maille (kg/m2/s) |
---|
| 1901 | C p : pressions aux milieux des couches (Pa) |
---|
| 1902 | C ph : pressions aux interfaces (Pa) |
---|
| 1903 | C ppi : (p/p_0)**kapa (adim) |
---|
| 1904 | C dtime: increment temporel (s) |
---|
| 1905 | c |
---|
| 1906 | C Variables internes : |
---|
| 1907 | c |
---|
| 1908 | c rhow : masse volumique de la poche froide |
---|
| 1909 | C rho : environment density at P levels |
---|
| 1910 | C rhoh : environment density at Ph levels |
---|
| 1911 | C te : environment temperature | may change within |
---|
| 1912 | C qe : environment humidity | sub-time-stepping |
---|
| 1913 | C the : environment potential temperature |
---|
| 1914 | C thu : potential temperature in undisturbed area |
---|
| 1915 | C tu : temperature in undisturbed area |
---|
| 1916 | C qu : humidity in undisturbed area |
---|
| 1917 | C dp_omgb: vertical gradient og LS omega |
---|
| 1918 | C omgbw : wake average vertical omega |
---|
| 1919 | C dp_omgbw: vertical gradient of omgbw |
---|
| 1920 | C omgbdq : flux of Delta_q transported by LS omega |
---|
| 1921 | C dth : potential temperature diff. wake-undist. |
---|
| 1922 | C th1 : first pot. temp. for vertical advection (=thu) |
---|
| 1923 | C th2 : second pot. temp. for vertical advection (=thw) |
---|
| 1924 | C q1 : first humidity for vertical advection |
---|
| 1925 | C q2 : second humidity for vertical advection |
---|
| 1926 | C d_deltatw : terme de redistribution pour deltatw |
---|
| 1927 | C d_deltaqw : terme de redistribution pour deltaqw |
---|
| 1928 | C deltatw0 : deltatw initial |
---|
| 1929 | C deltaqw0 : deltaqw initial |
---|
| 1930 | C hw0 : hw initial |
---|
| 1931 | C sigmaw0: sigmaw initial |
---|
| 1932 | C amflux : horizontal mass flux through wake boundary |
---|
| 1933 | C wdens : number of wakes per unit area (3D) or per |
---|
| 1934 | C unit length (2D) |
---|
| 1935 | C Tgw : 1 sur la période de onde de gravité |
---|
| 1936 | c Cgw : vitesse de propagation de onde de gravité |
---|
| 1937 | c LL : distance entre 2 poches |
---|
| 1938 | |
---|
| 1939 | c------------------------------------------------------------------------- |
---|
| 1940 | c Déclaration de variables |
---|
| 1941 | c------------------------------------------------------------------------- |
---|
| 1942 | |
---|
| 1943 | #include "dimensions.h" |
---|
[940] | 1944 | cccc#include "dimphy.h" |
---|
[879] | 1945 | #include "YOMCST.h" |
---|
| 1946 | #include "cvthermo.h" |
---|
[953] | 1947 | #include "iniprint.h" |
---|
[879] | 1948 | |
---|
| 1949 | c Arguments en entree |
---|
| 1950 | c-------------------- |
---|
| 1951 | |
---|
| 1952 | REAL p(klev),ph(klev+1),ppi(klev) |
---|
| 1953 | REAL dtime |
---|
| 1954 | REAL te0(klev),qe0(klev) |
---|
| 1955 | REAL omgb(klev+1) |
---|
| 1956 | REAL dtdwn(klev), dqdwn(klev) |
---|
| 1957 | REAL wdtPBL(klev),wdqPBL(klev) |
---|
| 1958 | REAL udtPBL(klev),udqPBL(klev) |
---|
| 1959 | REAL amdwn(klev), amup(klev) |
---|
| 1960 | REAL dta(klev), dqa(klev) |
---|
| 1961 | REAL sigd_con |
---|
| 1962 | |
---|
| 1963 | c Sorties |
---|
| 1964 | c-------- |
---|
| 1965 | |
---|
| 1966 | REAL deltatw(klev), deltaqw(klev), dth(klev) |
---|
| 1967 | REAL tu(klev), qu(klev) |
---|
| 1968 | REAL dtls(klev), dqls(klev) |
---|
| 1969 | REAL dtKE(klev), dqKE(klev) |
---|
| 1970 | REAL dtPBL(klev), dqPBL(klev) |
---|
| 1971 | REAL spread(klev) |
---|
| 1972 | REAL d_deltatgw(klev) |
---|
| 1973 | REAL d_deltatw2(klev), d_deltaqw2(klev) |
---|
| 1974 | REAL omgbdth(klev+1), omg(klev+1) |
---|
| 1975 | REAL dp_omgb(klev), dp_deltomg(klev) |
---|
| 1976 | REAL d_deltat_gw(klev) |
---|
| 1977 | REAL hw, sigmaw, wape, fip, gfl, Cstar |
---|
| 1978 | INTEGER ktopw |
---|
| 1979 | |
---|
| 1980 | c Variables internes |
---|
| 1981 | c------------------- |
---|
| 1982 | |
---|
| 1983 | c Variables à fixer |
---|
| 1984 | REAL ALON |
---|
| 1985 | REAL coefgw |
---|
| 1986 | REAL wdens0, wdens |
---|
| 1987 | REAL stark |
---|
| 1988 | REAL alpk |
---|
| 1989 | REAL delta_t_min |
---|
| 1990 | REAL Pupper |
---|
| 1991 | INTEGER nsub |
---|
| 1992 | REAL dtimesub |
---|
| 1993 | REAL sigmad, hwmin |
---|
| 1994 | |
---|
| 1995 | c Variables de sauvegarde |
---|
| 1996 | REAL deltatw0(klev) |
---|
| 1997 | REAL deltaqw0(klev) |
---|
| 1998 | REAL te(klev), qe(klev) |
---|
| 1999 | REAL sigmaw0, sigmaw1 |
---|
| 2000 | |
---|
| 2001 | c Variables pour les GW |
---|
| 2002 | REAL LL |
---|
| 2003 | REAL N2(klev) |
---|
| 2004 | REAL Cgw(klev) |
---|
| 2005 | REAL Tgw(klev) |
---|
| 2006 | |
---|
| 2007 | c Variables liées au calcul de hw |
---|
| 2008 | REAL ptop_provis, ptop, ptop_new |
---|
| 2009 | REAL sum_dth |
---|
| 2010 | REAL dthmin |
---|
| 2011 | REAL z, dz, hw0 |
---|
| 2012 | INTEGER ktop, kupper |
---|
| 2013 | |
---|
| 2014 | c Autres variables internes |
---|
| 2015 | INTEGER isubstep, k |
---|
| 2016 | |
---|
| 2017 | REAL sum_thu, sum_tu, sum_qu,sum_thvu |
---|
| 2018 | REAL sum_dq, sum_rho |
---|
| 2019 | REAL sum_dtdwn, sum_dqdwn |
---|
| 2020 | REAL av_thu, av_tu, av_qu, av_thvu |
---|
| 2021 | REAL av_dth, av_dq, av_rho |
---|
| 2022 | REAL av_dtdwn, av_dqdwn |
---|
| 2023 | |
---|
| 2024 | REAL rho(klev), rhoh(klev+1), rhow(klev) |
---|
| 2025 | REAL rhow_moyen(klev) |
---|
| 2026 | REAL zh(klev), zhh(klev+1) |
---|
| 2027 | REAL epaisseur1(klev), epaisseur2(klev) |
---|
| 2028 | |
---|
| 2029 | REAL the(klev), thu(klev) |
---|
| 2030 | |
---|
| 2031 | REAL d_deltatw(klev), d_deltaqw(klev) |
---|
| 2032 | |
---|
| 2033 | REAL omgbw(klev+1), omgtop |
---|
| 2034 | REAL dp_omgbw(klev) |
---|
| 2035 | REAL ztop, dztop |
---|
| 2036 | REAL alpha_up(klev) |
---|
| 2037 | |
---|
| 2038 | REAL RRe1, RRe2, RRd1, RRd2 |
---|
| 2039 | REAL Th1(klev), Th2(klev), q1(klev), q2(klev) |
---|
| 2040 | REAL D_Th1(klev), D_Th2(klev), D_dth(klev) |
---|
| 2041 | REAL D_q1(klev), D_q2(klev), D_dq(klev) |
---|
| 2042 | REAL omgbdq(klev) |
---|
| 2043 | |
---|
| 2044 | REAL ff, gg |
---|
| 2045 | REAL wape2, Cstar2, heff |
---|
| 2046 | |
---|
| 2047 | REAL Crep(klev) |
---|
| 2048 | REAL Crep_upper, Crep_sol |
---|
| 2049 | |
---|
| 2050 | C------------------------------------------------------------------------- |
---|
| 2051 | c Initialisations |
---|
| 2052 | c------------------------------------------------------------------------- |
---|
| 2053 | |
---|
| 2054 | c print*, 'wake initialisations' |
---|
| 2055 | |
---|
| 2056 | c Essais d'initialisation avec sigmaw = 0.02 et hw = 10. |
---|
| 2057 | c------------------------------------------------------------------------- |
---|
| 2058 | |
---|
| 2059 | DATA sigmad, hwmin /.02,10./ |
---|
| 2060 | |
---|
| 2061 | C Longueur de maille (en m) |
---|
| 2062 | c------------------------------------------------------------------------- |
---|
| 2063 | |
---|
| 2064 | c ALON = 3.e5 |
---|
| 2065 | ALON = 1.e6 |
---|
| 2066 | |
---|
| 2067 | |
---|
| 2068 | C Configuration de coefgw,stark,wdens (22/02/06 by YU Jingmei) |
---|
| 2069 | c |
---|
| 2070 | c coefgw : Coefficient pour les ondes de gravité |
---|
| 2071 | c stark : Coefficient k dans Cstar=k*sqrt(2*WAPE) |
---|
| 2072 | c wdens : Densité de poche froide par maille |
---|
| 2073 | c------------------------------------------------------------------------- |
---|
| 2074 | |
---|
| 2075 | coefgw=10 |
---|
| 2076 | c coefgw=1 |
---|
| 2077 | c wdens0 = 1.0/(alon**2) |
---|
| 2078 | wdens = 1.0/(alon**2) |
---|
| 2079 | stark = 0.50 |
---|
| 2080 | cCRtest |
---|
| 2081 | alpk=0.1 |
---|
| 2082 | c alpk = 1.0 |
---|
| 2083 | c alpk = 0.5 |
---|
| 2084 | c alpk = 0.05 |
---|
| 2085 | Crep_upper=0.9 |
---|
| 2086 | Crep_sol=1.0 |
---|
| 2087 | |
---|
| 2088 | |
---|
| 2089 | C Minimum value for |T_wake - T_undist|. Used for wake top definition |
---|
| 2090 | c------------------------------------------------------------------------- |
---|
| 2091 | |
---|
| 2092 | delta_t_min = 0.2 |
---|
| 2093 | |
---|
| 2094 | |
---|
| 2095 | C 1. - Save initial values and initialize tendencies |
---|
| 2096 | C -------------------------------------------------- |
---|
| 2097 | |
---|
| 2098 | DO k=1,klev |
---|
| 2099 | deltatw0(k) = deltatw(k) |
---|
| 2100 | deltaqw0(k)= deltaqw(k) |
---|
| 2101 | te(k) = te0(k) |
---|
| 2102 | qe(k) = qe0(k) |
---|
| 2103 | dtls(k) = 0. |
---|
| 2104 | dqls(k) = 0. |
---|
| 2105 | d_deltat_gw(k)=0. |
---|
| 2106 | d_deltatw2(k)=0. |
---|
| 2107 | d_deltaqw2(k)=0. |
---|
| 2108 | ENDDO |
---|
| 2109 | c sigmaw1=sigmaw |
---|
| 2110 | c IF (sigd_con.GT.sigmaw1) THEN |
---|
| 2111 | c print*, 'sigmaw,sigd_con', sigmaw, sigd_con |
---|
| 2112 | c ENDIF |
---|
| 2113 | sigmaw = max(sigmaw,sigd_con) |
---|
| 2114 | sigmaw = max(sigmaw,sigmad) |
---|
| 2115 | sigmaw = min(sigmaw,0.99) |
---|
| 2116 | sigmaw0 = sigmaw |
---|
| 2117 | c wdens=wdens0/(10.*sigmaw) |
---|
| 2118 | c IF (sigd_con.GT.sigmaw1) THEN |
---|
| 2119 | c print*, 'sigmaw1,sigd1', sigmaw, sigd_con |
---|
| 2120 | c ENDIF |
---|
| 2121 | |
---|
| 2122 | C 2. - Prognostic part |
---|
| 2123 | C ========================================================= |
---|
| 2124 | |
---|
| 2125 | c print *, 'prognostic wake computation' |
---|
| 2126 | |
---|
| 2127 | |
---|
| 2128 | C 2.1 - Undisturbed area and Wake integrals |
---|
| 2129 | C --------------------------------------------------------- |
---|
| 2130 | |
---|
| 2131 | z = 0. |
---|
| 2132 | ktop=0 |
---|
| 2133 | kupper = 0 |
---|
| 2134 | sum_thu = 0. |
---|
| 2135 | sum_tu = 0. |
---|
| 2136 | sum_qu = 0. |
---|
| 2137 | sum_thvu = 0. |
---|
| 2138 | sum_dth = 0. |
---|
| 2139 | sum_dq = 0. |
---|
| 2140 | sum_rho = 0. |
---|
| 2141 | sum_dtdwn = 0. |
---|
| 2142 | sum_dqdwn = 0. |
---|
| 2143 | |
---|
| 2144 | av_thu = 0. |
---|
| 2145 | av_tu =0. |
---|
| 2146 | av_qu =0. |
---|
| 2147 | av_thvu = 0. |
---|
| 2148 | av_dth = 0. |
---|
| 2149 | av_dq = 0. |
---|
| 2150 | av_rho =0. |
---|
| 2151 | av_dtdwn =0. |
---|
| 2152 | av_dqdwn = 0. |
---|
| 2153 | |
---|
| 2154 | C Potential temperatures and humidity |
---|
| 2155 | c---------------------------------------------------------- |
---|
| 2156 | |
---|
| 2157 | DO k =1,klev |
---|
| 2158 | rho(k) = p(k)/(rd*te(k)) |
---|
| 2159 | IF(k .eq. 1) THEN |
---|
| 2160 | rhoh(k) = ph(k)/(rd*te(k)) |
---|
| 2161 | zhh(k)=0 |
---|
| 2162 | ELSE |
---|
| 2163 | rhoh(k) = ph(k)*2./(rd*(te(k)+te(k-1))) |
---|
| 2164 | zhh(k)=(ph(k)-ph(k-1))/(-rhoh(k)*RG)+zhh(k-1) |
---|
| 2165 | ENDIF |
---|
| 2166 | the(k) = te(k)/ppi(k) |
---|
| 2167 | thu(k) = (te(k) - deltatw(k)*sigmaw)/ppi(k) |
---|
| 2168 | tu(k) = te(k) - deltatw(k)*sigmaw |
---|
| 2169 | qu(k) = qe(k) - deltaqw(k)*sigmaw |
---|
| 2170 | rhow(k) = p(k)/(rd*(te(k)+deltatw(k))) |
---|
| 2171 | dth(k) = deltatw(k)/ppi(k) |
---|
| 2172 | LL = (1-sqrt(sigmaw))/sqrt(wdens) |
---|
| 2173 | ENDDO |
---|
| 2174 | |
---|
| 2175 | DO k = 1, klev-1 |
---|
| 2176 | IF(k.eq.1) THEN |
---|
| 2177 | N2(k)=0 |
---|
| 2178 | ELSE |
---|
| 2179 | N2(k)=max(0.,-RG**2/the(k)*rho(k)*(the(k+1)-the(k-1)) |
---|
| 2180 | $ /(p(k+1)-p(k-1))) |
---|
| 2181 | ENDIF |
---|
| 2182 | ZH(k)=(zhh(k)+zhh(k+1))/2 |
---|
| 2183 | |
---|
| 2184 | Cgw(k)=sqrt(N2(k))*ZH(k) |
---|
| 2185 | Tgw(k)=coefgw*Cgw(k)/LL |
---|
| 2186 | ENDDO |
---|
| 2187 | |
---|
| 2188 | N2(klev)=0 |
---|
| 2189 | ZH(klev)=0 |
---|
| 2190 | Cgw(klev)=0 |
---|
| 2191 | Tgw(klev)=0 |
---|
| 2192 | |
---|
| 2193 | c Calcul de la masse volumique moyenne de la colonne |
---|
| 2194 | c----------------------------------------------------------------- |
---|
| 2195 | |
---|
| 2196 | DO k=1,klev |
---|
| 2197 | epaisseur1(k)=0. |
---|
| 2198 | epaisseur2(k)=0. |
---|
| 2199 | ENDDO |
---|
| 2200 | |
---|
| 2201 | epaisseur1(1)= -(Ph(2)-Ph(1))/(rho(1)*rg)+1. |
---|
| 2202 | epaisseur2(1)= -(Ph(2)-Ph(1))/(rho(1)*rg)+1. |
---|
| 2203 | rhow_moyen(1) = rhow(1) |
---|
| 2204 | |
---|
| 2205 | DO k = 2, klev |
---|
| 2206 | epaisseur1(k)= -(Ph(k+1)-Ph(k))/(rho(k)*rg) +1. |
---|
| 2207 | epaisseur2(k)=epaisseur2(k-1)+epaisseur1(k) |
---|
| 2208 | rhow_moyen(k) = (rhow_moyen(k-1)*epaisseur2(k-1)+ |
---|
| 2209 | $ rhow(k)*epaisseur1(k))/epaisseur2(k) |
---|
| 2210 | ENDDO |
---|
| 2211 | |
---|
| 2212 | |
---|
| 2213 | C Choose an integration bound well above wake top |
---|
| 2214 | c----------------------------------------------------------------- |
---|
| 2215 | |
---|
| 2216 | c Pupper = 50000. ! melting level |
---|
| 2217 | Pupper = 60000. |
---|
| 2218 | c Pupper = 70000. |
---|
| 2219 | |
---|
| 2220 | |
---|
| 2221 | C Determine Wake top pressure (Ptop) from buoyancy integral |
---|
| 2222 | C----------------------------------------------------------------- |
---|
| 2223 | |
---|
| 2224 | c-1/ Pressure of the level where dth becomes less than delta_t_min. |
---|
| 2225 | |
---|
| 2226 | Ptop_provis=ph(1) |
---|
| 2227 | DO k= 2,klev |
---|
| 2228 | IF (dth(k) .GT. -delta_t_min .and. |
---|
| 2229 | $ dth(k-1).LT. -delta_t_min) THEN |
---|
| 2230 | Ptop_provis = ((dth(k)+delta_t_min)*p(k-1) |
---|
| 2231 | $ - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1)) |
---|
| 2232 | GO TO 25 |
---|
| 2233 | ENDIF |
---|
| 2234 | ENDDO |
---|
| 2235 | 25 CONTINUE |
---|
| 2236 | |
---|
| 2237 | c-2/ dth integral |
---|
| 2238 | |
---|
| 2239 | sum_dth = 0. |
---|
| 2240 | dthmin = -delta_t_min |
---|
| 2241 | z = 0. |
---|
| 2242 | |
---|
| 2243 | DO k = 1,klev |
---|
| 2244 | dz = -(max(ph(k+1),Ptop_provis)-Ph(k))/(rho(k)*rg) |
---|
| 2245 | IF (dz .le. 0) GO TO 40 |
---|
| 2246 | z = z+dz |
---|
| 2247 | sum_dth = sum_dth + dth(k)*dz |
---|
| 2248 | dthmin = min(dthmin,dth(k)) |
---|
| 2249 | ENDDO |
---|
| 2250 | 40 CONTINUE |
---|
| 2251 | |
---|
| 2252 | c-3/ height of triangle with area= sum_dth and base = dthmin |
---|
| 2253 | |
---|
| 2254 | hw0 = 2.*sum_dth/min(dthmin,-0.5) |
---|
| 2255 | hw0 = max(hwmin,hw0) |
---|
| 2256 | |
---|
| 2257 | c-4/ now, get Ptop |
---|
| 2258 | |
---|
| 2259 | z = 0. |
---|
| 2260 | ptop = ph(1) |
---|
| 2261 | |
---|
| 2262 | DO k = 1,klev |
---|
| 2263 | dz = min(-(ph(k+1)-Ph(k))/(rho(k)*rg),hw0-z) |
---|
| 2264 | IF (dz .le. 0) GO TO 45 |
---|
| 2265 | z = z+dz |
---|
| 2266 | Ptop = Ph(k)-rho(k)*rg*dz |
---|
| 2267 | ENDDO |
---|
| 2268 | 45 CONTINUE |
---|
| 2269 | |
---|
| 2270 | |
---|
| 2271 | C-5/ Determination de ktop et kupper |
---|
| 2272 | |
---|
| 2273 | DO k=klev,1,-1 |
---|
| 2274 | IF (ph(k+1) .lt. ptop) ktop=k |
---|
| 2275 | IF (ph(k+1) .lt. pupper) kupper=k |
---|
| 2276 | ENDDO |
---|
| 2277 | |
---|
| 2278 | c-6/ Correct ktop and ptop |
---|
| 2279 | |
---|
| 2280 | Ptop_new=ptop |
---|
| 2281 | DO k= ktop,2,-1 |
---|
| 2282 | IF (dth(k) .GT. -delta_t_min .and. |
---|
| 2283 | $ dth(k-1).LT. -delta_t_min) THEN |
---|
| 2284 | Ptop_new = ((dth(k)+delta_t_min)*p(k-1) |
---|
| 2285 | $ - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1)) |
---|
| 2286 | GO TO 225 |
---|
| 2287 | ENDIF |
---|
| 2288 | ENDDO |
---|
| 2289 | 225 CONTINUE |
---|
| 2290 | |
---|
| 2291 | ptop = ptop_new |
---|
| 2292 | |
---|
| 2293 | DO k=klev,1,-1 |
---|
| 2294 | IF (ph(k+1) .lt. ptop) ktop=k |
---|
| 2295 | ENDDO |
---|
| 2296 | |
---|
| 2297 | c Set deltatw & deltaqw to 0 above kupper |
---|
| 2298 | c----------------------------------------------------------- |
---|
| 2299 | |
---|
| 2300 | DO k = kupper,klev |
---|
| 2301 | deltatw(k) = 0. |
---|
| 2302 | deltaqw(k) = 0. |
---|
| 2303 | ENDDO |
---|
| 2304 | |
---|
| 2305 | |
---|
| 2306 | C Vertical gradient of LS omega |
---|
| 2307 | C------------------------------------------------------------ |
---|
| 2308 | |
---|
| 2309 | DO k = 1,kupper |
---|
| 2310 | dp_omgb(k) = (omgb(k+1) - omgb(k))/(ph(k+1)-ph(k)) |
---|
| 2311 | ENDDO |
---|
| 2312 | |
---|
| 2313 | |
---|
| 2314 | C Integrals (and wake top level number) |
---|
| 2315 | C ----------------------------------------------------------- |
---|
| 2316 | |
---|
| 2317 | C Initialize sum_thvu to 1st level virt. pot. temp. |
---|
| 2318 | |
---|
| 2319 | z = 1. |
---|
| 2320 | dz = 1. |
---|
| 2321 | sum_thvu = thu(1)*(1.+eps*qu(1))*dz |
---|
| 2322 | sum_dth = 0. |
---|
| 2323 | |
---|
| 2324 | DO k = 1,klev |
---|
| 2325 | dz = -(max(ph(k+1),Ptop)-Ph(k))/(rho(k)*rg) |
---|
| 2326 | IF (dz .LE. 0) GO TO 50 |
---|
| 2327 | z = z+dz |
---|
| 2328 | sum_thu = sum_thu + thu(k)*dz |
---|
| 2329 | sum_tu = sum_tu + tu(k)*dz |
---|
| 2330 | sum_qu = sum_qu + qu(k)*dz |
---|
| 2331 | sum_thvu = sum_thvu + thu(k)*(1.+eps*qu(k))*dz |
---|
| 2332 | sum_dth = sum_dth + dth(k)*dz |
---|
| 2333 | sum_dq = sum_dq + deltaqw(k)*dz |
---|
| 2334 | sum_rho = sum_rho + rhow(k)*dz |
---|
| 2335 | sum_dtdwn = sum_dtdwn + dtdwn(k)*dz |
---|
| 2336 | sum_dqdwn = sum_dqdwn + dqdwn(k)*dz |
---|
| 2337 | ENDDO |
---|
| 2338 | 50 CONTINUE |
---|
| 2339 | |
---|
| 2340 | hw0 = z |
---|
| 2341 | |
---|
| 2342 | C 2.1 - WAPE and mean forcing computation |
---|
| 2343 | C------------------------------------------------------------- |
---|
| 2344 | |
---|
| 2345 | C Means |
---|
| 2346 | |
---|
| 2347 | av_thu = sum_thu/hw0 |
---|
| 2348 | av_tu = sum_tu/hw0 |
---|
| 2349 | av_qu = sum_qu/hw0 |
---|
| 2350 | av_thvu = sum_thvu/hw0 |
---|
| 2351 | c av_thve = sum_thve/hw0 |
---|
| 2352 | av_dth = sum_dth/hw0 |
---|
| 2353 | av_dq = sum_dq/hw0 |
---|
| 2354 | av_rho = sum_rho/hw0 |
---|
| 2355 | av_dtdwn = sum_dtdwn/hw0 |
---|
| 2356 | av_dqdwn = sum_dqdwn/hw0 |
---|
| 2357 | |
---|
| 2358 | wape = - rg*hw0*(av_dth |
---|
| 2359 | $ + eps*(av_thu*av_dq+av_dth*av_qu+av_dth*av_dq ))/av_thvu |
---|
| 2360 | |
---|
| 2361 | C 2.2 Prognostic variable update |
---|
| 2362 | C ------------------------------------------------------------ |
---|
| 2363 | |
---|
| 2364 | C Filter out bad wakes |
---|
| 2365 | |
---|
| 2366 | IF ( wape .LT. 0.) THEN |
---|
[953] | 2367 | if(prt_level.ge.10) print*,'wape<0' |
---|
[879] | 2368 | wape = 0. |
---|
| 2369 | hw = hwmin |
---|
| 2370 | sigmaw = max(sigmad,sigd_con) |
---|
| 2371 | fip = 0. |
---|
| 2372 | DO k = 1,klev |
---|
| 2373 | deltatw(k) = 0. |
---|
| 2374 | deltaqw(k) = 0. |
---|
| 2375 | dth(k) = 0. |
---|
| 2376 | ENDDO |
---|
| 2377 | ELSE |
---|
[953] | 2378 | if(prt_level.ge.10) print*,'wape>0' |
---|
[879] | 2379 | Cstar = stark*sqrt(2.*wape) |
---|
| 2380 | ENDIF |
---|
| 2381 | |
---|
| 2382 | C------------------------------------------------------------------ |
---|
| 2383 | C Sub-time-stepping |
---|
| 2384 | C------------------------------------------------------------------ |
---|
| 2385 | |
---|
| 2386 | c nsub=36 |
---|
| 2387 | nsub=10 |
---|
| 2388 | dtimesub=dtime/nsub |
---|
| 2389 | |
---|
| 2390 | c------------------------------------------------------------ |
---|
| 2391 | DO isubstep = 1,nsub |
---|
| 2392 | c------------------------------------------------------------ |
---|
| 2393 | |
---|
| 2394 | c print*,'---------------','substep=',isubstep,'-------------' |
---|
| 2395 | |
---|
| 2396 | c Evolution of sigmaw |
---|
| 2397 | |
---|
| 2398 | |
---|
| 2399 | gfl = 2.*sqrt(3.14*wdens*sigmaw) |
---|
| 2400 | |
---|
| 2401 | sigmaw =sigmaw + gfl*Cstar*dtimesub |
---|
| 2402 | sigmaw =min(sigmaw,0.99) !!!!!!!! |
---|
| 2403 | c wdens = wdens0/(10.*sigmaw) |
---|
| 2404 | c sigmaw =max(sigmaw,sigd_con) |
---|
| 2405 | c sigmaw =max(sigmaw,sigmad) |
---|
| 2406 | |
---|
| 2407 | c calcul de la difference de vitesse verticale poche - zone non perturbee |
---|
| 2408 | |
---|
| 2409 | z= 0. |
---|
| 2410 | dp_deltomg(1:klev)=0. |
---|
| 2411 | omg(1:klev+1)=0. |
---|
| 2412 | |
---|
| 2413 | omg(1) = 0. |
---|
| 2414 | dp_deltomg(1) = -(gfl*Cstar)/(sigmaw * (1-sigmaw)) |
---|
| 2415 | |
---|
| 2416 | DO k=2,ktop |
---|
| 2417 | dz = -(Ph(k)-Ph(k-1))/(rho(k-1)*rg) |
---|
| 2418 | z = z+dz |
---|
| 2419 | dp_deltomg(k)= dp_deltomg(1) |
---|
| 2420 | omg(k)= dp_deltomg(1)*z |
---|
| 2421 | ENDDO |
---|
| 2422 | |
---|
| 2423 | dztop=-(Ptop-Ph(ktop))/(rho(ktop)*rg) |
---|
| 2424 | ztop = z+dztop |
---|
| 2425 | omgtop=dp_deltomg(1)*ztop |
---|
| 2426 | |
---|
| 2427 | |
---|
| 2428 | c Conversion de la vitesse verticale de m/s a Pa/s |
---|
| 2429 | |
---|
| 2430 | omgtop = -rho(ktop)*rg*omgtop |
---|
| 2431 | dp_deltomg(1) = omgtop/(ptop-ph(1)) |
---|
| 2432 | |
---|
| 2433 | DO k = 1,ktop |
---|
| 2434 | omg(k) = - rho(k)*rg*omg(k) |
---|
| 2435 | dp_deltomg(k) = dp_deltomg(1) |
---|
| 2436 | ENDDO |
---|
| 2437 | |
---|
| 2438 | c raccordement lineaire de omg de ptop a pupper |
---|
| 2439 | |
---|
| 2440 | IF (kupper .GT. ktop) THEN |
---|
| 2441 | omg(kupper+1) = - Rg*amdwn(kupper+1)/sigmaw |
---|
| 2442 | $ + Rg*amup(kupper+1)/(1.-sigmaw) |
---|
| 2443 | dp_deltomg(kupper) = (omgtop-omg(kupper+1))/(Ptop-Pupper) |
---|
| 2444 | DO k=ktop+1,kupper |
---|
| 2445 | dp_deltomg(k) = dp_deltomg(kupper) |
---|
| 2446 | omg(k) = omgtop+(ph(k)-Ptop)*dp_deltomg(kupper) |
---|
| 2447 | ENDDO |
---|
| 2448 | ENDIF |
---|
| 2449 | |
---|
| 2450 | c Compute wake average vertical velocity omgbw |
---|
| 2451 | |
---|
| 2452 | DO k = 1,klev+1 |
---|
| 2453 | omgbw(k) = omgb(k)+(1.-sigmaw)*omg(k) |
---|
| 2454 | ENDDO |
---|
| 2455 | |
---|
| 2456 | c and its vertical gradient dp_omgbw |
---|
| 2457 | |
---|
| 2458 | DO k = 1,klev |
---|
| 2459 | dp_omgbw(k) = (omgbw(k+1)-omgbw(k))/(ph(k+1)-ph(k)) |
---|
| 2460 | ENDDO |
---|
| 2461 | |
---|
| 2462 | |
---|
| 2463 | c Upstream coefficients for omgb velocity |
---|
| 2464 | c-- (alpha_up(k) is the coefficient of the value at level k) |
---|
| 2465 | c-- (1-alpha_up(k) is the coefficient of the value at level k-1) |
---|
| 2466 | |
---|
| 2467 | DO k = 1,klev |
---|
| 2468 | alpha_up(k) = 0. |
---|
| 2469 | IF (omgb(k) .GT. 0.) alpha_up(k) = 1. |
---|
| 2470 | ENDDO |
---|
| 2471 | |
---|
| 2472 | c Matrix expressing [The,deltatw] from [Th1,Th2] |
---|
| 2473 | |
---|
| 2474 | RRe1 = 1.-sigmaw |
---|
| 2475 | RRe2 = sigmaw |
---|
| 2476 | RRd1 = -1. |
---|
| 2477 | RRd2 = 1. |
---|
| 2478 | |
---|
| 2479 | c Get [Th1,Th2], dth and [q1,q2] |
---|
| 2480 | |
---|
| 2481 | DO k = 1,kupper+1 |
---|
| 2482 | dth(k) = deltatw(k)/ppi(k) |
---|
| 2483 | Th1(k) = the(k) - sigmaw *dth(k) ! undisturbed area |
---|
| 2484 | Th2(k) = the(k) + (1.-sigmaw)*dth(k) ! wake |
---|
| 2485 | q1(k) = qe(k) - sigmaw *deltaqw(k) ! undisturbed area |
---|
| 2486 | q2(k) = qe(k) + (1.-sigmaw)*deltaqw(k) ! wake |
---|
| 2487 | ENDDO |
---|
| 2488 | |
---|
| 2489 | D_Th1(1) = 0. |
---|
| 2490 | D_Th2(1) = 0. |
---|
| 2491 | D_dth(1) = 0. |
---|
| 2492 | D_q1(1) = 0. |
---|
| 2493 | D_q2(1) = 0. |
---|
| 2494 | D_dq(1) = 0. |
---|
| 2495 | |
---|
| 2496 | DO k = 2,kupper+1 ! loop on interfaces |
---|
| 2497 | D_Th1(k) = Th1(k-1)-Th1(k) |
---|
| 2498 | D_Th2(k) = Th2(k-1)-Th2(k) |
---|
| 2499 | D_dth(k) = dth(k-1)-dth(k) |
---|
| 2500 | D_q1(k) = q1(k-1)-q1(k) |
---|
| 2501 | D_q2(k) = q2(k-1)-q2(k) |
---|
| 2502 | D_dq(k) = deltaqw(k-1)-deltaqw(k) |
---|
| 2503 | ENDDO |
---|
| 2504 | |
---|
| 2505 | omgbdth(1) = 0. |
---|
| 2506 | omgbdq(1) = 0. |
---|
| 2507 | |
---|
| 2508 | DO k = 2,kupper+1 ! loop on interfaces |
---|
| 2509 | omgbdth(k) = omgb(k)*( dth(k-1) - dth(k)) |
---|
| 2510 | omgbdq(k) = omgb(k)*(deltaqw(k-1) - deltaqw(k)) |
---|
| 2511 | ENDDO |
---|
| 2512 | |
---|
| 2513 | |
---|
| 2514 | c----------------------------------------------------------------- |
---|
| 2515 | DO k=1,kupper-1 |
---|
| 2516 | c----------------------------------------------------------------- |
---|
| 2517 | c |
---|
| 2518 | c Compute redistribution (advective) term |
---|
| 2519 | c |
---|
| 2520 | d_deltatw(k) = |
---|
| 2521 | $ dtimesub/(Ph(k)-Ph(k+1))*( |
---|
| 2522 | $ RRd1*omg(k )*sigmaw *D_Th1(k) |
---|
| 2523 | $ -RRd2*omg(k+1)*(1.-sigmaw)*D_Th2(k+1) |
---|
| 2524 | $ -(1.-alpha_up(k))*omgbdth(k) - alpha_up(k+1)*omgbdth(k+1) |
---|
| 2525 | $ )*ppi(k) |
---|
| 2526 | c print*,'d_deltatw=',d_deltatw(k) |
---|
| 2527 | c |
---|
| 2528 | d_deltaqw(k) = |
---|
| 2529 | $ dtimesub/(Ph(k)-Ph(k+1))*( |
---|
| 2530 | $ RRd1*omg(k )*sigmaw *D_q1(k) |
---|
| 2531 | $ -RRd2*omg(k+1)*(1.-sigmaw)*D_q2(k+1) |
---|
| 2532 | $ -(1.-alpha_up(k))*omgbdq(k) - alpha_up(k+1)*omgbdq(k+1) |
---|
| 2533 | $ ) |
---|
| 2534 | c print*,'d_deltaqw=',d_deltaqw(k) |
---|
| 2535 | c |
---|
| 2536 | c and increment large scale tendencies |
---|
| 2537 | c |
---|
| 2538 | dtls(k) = dtls(k) + |
---|
| 2539 | $ dtimesub*( |
---|
| 2540 | $ ( RRe1*omg(k )*sigmaw *D_Th1(k) |
---|
| 2541 | $ -RRe2*omg(k+1)*(1.-sigmaw)*D_Th2(k+1) ) |
---|
| 2542 | $ /(Ph(k)-Ph(k+1)) |
---|
| 2543 | $ -sigmaw*(1.-sigmaw)*dth(k)*dp_deltomg(k) |
---|
| 2544 | $ )*ppi(k) |
---|
| 2545 | c print*,'dtls=',dtls(k) |
---|
| 2546 | c |
---|
| 2547 | dqls(k) = dqls(k) + |
---|
| 2548 | $ dtimesub*( |
---|
| 2549 | $ ( RRe1*omg(k )*sigmaw *D_q1(k) |
---|
| 2550 | $ -RRe2*omg(k+1)*(1.-sigmaw)*D_q2(k+1) ) |
---|
| 2551 | $ /(Ph(k)-Ph(k+1)) |
---|
| 2552 | $ -sigmaw*(1.-sigmaw)*deltaqw(k)*dp_deltomg(k) |
---|
| 2553 | $ ) |
---|
| 2554 | c print*,'dqls=',dqls(k) |
---|
| 2555 | |
---|
| 2556 | c------------------------------------------------------------------- |
---|
| 2557 | ENDDO |
---|
| 2558 | c------------------------------------------------------------------ |
---|
| 2559 | |
---|
| 2560 | C Increment state variables |
---|
| 2561 | |
---|
| 2562 | DO k = 1,kupper-1 |
---|
| 2563 | |
---|
| 2564 | c Coefficient de répartition |
---|
| 2565 | |
---|
| 2566 | Crep(k)=Crep_sol*(ph(kupper)-ph(k))/(ph(kupper)-ph(1)) |
---|
| 2567 | Crep(k)=Crep(k)+Crep_upper*(ph(1)-ph(k))/(p(1)-ph(kupper)) |
---|
| 2568 | |
---|
| 2569 | |
---|
| 2570 | c Reintroduce compensating subsidence term. |
---|
| 2571 | |
---|
| 2572 | c dtKE(k)=(dtdwn(k)*Crep(k))/sigmaw |
---|
| 2573 | c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)) |
---|
| 2574 | c . /(1-sigmaw) |
---|
| 2575 | c dqKE(k)=(dqdwn(k)*Crep(k))/sigmaw |
---|
| 2576 | c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)) |
---|
| 2577 | c . /(1-sigmaw) |
---|
| 2578 | c |
---|
| 2579 | c dtKE(k)=(dtdwn(k)*Crep(k)+(1-Crep(k))*dta(k))/sigmaw |
---|
| 2580 | c dtKE(k)=dtKE(k)-(dtdwn(k)*(1-Crep(k))+dta(k)*Crep(k)) |
---|
| 2581 | c . /(1-sigmaw) |
---|
| 2582 | c dqKE(k)=(dqdwn(k)*Crep(k)+(1-Crep(k))*dqa(k))/sigmaw |
---|
| 2583 | c dqKE(k)=dqKE(k)-(dqdwn(k)*(1-Crep(k))+dqa(k)*Crep(k)) |
---|
| 2584 | c . /(1-sigmaw) |
---|
| 2585 | |
---|
| 2586 | dtKE(k)=(dtdwn(k)/sigmaw - dta(k)/(1.-sigmaw)) |
---|
| 2587 | dqKE(k)=(dqdwn(k)/sigmaw - dqa(k)/(1.-sigmaw)) |
---|
| 2588 | c print*,'dtKE=',dtKE(k) |
---|
| 2589 | c print*,'dqKE=',dqKE(k) |
---|
| 2590 | c |
---|
| 2591 | dtPBL(k)=(wdtPBL(k)/sigmaw - udtPBL(k)/(1.-sigmaw)) |
---|
| 2592 | dqPBL(k)=(wdqPBL(k)/sigmaw - udqPBL(k)/(1.-sigmaw)) |
---|
| 2593 | c |
---|
| 2594 | spread(k) = (1.-sigmaw)*dp_deltomg(k)+gfl*Cstar/sigmaw |
---|
| 2595 | c print*,'spread=',spread(k) |
---|
| 2596 | |
---|
| 2597 | |
---|
| 2598 | c ajout d'un effet onde de gravité -Tgw(k)*deltatw(k) 03/02/06 YU Jingmei |
---|
| 2599 | |
---|
| 2600 | d_deltat_gw(k)=d_deltat_gw(k)-Tgw(k)*deltatw(k)* dtimesub |
---|
| 2601 | c print*,'d_delta_gw=',d_deltat_gw(k) |
---|
| 2602 | ff=d_deltatw(k)/dtimesub |
---|
| 2603 | |
---|
| 2604 | c Sans GW |
---|
| 2605 | c |
---|
| 2606 | c deltatw(k)=deltatw(k)+dtimesub*(ff+dtKE(k)-spread(k)*deltatw(k)) |
---|
| 2607 | c |
---|
| 2608 | c GW formule 1 |
---|
| 2609 | c |
---|
| 2610 | c deltatw(k) = deltatw(k)+dtimesub* |
---|
| 2611 | c $ (ff+dtKE(k) - spread(k)*deltatw(k)-Tgw(k)*deltatw(k)) |
---|
| 2612 | c |
---|
| 2613 | c GW formule 2 |
---|
| 2614 | |
---|
| 2615 | IF (dtimesub*Tgw(k).lt.1.e-10) THEN |
---|
| 2616 | deltatw(k) = deltatw(k)+dtimesub* |
---|
| 2617 | $ (ff+dtKE(k)+dtPBL(k) |
---|
| 2618 | $ - spread(k)*deltatw(k)-Tgw(k)*deltatw(k)) |
---|
| 2619 | ELSE |
---|
| 2620 | deltatw(k) = deltatw(k)+1/Tgw(k)*(1-exp(-dtimesub*Tgw(k)))* |
---|
| 2621 | $ (ff+dtKE(k)+dtPBL(k) |
---|
| 2622 | $ - spread(k)*deltatw(k)-Tgw(k)*deltatw(k)) |
---|
| 2623 | ENDIF |
---|
| 2624 | |
---|
| 2625 | dth(k) = deltatw(k)/ppi(k) |
---|
| 2626 | |
---|
| 2627 | gg=d_deltaqw(k)/dtimesub |
---|
| 2628 | |
---|
| 2629 | deltaqw(k) = deltaqw(k) + |
---|
| 2630 | $ dtimesub*(gg+ dqKE(k)+dqPBL(k) - spread(k)*deltaqw(k)) |
---|
| 2631 | |
---|
| 2632 | d_deltatw2(k)=d_deltatw2(k)+d_deltatw(k) |
---|
| 2633 | d_deltaqw2(k)=d_deltaqw2(k)+d_deltaqw(k) |
---|
| 2634 | ENDDO |
---|
| 2635 | |
---|
| 2636 | C And update large scale variables |
---|
| 2637 | |
---|
| 2638 | DO k = 1,kupper |
---|
| 2639 | te(k) = te0(k) + dtls(k) |
---|
| 2640 | qe(k) = qe0(k) + dqls(k) |
---|
| 2641 | the(k) = te(k)/ppi(k) |
---|
| 2642 | ENDDO |
---|
| 2643 | |
---|
| 2644 | c Determine Ptop from buoyancy integral |
---|
| 2645 | c---------------------------------------------------------------------- |
---|
| 2646 | |
---|
| 2647 | c-1/ Pressure of the level where dth changes sign. |
---|
| 2648 | |
---|
| 2649 | Ptop_provis=ph(1) |
---|
| 2650 | |
---|
| 2651 | DO k= 2,klev |
---|
| 2652 | IF (dth(k) .GT. -delta_t_min .and. |
---|
| 2653 | $ dth(k-1).LT. -delta_t_min) THEN |
---|
| 2654 | Ptop_provis = ((dth(k)+delta_t_min)*p(k-1) |
---|
| 2655 | $ - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1)) |
---|
| 2656 | GO TO 65 |
---|
| 2657 | ENDIF |
---|
| 2658 | ENDDO |
---|
| 2659 | 65 CONTINUE |
---|
| 2660 | |
---|
| 2661 | c-2/ dth integral |
---|
| 2662 | |
---|
| 2663 | sum_dth = 0. |
---|
| 2664 | dthmin = -delta_t_min |
---|
| 2665 | z = 0. |
---|
| 2666 | |
---|
| 2667 | DO k = 1,klev |
---|
| 2668 | dz = -(max(ph(k+1),Ptop_provis)-Ph(k))/(rho(k)*rg) |
---|
| 2669 | IF (dz .le. 0) GO TO 70 |
---|
| 2670 | z = z+dz |
---|
| 2671 | sum_dth = sum_dth + dth(k)*dz |
---|
| 2672 | dthmin = min(dthmin,dth(k)) |
---|
| 2673 | ENDDO |
---|
| 2674 | 70 CONTINUE |
---|
| 2675 | |
---|
| 2676 | c-3/ height of triangle with area= sum_dth and base = dthmin |
---|
| 2677 | |
---|
| 2678 | hw = 2.*sum_dth/min(dthmin,-0.5) |
---|
| 2679 | hw = max(hwmin,hw) |
---|
| 2680 | |
---|
| 2681 | c-4/ now, get Ptop |
---|
| 2682 | |
---|
| 2683 | ktop = 0 |
---|
| 2684 | z=0. |
---|
| 2685 | |
---|
| 2686 | DO k = 1,klev |
---|
| 2687 | dz = min(-(ph(k+1)-Ph(k))/(rho(k)*rg),hw-z) |
---|
| 2688 | IF (dz .le. 0) GO TO 75 |
---|
| 2689 | z = z+dz |
---|
| 2690 | Ptop = Ph(k)-rho(k)*rg*dz |
---|
| 2691 | ktop = k |
---|
| 2692 | ENDDO |
---|
| 2693 | 75 CONTINUE |
---|
| 2694 | |
---|
| 2695 | c-5/Correct ktop and ptop |
---|
| 2696 | |
---|
| 2697 | Ptop_new=ptop |
---|
| 2698 | |
---|
| 2699 | DO k= ktop,2,-1 |
---|
| 2700 | IF (dth(k) .GT. -delta_t_min .and. |
---|
| 2701 | $ dth(k-1).LT. -delta_t_min) THEN |
---|
| 2702 | Ptop_new = ((dth(k)+delta_t_min)*p(k-1) |
---|
| 2703 | $ - (dth(k-1)+delta_t_min)*p(k)) /(dth(k) - dth(k-1)) |
---|
| 2704 | GO TO 275 |
---|
| 2705 | ENDIF |
---|
| 2706 | ENDDO |
---|
| 2707 | 275 CONTINUE |
---|
| 2708 | |
---|
| 2709 | ptop = ptop_new |
---|
| 2710 | |
---|
| 2711 | DO k=klev,1,-1 |
---|
| 2712 | IF (ph(k+1) .LT. ptop) ktop=k |
---|
| 2713 | ENDDO |
---|
| 2714 | |
---|
| 2715 | c-6/ Set deltatw & deltaqw to 0 above kupper |
---|
| 2716 | |
---|
| 2717 | DO k = kupper,klev |
---|
| 2718 | deltatw(k) = 0. |
---|
| 2719 | deltaqw(k) = 0. |
---|
| 2720 | ENDDO |
---|
| 2721 | |
---|
| 2722 | c------------------------------------------------------------------ |
---|
| 2723 | ENDDO ! end sub-timestep loop |
---|
| 2724 | C ----------------------------------------------------------------- |
---|
| 2725 | |
---|
| 2726 | c Get back to tendencies per second |
---|
| 2727 | |
---|
| 2728 | DO k = 1,kupper-1 |
---|
| 2729 | dtls(k) = dtls(k)/dtime |
---|
| 2730 | dqls(k) = dqls(k)/dtime |
---|
| 2731 | d_deltatw2(k)=d_deltatw2(k)/dtime |
---|
| 2732 | d_deltaqw2(k)=d_deltaqw2(k)/dtime |
---|
| 2733 | d_deltat_gw(k) = d_deltat_gw(k)/dtime |
---|
| 2734 | ENDDO |
---|
| 2735 | |
---|
| 2736 | C 2.1 - Undisturbed area and Wake integrals |
---|
| 2737 | C --------------------------------------------------------- |
---|
| 2738 | |
---|
| 2739 | z = 0. |
---|
| 2740 | sum_thu = 0. |
---|
| 2741 | sum_tu = 0. |
---|
| 2742 | sum_qu = 0. |
---|
| 2743 | sum_thvu = 0. |
---|
| 2744 | sum_dth = 0. |
---|
| 2745 | sum_dq = 0. |
---|
| 2746 | sum_rho = 0. |
---|
| 2747 | sum_dtdwn = 0. |
---|
| 2748 | sum_dqdwn = 0. |
---|
| 2749 | |
---|
| 2750 | av_thu = 0. |
---|
| 2751 | av_tu =0. |
---|
| 2752 | av_qu =0. |
---|
| 2753 | av_thvu = 0. |
---|
| 2754 | av_dth = 0. |
---|
| 2755 | av_dq = 0. |
---|
| 2756 | av_rho =0. |
---|
| 2757 | av_dtdwn =0. |
---|
| 2758 | av_dqdwn = 0. |
---|
| 2759 | |
---|
| 2760 | C Potential temperatures and humidity |
---|
| 2761 | c---------------------------------------------------------- |
---|
| 2762 | |
---|
| 2763 | DO k =1,klev |
---|
| 2764 | rho(k) = p(k)/(rd*te(k)) |
---|
| 2765 | IF(k .eq. 1) THEN |
---|
| 2766 | rhoh(k) = ph(k)/(rd*te(k)) |
---|
| 2767 | zhh(k)=0 |
---|
| 2768 | ELSE |
---|
| 2769 | rhoh(k) = ph(k)*2./(rd*(te(k)+te(k-1))) |
---|
| 2770 | zhh(k)=(ph(k)-ph(k-1))/(-rhoh(k)*RG)+zhh(k-1) |
---|
| 2771 | ENDIF |
---|
| 2772 | the(k) = te(k)/ppi(k) |
---|
| 2773 | thu(k) = (te(k) - deltatw(k)*sigmaw)/ppi(k) |
---|
| 2774 | tu(k) = te(k) - deltatw(k)*sigmaw |
---|
| 2775 | qu(k) = qe(k) - deltaqw(k)*sigmaw |
---|
| 2776 | rhow(k) = p(k)/(rd*(te(k)+deltatw(k))) |
---|
| 2777 | dth(k) = deltatw(k)/ppi(k) |
---|
| 2778 | |
---|
| 2779 | ENDDO |
---|
| 2780 | |
---|
| 2781 | C Integrals (and wake top level number) |
---|
| 2782 | C ----------------------------------------------------------- |
---|
| 2783 | |
---|
| 2784 | C Initialize sum_thvu to 1st level virt. pot. temp. |
---|
| 2785 | |
---|
| 2786 | z = 1. |
---|
| 2787 | dz = 1. |
---|
| 2788 | sum_thvu = thu(1)*(1.+eps*qu(1))*dz |
---|
| 2789 | sum_dth = 0. |
---|
| 2790 | |
---|
| 2791 | DO k = 1,klev |
---|
| 2792 | dz = -(max(ph(k+1),Ptop)-Ph(k))/(rho(k)*rg) |
---|
| 2793 | |
---|
| 2794 | IF (dz .LE. 0) GO TO 51 |
---|
| 2795 | z = z+dz |
---|
| 2796 | sum_thu = sum_thu + thu(k)*dz |
---|
| 2797 | sum_tu = sum_tu + tu(k)*dz |
---|
| 2798 | sum_qu = sum_qu + qu(k)*dz |
---|
| 2799 | sum_thvu = sum_thvu + thu(k)*(1.+eps*qu(k))*dz |
---|
| 2800 | sum_dth = sum_dth + dth(k)*dz |
---|
| 2801 | sum_dq = sum_dq + deltaqw(k)*dz |
---|
| 2802 | sum_rho = sum_rho + rhow(k)*dz |
---|
| 2803 | sum_dtdwn = sum_dtdwn + dtdwn(k)*dz |
---|
| 2804 | sum_dqdwn = sum_dqdwn + dqdwn(k)*dz |
---|
| 2805 | ENDDO |
---|
| 2806 | 51 CONTINUE |
---|
| 2807 | |
---|
| 2808 | hw0 = z |
---|
| 2809 | |
---|
| 2810 | C 2.1 - WAPE and mean forcing computation |
---|
| 2811 | C------------------------------------------------------------- |
---|
| 2812 | |
---|
| 2813 | C Means |
---|
| 2814 | |
---|
| 2815 | av_thu = sum_thu/hw0 |
---|
| 2816 | av_tu = sum_tu/hw0 |
---|
| 2817 | av_qu = sum_qu/hw0 |
---|
| 2818 | av_thvu = sum_thvu/hw0 |
---|
| 2819 | av_dth = sum_dth/hw0 |
---|
| 2820 | av_dq = sum_dq/hw0 |
---|
| 2821 | av_rho = sum_rho/hw0 |
---|
| 2822 | av_dtdwn = sum_dtdwn/hw0 |
---|
| 2823 | av_dqdwn = sum_dqdwn/hw0 |
---|
| 2824 | |
---|
| 2825 | wape2 = - rg*hw0*(av_dth |
---|
| 2826 | $ + eps*(av_thu*av_dq+av_dth*av_qu+av_dth*av_dq ))/av_thvu |
---|
| 2827 | |
---|
| 2828 | |
---|
| 2829 | C 2.2 Prognostic variable update |
---|
| 2830 | C ------------------------------------------------------------ |
---|
| 2831 | |
---|
| 2832 | C Filter out bad wakes |
---|
| 2833 | |
---|
| 2834 | IF ( wape2 .LT. 0.) THEN |
---|
[953] | 2835 | if(prt_level.ge.10) print*,'wape2<0' |
---|
[879] | 2836 | wape2 = 0. |
---|
| 2837 | hw = hwmin |
---|
| 2838 | sigmaw = max(sigmad,sigd_con) |
---|
| 2839 | fip = 0. |
---|
| 2840 | DO k = 1,klev |
---|
| 2841 | deltatw(k) = 0. |
---|
| 2842 | deltaqw(k) = 0. |
---|
| 2843 | dth(k) = 0. |
---|
| 2844 | ENDDO |
---|
| 2845 | ELSE |
---|
[953] | 2846 | if(prt_level.ge.10) print*,'wape2>0' |
---|
[879] | 2847 | Cstar2 = stark*sqrt(2.*wape2) |
---|
| 2848 | |
---|
| 2849 | ENDIF |
---|
| 2850 | |
---|
| 2851 | ktopw = ktop |
---|
| 2852 | |
---|
| 2853 | IF (ktopw .gt. 0) then |
---|
| 2854 | |
---|
| 2855 | Cjyg1 Utilisation d'un h_efficace constant ( ~ feeding layer) |
---|
| 2856 | ccc heff = 600. |
---|
| 2857 | C Utilisation de la hauteur hw |
---|
| 2858 | cc heff = 0.7*hw |
---|
| 2859 | heff = hw |
---|
| 2860 | |
---|
| 2861 | FIP = 0.5*rho(ktopw)*Cstar2**3*heff*2*sqrt(sigmaw*wdens*3.14) |
---|
| 2862 | FIP = alpk * FIP |
---|
| 2863 | Cjyg2 |
---|
| 2864 | ELSE |
---|
| 2865 | FIP = 0. |
---|
| 2866 | ENDIF |
---|
| 2867 | |
---|
| 2868 | |
---|
| 2869 | C Limitation de sigmaw |
---|
| 2870 | c |
---|
| 2871 | C sécurité : si le wake occuppe plus de 90 % de la surface de la maille, |
---|
| 2872 | C alors il disparait en se mélangeant à la partie undisturbed |
---|
| 2873 | |
---|
[1146] | 2874 | ! correction NICOLAS . ((wape.ge.wape2).and.(wape2.le.1.0))) THEN |
---|
[879] | 2875 | IF ((sigmaw.GT.0.9).or. |
---|
[1059] | 2876 | . ((wape.ge.wape2).and.(wape2.le.1.0)).or.(ktopw.le.2)) THEN |
---|
| 2877 | cIM cf NR/JYG 251108 . ((wape.ge.wape2).and.(wape2.le.1.0))) THEN |
---|
[879] | 2878 | c IF (sigmaw.GT.0.9) THEN |
---|
| 2879 | DO k = 1,klev |
---|
| 2880 | dtls(k) = 0. |
---|
| 2881 | dqls(k) = 0. |
---|
| 2882 | deltatw(k) = 0. |
---|
| 2883 | deltaqw(k) = 0. |
---|
| 2884 | ENDDO |
---|
| 2885 | wape = 0. |
---|
| 2886 | hw = hwmin |
---|
| 2887 | sigmaw = sigmad |
---|
| 2888 | fip = 0. |
---|
| 2889 | ENDIF |
---|
| 2890 | |
---|
| 2891 | RETURN |
---|
| 2892 | END |
---|
| 2893 | |
---|
| 2894 | |
---|
| 2895 | |
---|