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