| 1 | MODULE rocketduststorm_mod |
|---|
| 2 | |
|---|
| 3 | IMPLICIT NONE |
|---|
| 4 | |
|---|
| 5 | REAL, SAVE, ALLOCATABLE :: dustliftday(:) ! dust lifting rate (s-1) |
|---|
| 6 | REAL, SAVE, ALLOCATABLE :: alpha_hmons(:) ! slope winds lifting mesh fraction |
|---|
| 7 | |
|---|
| 8 | CONTAINS |
|---|
| 9 | |
|---|
| 10 | !======================================================================= |
|---|
| 11 | ! ROCKET DUST STORM - vertical transport and detrainment |
|---|
| 12 | !======================================================================= |
|---|
| 13 | ! calculating the vertical flux |
|---|
| 14 | ! calling vl_storm : transport scheme of stormdust tracers |
|---|
| 15 | ! detrainement of stormdust into nomal background dust |
|---|
| 16 | ! ----------------------------------------------------------------------- |
|---|
| 17 | ! Authors: C. Wang; F. Forget; T. Bertrand |
|---|
| 18 | ! Institution: Laboratoire de Meteorologie Dynamique (LMD) Paris, France |
|---|
| 19 | ! ----------------------------------------------------------------------- |
|---|
| 20 | |
|---|
| 21 | SUBROUTINE rocketduststorm(ngrid,nlayer,nq,ptime,ptimestep, & |
|---|
| 22 | pq,pdqfi,pt,pdtfi,pplev,pplay,pzlev, & |
|---|
| 23 | pzlay,pdtsw,pdtlw, & |
|---|
| 24 | ! input for radiative transfer |
|---|
| 25 | clearatm,icount,zday,zls, & |
|---|
| 26 | tsurf,igout,totstormfract, & |
|---|
| 27 | ! input sub-grid scale cloud |
|---|
| 28 | clearsky,totcloudfrac, & |
|---|
| 29 | ! output |
|---|
| 30 | pdqrds,wspeed,dsodust,dsords) |
|---|
| 31 | |
|---|
| 32 | use tracer_mod, only: igcm_stormdust_mass,igcm_stormdust_number & |
|---|
| 33 | ,igcm_dust_mass,igcm_dust_number & |
|---|
| 34 | ,rho_dust |
|---|
| 35 | USE comcstfi_h, only: r,g,cpp,rcp |
|---|
| 36 | use dimradmars_mod, only: albedo,naerkind |
|---|
| 37 | use comsaison_h, only: dist_sol,mu0,fract |
|---|
| 38 | use surfdat_h, only: emis,co2ice,zmea, zstd, zsig, hmons |
|---|
| 39 | use planetwide_mod, only: planetwide_maxval,planetwide_minval |
|---|
| 40 | ! use rocketstorm_h, only: rdsinject |
|---|
| 41 | use callradite_mod |
|---|
| 42 | implicit none |
|---|
| 43 | |
|---|
| 44 | !-------------------------------------------------------- |
|---|
| 45 | ! Input Variables |
|---|
| 46 | !-------------------------------------------------------- |
|---|
| 47 | |
|---|
| 48 | INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points |
|---|
| 49 | INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points |
|---|
| 50 | INTEGER, INTENT(IN) :: nq ! number of tracer species |
|---|
| 51 | REAL, INTENT(IN) :: ptime |
|---|
| 52 | REAL, INTENT(IN) :: ptimestep |
|---|
| 53 | |
|---|
| 54 | REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) ! advected field nq |
|---|
| 55 | REAL, INTENT(IN) :: pdqfi(ngrid,nlayer,nq)! tendancy field pq |
|---|
| 56 | REAL, INTENT(IN) :: pt(ngrid,nlayer) ! temperature at mid-layer (K) |
|---|
| 57 | REAL, INTENT(IN) :: pdtfi(ngrid,nlayer) ! tendancy temperature at mid-layer (K) |
|---|
| 58 | |
|---|
| 59 | REAL, INTENT(IN) :: pplay(ngrid,nlayer) ! pressure at middle of the layers |
|---|
| 60 | REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) ! pressure at intermediate levels |
|---|
| 61 | REAL, INTENT(IN) :: pzlay(ngrid,nlayer) ! altitude at the middle of the layers |
|---|
| 62 | REAL, INTENT(IN) :: pzlev(ngrid,nlayer+1) ! altitude at layer boundaries |
|---|
| 63 | |
|---|
| 64 | REAL, INTENT(IN) :: pdtsw(ngrid,nlayer) ! (K/s) env |
|---|
| 65 | REAL, INTENT(IN) :: pdtlw(ngrid,nlayer) ! (K/s) env |
|---|
| 66 | |
|---|
| 67 | ! input for second radiative transfer |
|---|
| 68 | logical, INTENT(IN) :: clearatm |
|---|
| 69 | INTEGER, INTENT(INOUT) :: icount |
|---|
| 70 | real, intent(in) :: zday |
|---|
| 71 | real, intent(in) :: zls |
|---|
| 72 | real, intent(in) :: tsurf(ngrid) |
|---|
| 73 | integer, intent(in) :: igout |
|---|
| 74 | real, intent(in) :: totstormfract(ngrid) |
|---|
| 75 | ! sbgrid scale water ice clouds |
|---|
| 76 | logical, intent(in) :: clearsky |
|---|
| 77 | real, intent(in) :: totcloudfrac(ngrid) |
|---|
| 78 | |
|---|
| 79 | !-------------------------------------------------------- |
|---|
| 80 | ! Output Variables |
|---|
| 81 | !-------------------------------------------------------- |
|---|
| 82 | |
|---|
| 83 | REAL, INTENT(OUT) :: pdqrds(ngrid,nlayer,nq) ! tendancy field for dust when detraining |
|---|
| 84 | REAL, INTENT(OUT) :: wspeed(ngrid,nlayer+1) ! vertical speed within the rocket dust storm |
|---|
| 85 | REAL, INTENT(OUT) :: dsodust(ngrid,nlayer) ! density scaled opacity of env. dust |
|---|
| 86 | REAL, INTENT(OUT) :: dsords(ngrid,nlayer) ! density scaled opacity of storm dust |
|---|
| 87 | |
|---|
| 88 | !-------------------------------------------------------- |
|---|
| 89 | ! Local variables |
|---|
| 90 | !-------------------------------------------------------- |
|---|
| 91 | INTEGER l,ig,tsub,iq,ll |
|---|
| 92 | ! chao local variables from callradite.F |
|---|
| 93 | REAL zdtlw1(ngrid,nlayer) ! (K/s) storm |
|---|
| 94 | REAL zdtsw1(ngrid,nlayer) ! (K/s) storm |
|---|
| 95 | REAL zt(ngrid,nlayer) ! actual temperature at mid-layer (K) |
|---|
| 96 | REAL zdtvert(nlayer) ! dT/dz , lapse rate |
|---|
| 97 | REAL ztlev(nlayer) ! temperature at intermediate levels l+1/2 without last level |
|---|
| 98 | |
|---|
| 99 | REAL zdtlw1_lev(nlayer),zdtsw1_lev(nlayer) ! rad. heating rate at intermediate levels l+1/2 |
|---|
| 100 | REAL zdtlw_lev(nlayer),zdtsw_lev(nlayer) ! rad. heating rate at intermediate levels l+1/2 |
|---|
| 101 | |
|---|
| 102 | REAL zqstorm_mass(ngrid,nlayer) ! tracer pq mass intermediate |
|---|
| 103 | REAL zqstorm_mass_col(nlayer) |
|---|
| 104 | REAL zqstorm_number(ngrid,nlayer) ! tracer field pq number intermediate |
|---|
| 105 | REAL zqstorm_number_col(nlayer) |
|---|
| 106 | |
|---|
| 107 | REAL zqi_mass(ngrid,nlayer) ! tracer pq mass intermediate |
|---|
| 108 | REAL zqi_number(ngrid,nlayer) ! tracer pq mass intermediate |
|---|
| 109 | REAL zdqvlstorm_mass(ngrid,nlayer) ! tendancy pdq mass after vertical transport |
|---|
| 110 | REAL zdqvlstorm_number(ngrid,nlayer) ! tendancy pdq number after vertical transport |
|---|
| 111 | REAL zdqdetstorm_mass(ngrid,nlayer) ! tendancy field pq mass after detrainment only |
|---|
| 112 | REAL zdqdetstorm_number(ngrid,nlayer) ! tendancy field pq number after detrainment only |
|---|
| 113 | |
|---|
| 114 | REAL zdqenv_mass(ngrid,nlayer) ! tendancy pdq mass from dust-> |
|---|
| 115 | ! stormdust in slp |
|---|
| 116 | REAL zdqenv_number(ngrid,nlayer) ! tendancy pdq number from dust-> |
|---|
| 117 | ! stormdust in slp |
|---|
| 118 | |
|---|
| 119 | REAL masse(nlayer) |
|---|
| 120 | REAL zq(ngrid,nlayer,nq) ! updated tracers |
|---|
| 121 | |
|---|
| 122 | REAL wrds(nlayer) ! vertical flux within the rocket dust storm |
|---|
| 123 | REAL wqrdsmass(nlayer+1) ! flux mass from vl_storm |
|---|
| 124 | REAL wqrdsnumber(nlayer+1) ! flux number from vl_storm |
|---|
| 125 | |
|---|
| 126 | INTEGER nsubtimestep !number of subtimestep when calling vl_storm |
|---|
| 127 | REAL subtimestep !ptimestep/nsubtimestep |
|---|
| 128 | REAL dtmax !considered time needed for dust to cross one layer |
|---|
| 129 | !minimal value over a column |
|---|
| 130 | logical storm(ngrid) !logical : true if you have some storm dust in the column |
|---|
| 131 | ! real slope(ngrid) !logical : true if you don't have storm and have |
|---|
| 132 | !a slope |
|---|
| 133 | ! real wslplev(ngrid,nlayer) |
|---|
| 134 | ! real wslp(ngrid,nlayer) |
|---|
| 135 | REAL coefdetrain !coefficient for detrainment : % of stormdust detrained |
|---|
| 136 | |
|---|
| 137 | REAL,PARAMETER:: coefmin =0.025 !C 0<c<1 Minimum fraction of stormdust detrained |
|---|
| 138 | REAL,PARAMETER:: detrainlim =0.1!0.25 !L stormdust detrained if wspeed < detrainlim |
|---|
| 139 | REAL,PARAMETER:: wlim =10. ! maximum vertical speed of rocket storms (m/s) |
|---|
| 140 | REAL,PARAMETER:: secu=3. !coefficient on wspeed to avoid dust crossing many layers during subtimestep |
|---|
| 141 | |
|---|
| 142 | ! terms for buoyancy and W^2 in equation: |
|---|
| 143 | ! w*dw/dz = k1*g*(T'-T)/T - k2*w^2 |
|---|
| 144 | real,parameter:: k1=0.00001 |
|---|
| 145 | real,parameter:: k2=0.25 |
|---|
| 146 | real,parameter:: mu0lim=0.1 |
|---|
| 147 | |
|---|
| 148 | ! diagnose |
|---|
| 149 | REAL detrainment(ngrid,nlayer) |
|---|
| 150 | real lapserate(ngrid,nlayer) |
|---|
| 151 | real deltahr(ngrid,nlayer+1) |
|---|
| 152 | real stormdust_m0(ngrid,nlayer) |
|---|
| 153 | real stormdust_m1(ngrid,nlayer) |
|---|
| 154 | real stormdust_m2(ngrid,nlayer) |
|---|
| 155 | |
|---|
| 156 | real hmax,hmin |
|---|
| 157 | ! real zh(ngrid,nlayer) |
|---|
| 158 | |
|---|
| 159 | logical,save :: firstcall=.true. |
|---|
| 160 | real alpha(ngrid) ! scale of the vertical motion (applicable for |
|---|
| 161 | ! rds and slp |
|---|
| 162 | ! variables for radiative transfer |
|---|
| 163 | REAL fluxsurf_lw1(ngrid) |
|---|
| 164 | REAL fluxsurf_sw1(ngrid,2) |
|---|
| 165 | REAL fluxtop_lw1(ngrid) |
|---|
| 166 | REAL fluxtop_sw1(ngrid,2) |
|---|
| 167 | REAL tauref(ngrid) |
|---|
| 168 | REAL tau(ngrid,naerkind) |
|---|
| 169 | REAL aerosol(ngrid,nlayer,naerkind) |
|---|
| 170 | REAL tauscaling(ngrid) |
|---|
| 171 | REAL taucloudtes(ngrid) |
|---|
| 172 | REAL rdust(ngrid,nlayer) |
|---|
| 173 | REAL rstormdust(ngrid,nlayer) |
|---|
| 174 | REAL rice(ngrid,nlayer) |
|---|
| 175 | REAL nuice(ngrid,nlayer) |
|---|
| 176 | |
|---|
| 177 | !variables related to slope,reference layer... |
|---|
| 178 | ! integer lref(ngrid),llref ! the reference layer of slopewind |
|---|
| 179 | ! real buoyt(nlayer) ! buoyancy term when there is a subgrid mountain |
|---|
| 180 | real slpbg(ngrid) ! temperature difference at half height of a mountain |
|---|
| 181 | |
|---|
| 182 | real zqdustslp(ngrid,nlayer),zndustslp(ngrid,nlayer) |
|---|
| 183 | real zqstormdustslp(ngrid,nlayer),znstormdustslp(ngrid,nlayer) |
|---|
| 184 | |
|---|
| 185 | real rdsdustqvl0(ngrid,nlayer) |
|---|
| 186 | real rdsdustqvl1(ngrid,nlayer) |
|---|
| 187 | ! real q2rds(ngrid,nlayer) |
|---|
| 188 | |
|---|
| 189 | !for second formule of wslp |
|---|
| 190 | real wtemp(ngrid,nlayer) ! a intermediate variable for wspeed |
|---|
| 191 | |
|---|
| 192 | !merge rds and slp |
|---|
| 193 | real newzt(ngrid,nlayer) !temperature with perturbation (integrated from |
|---|
| 194 | ! vetical motion) |
|---|
| 195 | real w0(ngrid) !prescribed slope winds at the first layer of atmosphere |
|---|
| 196 | ! real w1(ngrid) !prescribed slope winds at the first layer of atmosphere |
|---|
| 197 | ! real ztb1(ngrid) |
|---|
| 198 | real wadiabatic(ngrid,nlayer) !for diagnosis |
|---|
| 199 | ! real czt(nlayer),czlay(nlayer),czlev(nlayer+1) !temporary variables |
|---|
| 200 | real tnew(nlayer) !interpolated temperature profile above the top of Mons |
|---|
| 201 | real envtemp(nlayer) !interpolated background temperature profile |
|---|
| 202 | ! as the same height as tnew |
|---|
| 203 | real envt(ngrid,nlayer) ! output,for diagnosing |
|---|
| 204 | integer scheme(ngrid) ! diagnose |
|---|
| 205 | |
|---|
| 206 | ! Chao: for checking conservation of dust |
|---|
| 207 | ! real totdust0(ngrid) |
|---|
| 208 | ! real totdust1(ngrid) |
|---|
| 209 | |
|---|
| 210 | |
|---|
| 211 | ! ********************************************************************** |
|---|
| 212 | ! ********************************************************************** |
|---|
| 213 | ! Detached dust layers parametrization: two processes are included |
|---|
| 214 | ! 1) rocket dust storm |
|---|
| 215 | ! The radiative warming due to the presence of storm dust is |
|---|
| 216 | ! balanced by the adiabatic cooling. The tracer "storm dust" |
|---|
| 217 | ! is transported by the upward/downward flow. |
|---|
| 218 | ! 2) daytime slope winds |
|---|
| 219 | ! The daytime thermally driven upslope wind blows dust from the |
|---|
| 220 | ! bottom to the top of the mountain, upward flow keeps rising untill |
|---|
| 221 | ! the velocity becomes zero. Both the storm dust and environment dust |
|---|
| 222 | ! will be transported. |
|---|
| 223 | ! ********************************************************************** |
|---|
| 224 | ! ********************************************************************** |
|---|
| 225 | !! 1. Radiative transfer in storm dust |
|---|
| 226 | !! 2. Compute vertical velocity for storm dust |
|---|
| 227 | !! case 1 storm = false and nighttime: nothing to do |
|---|
| 228 | !! case 2 daytime slope wind scheme: (mu0(ig) > mu0lim and with storm=false) |
|---|
| 229 | !! case 3 rocket dust storm (storm=true) |
|---|
| 230 | !! 3. Vertical transport |
|---|
| 231 | !! 4. Detrainment |
|---|
| 232 | ! ********************************************************************** |
|---|
| 233 | ! ********************************************************************** |
|---|
| 234 | |
|---|
| 235 | |
|---|
| 236 | !! ********************************************************************** |
|---|
| 237 | !! Firstcall: Evaluate slope wind mesh fraction |
|---|
| 238 | IF (firstcall) then |
|---|
| 239 | call planetwide_maxval(hmons,hmax ) |
|---|
| 240 | call planetwide_minval(hmons,hmin ) |
|---|
| 241 | do ig=1,ngrid |
|---|
| 242 | ! It's hard to know the exact the scale of upward flow, |
|---|
| 243 | ! we assume that the maximun is 10% of the mesh area. |
|---|
| 244 | alpha_hmons(ig)= 0.1*(hmons(ig)-hmin)/(hmax-hmin) |
|---|
| 245 | enddo |
|---|
| 246 | firstcall = .false. |
|---|
| 247 | ENDIF !firstcall |
|---|
| 248 | |
|---|
| 249 | ! ********************************************************************** |
|---|
| 250 | ! 0. Initializations |
|---|
| 251 | ! ********************************************************************** |
|---|
| 252 | storm(:)=.false. |
|---|
| 253 | pdqrds(:,:,:) = 0. |
|---|
| 254 | zdqdetstorm_mass(:,:)=0. |
|---|
| 255 | zdqdetstorm_number(:,:)=0. |
|---|
| 256 | wspeed(:,:)=0. |
|---|
| 257 | detrainment(:,:)=0. |
|---|
| 258 | zqstorm_mass_col(:)=0. |
|---|
| 259 | zqstorm_number_col(:)=0. |
|---|
| 260 | lapserate(:,:)=0. |
|---|
| 261 | deltahr(:,:)=0. |
|---|
| 262 | rdsdustqvl0(:,:)=0. |
|---|
| 263 | rdsdustqvl1(:,:)=0. |
|---|
| 264 | zqstormdustslp(:,:)=0. |
|---|
| 265 | znstormdustslp(:,:)=0. |
|---|
| 266 | zqdustslp(:,:)=0. |
|---|
| 267 | zndustslp(:,:)=0. |
|---|
| 268 | zq(:,:,:) = 0. |
|---|
| 269 | |
|---|
| 270 | w0(:)=0. |
|---|
| 271 | ! w1(:)=0. |
|---|
| 272 | ! ztb1(:)=0. |
|---|
| 273 | newzt(:,:)=0 |
|---|
| 274 | wtemp(:,:)=0. |
|---|
| 275 | wadiabatic(:,:)=0. |
|---|
| 276 | slpbg(:)=0. |
|---|
| 277 | ! buoyt(:)=0. |
|---|
| 278 | tnew(:)=0. |
|---|
| 279 | envtemp(:)=0. |
|---|
| 280 | envt(:,:)=0. |
|---|
| 281 | scheme(:)=0 |
|---|
| 282 | alpha(:)=0. |
|---|
| 283 | stormdust_m0(:,:)=0. |
|---|
| 284 | stormdust_m1(:,:)=0. |
|---|
| 285 | stormdust_m2(:,:)=0. |
|---|
| 286 | ! totdust0(:)=0. |
|---|
| 287 | ! totdust1(:)=0. |
|---|
| 288 | |
|---|
| 289 | !! no update for the stormdust tracer and temperature fields |
|---|
| 290 | !! because previous callradite was for background dust |
|---|
| 291 | zq(1:ngrid,1:nlayer,1:nq)=pq(1:ngrid,1:nlayer,1:nq) |
|---|
| 292 | zt(1:ngrid,1:nlayer)=pt(1:ngrid,1:nlayer) |
|---|
| 293 | |
|---|
| 294 | !! get actual q for stormdust and dust tracers |
|---|
| 295 | do l=1,nlayer |
|---|
| 296 | do ig=1, ngrid |
|---|
| 297 | zqi_mass(ig,l)=zq(ig,l,igcm_dust_mass) |
|---|
| 298 | zqi_number(ig,l)=zq(ig,l,igcm_dust_number) |
|---|
| 299 | zqstorm_mass(ig,l)=zq(ig,l,igcm_stormdust_mass) |
|---|
| 300 | zqstorm_number(ig,l)=zq(ig,l,igcm_stormdust_number) |
|---|
| 301 | !for diagnostics: |
|---|
| 302 | stormdust_m0(ig,l)=zqstorm_mass(ig,l) |
|---|
| 303 | enddo |
|---|
| 304 | enddo ! of do l=1,nlayer |
|---|
| 305 | |
|---|
| 306 | !! Check if there is a rocket dust storm |
|---|
| 307 | do ig=1,ngrid |
|---|
| 308 | storm(ig)=.false. |
|---|
| 309 | do l=1,nlayer |
|---|
| 310 | if (zqstorm_mass(ig,l)/zqi_mass(ig,l) .gt. 1.E-4) then |
|---|
| 311 | storm(ig)=.true. |
|---|
| 312 | exit |
|---|
| 313 | endif |
|---|
| 314 | enddo |
|---|
| 315 | enddo |
|---|
| 316 | |
|---|
| 317 | ! ********************************************************************* |
|---|
| 318 | ! 1. Call the second radiative transfer for stormdust, obtain the extra heating |
|---|
| 319 | ! ********************************************************************* |
|---|
| 320 | CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo, & |
|---|
| 321 | emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, & |
|---|
| 322 | zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, & |
|---|
| 323 | fluxtop_sw1,tauref,tau,aerosol,dsodust,tauscaling, & |
|---|
| 324 | taucloudtes,rdust,rice,nuice,co2ice,rstormdust, & |
|---|
| 325 | totstormfract,clearatm,dsords, & |
|---|
| 326 | clearsky,totcloudfrac) |
|---|
| 327 | |
|---|
| 328 | ! ********************************************************************** |
|---|
| 329 | ! 2. Compute vertical velocity for storm dust |
|---|
| 330 | ! ********************************************************************** |
|---|
| 331 | DO ig=1,ngrid |
|---|
| 332 | !! ********************************************************************** |
|---|
| 333 | !! 2.1 case 1: Nothing to do when no storm and no slope, or |
|---|
| 334 | !! no storm and not daytime |
|---|
| 335 | IF (((alpha_hmons(ig) .EQ. 0.) .AND. .NOT.(storm(ig))) & |
|---|
| 336 | .OR. ((mu0(ig) .LE. mu0lim) .AND. .NOT.(storm(ig))) ) then |
|---|
| 337 | scheme(ig)=1 |
|---|
| 338 | cycle |
|---|
| 339 | endif |
|---|
| 340 | |
|---|
| 341 | |
|---|
| 342 | ! whatever the situation is, we need the vertical velocity computed by |
|---|
| 343 | ! the rds scheme |
|---|
| 344 | zdtvert(1)=0. !This is the env. lapse rate |
|---|
| 345 | DO l=1,nlayer-1 |
|---|
| 346 | zdtvert(l+1)=(zt(ig,l+1)-zt(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l)) |
|---|
| 347 | lapserate(ig,l+1)=zdtvert(l+1) !for diagnosing |
|---|
| 348 | ENDDO |
|---|
| 349 | |
|---|
| 350 | ! computing heating rates gradient at boundraies of each layer |
|---|
| 351 | ! start from surface |
|---|
| 352 | zdtlw1_lev(1)=0. |
|---|
| 353 | zdtsw1_lev(1)=0. |
|---|
| 354 | zdtlw_lev(1)=0. |
|---|
| 355 | zdtsw_lev(1)=0. |
|---|
| 356 | ztlev(1)=zt(ig,1) |
|---|
| 357 | |
|---|
| 358 | DO l=1,nlayer-1 |
|---|
| 359 | ! Calculation for the dust storm fraction |
|---|
| 360 | zdtlw1_lev(l+1)=(zdtlw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ & |
|---|
| 361 | zdtlw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / & |
|---|
| 362 | (pzlay(ig,l+1)-pzlay(ig,l)) |
|---|
| 363 | |
|---|
| 364 | zdtsw1_lev(l+1)=(zdtsw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ & |
|---|
| 365 | zdtsw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / & |
|---|
| 366 | (pzlay(ig,l+1)-pzlay(ig,l)) |
|---|
| 367 | !MV18: calculation for the background dust fraction |
|---|
| 368 | zdtlw_lev(l+1)=(pdtlw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ & |
|---|
| 369 | pdtlw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / & |
|---|
| 370 | (pzlay(ig,l+1)-pzlay(ig,l)) |
|---|
| 371 | |
|---|
| 372 | zdtsw_lev(l+1)=(pdtsw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ & |
|---|
| 373 | pdtsw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / & |
|---|
| 374 | (pzlay(ig,l+1)-pzlay(ig,l)) |
|---|
| 375 | |
|---|
| 376 | ztlev(l+1)=(zt(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ & |
|---|
| 377 | zt(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / & |
|---|
| 378 | (pzlay(ig,l+1)-pzlay(ig,l)) |
|---|
| 379 | ENDDO |
|---|
| 380 | |
|---|
| 381 | DO l=1,nlayer |
|---|
| 382 | deltahr(ig,l)=(zdtlw1_lev(l)+zdtsw1_lev(l)) & |
|---|
| 383 | -(zdtlw_lev(l)+zdtsw_lev(l)) |
|---|
| 384 | wadiabatic(ig,l)=-deltahr(ig,l)/(g/cpp+ & |
|---|
| 385 | max(zdtvert(l),-0.99*g/cpp)) |
|---|
| 386 | |
|---|
| 387 | !limit vertical wind in case of lapse rate close to adiabat |
|---|
| 388 | wadiabatic(ig,l)=max(wadiabatic(ig,l),-wlim) |
|---|
| 389 | wadiabatic(ig,l)=min(wadiabatic(ig,l),wlim) |
|---|
| 390 | ENDDO |
|---|
| 391 | |
|---|
| 392 | !! ********************************************************************** |
|---|
| 393 | !! 2.2 case 2: daytime slope wind scheme |
|---|
| 394 | IF ((mu0(ig) .gt. mu0lim) .and. .not. storm(ig)) then |
|---|
| 395 | scheme(ig)=2 |
|---|
| 396 | alpha(ig) = alpha_hmons(ig) |
|---|
| 397 | ! interpolate the env. temperature above the mountain top |
|---|
| 398 | call intep_vtemp(nlayer,hmons(ig),zt(ig,:),pzlay(ig,:), & |
|---|
| 399 | envtemp,slpbg(ig)) |
|---|
| 400 | envt(ig,:)=envtemp(:) !for diagnosis |
|---|
| 401 | |
|---|
| 402 | ! second: estimate the vertical velocity at boundraies of each layer |
|---|
| 403 | !wspeed(ig,1)=0. ! at surface, already initialized |
|---|
| 404 | |
|---|
| 405 | !!!!!!!the first layer of atmosphere!!!!!!!!!! |
|---|
| 406 | IF (slpbg(ig) .gt. 0.) THEN !only postive buoyancy |
|---|
| 407 | ! if slpbg(ig) lt 0, means the slope flow is colder than env. (night or |
|---|
| 408 | ! early morning ?) !!!!!!!!!!! |
|---|
| 409 | ! ideal method |
|---|
| 410 | !w1(ig)=-sqrt(g*slpbg(ig)/zt(ig,1)*hmons(ig))*sin(zsig(ig)) |
|---|
| 411 | ! new scheme, simply proportional to temperature and |
|---|
| 412 | ! mountain height |
|---|
| 413 | w0(ig)=-1.5e-4*g*slpbg(ig)/zt(ig,1)*hmons(ig) |
|---|
| 414 | ! otherwise, w0(ig) =0. |
|---|
| 415 | wspeed(ig,2)=w0(ig) |
|---|
| 416 | ELSE |
|---|
| 417 | wspeed(ig,2)=wadiabatic(ig,2) !! for early morning ? |
|---|
| 418 | ENDIF |
|---|
| 419 | |
|---|
| 420 | ! prepare the integration, NOTE: if w is too small, may have artificials |
|---|
| 421 | IF (abs(wspeed(ig,2)) .lt. 0.01 ) & |
|---|
| 422 | wspeed(ig,2)=sign(0.01,wspeed(ig,2)) |
|---|
| 423 | |
|---|
| 424 | newzt(ig,1)= zt(ig,1) !temperature of the first layer atmosphere |
|---|
| 425 | ! above the mountain top (radiative |
|---|
| 426 | ! equilibrium on Mars) |
|---|
| 427 | |
|---|
| 428 | ! estimate the vertical velocities |
|---|
| 429 | ! if w0(ig) >= 0, means downward motion, no upslope winds, we switch to |
|---|
| 430 | ! assume that the extra heating integrally convert to |
|---|
| 431 | ! vertical motion. |
|---|
| 432 | if ( w0(ig) .ge. 0 ) then !! normal, it is impossible, |
|---|
| 433 | !! because mu(ig)>0.1 here |
|---|
| 434 | do l=3,nlayer |
|---|
| 435 | wspeed(ig,l)=wadiabatic(ig,l) |
|---|
| 436 | enddo |
|---|
| 437 | else |
|---|
| 438 | ! estimate the velocities by taking into account the heating due |
|---|
| 439 | ! to storm dust, the cooling due to vertical motion ... |
|---|
| 440 | !!!!!!!!!!!the simple scheme!!!!!!!!! |
|---|
| 441 | do l=2,nlayer-1 |
|---|
| 442 | !if superadiabatic layer |
|---|
| 443 | if ( zdtvert(l) .lt. -g/cpp) then |
|---|
| 444 | !case 1 |
|---|
| 445 | ! test, also decrease adiabatically ? |
|---|
| 446 | !newzt(ig,l)= & |
|---|
| 447 | ! zt(ig,l-1)-g/cpp*(pzlay(ig,l)-pzlay(ig,l-1)) |
|---|
| 448 | newzt(ig,l)=zt(ig,l) |
|---|
| 449 | !wspeed(ig,l+1)=wspeed(ig,l) |
|---|
| 450 | else |
|---|
| 451 | !not superadiabatic |
|---|
| 452 | newzt(ig,l)=newzt(ig,l-1)+(deltahr(ig,l)/ & |
|---|
| 453 | (-wspeed(ig,l))-g/cpp)* & |
|---|
| 454 | (pzlay(ig,l)-pzlay(ig,l-1)) |
|---|
| 455 | endif ! end of if superadiabatic or not |
|---|
| 456 | |
|---|
| 457 | !wtemp(ig,l+1)=wspeed(ig,l)**2+2.*g*(pzlev(ig,l+1) & |
|---|
| 458 | ! -pzlev(ig,l))*(k1*(newzt(ig,l) & |
|---|
| 459 | ! -envtemp(l))/envtemp(l)) |
|---|
| 460 | wtemp(ig,l+1)=(1.-2.*k1*(pzlev(ig,l+1)-pzlev(ig,l)))*& |
|---|
| 461 | wspeed(ig,l)**2+2.*k2*g*(pzlev(ig,l+1) & |
|---|
| 462 | -pzlev(ig,l))*((newzt(ig,l) & |
|---|
| 463 | -envtemp(l))/envtemp(l)) |
|---|
| 464 | |
|---|
| 465 | if (wtemp(ig,l+1) .gt. 0.) then |
|---|
| 466 | !case 2 |
|---|
| 467 | wspeed(ig,l+1)=-sqrt(wtemp(ig,l+1)) |
|---|
| 468 | |
|---|
| 469 | ! if |wspeed| < |wadiabatic| then go to wadiabatic |
|---|
| 470 | if (wspeed(ig,l+1) .gt. wadiabatic(ig,l+1)) then |
|---|
| 471 | do ll=l,nlayer-1 |
|---|
| 472 | newzt(ig,ll)=envtemp(ll) |
|---|
| 473 | wspeed(ig,ll+1)=wadiabatic(ig,ll+1) |
|---|
| 474 | enddo |
|---|
| 475 | exit |
|---|
| 476 | endif |
|---|
| 477 | |
|---|
| 478 | ! avoid artificials |
|---|
| 479 | if (abs(wspeed(ig,l+1)) .lt. 0.01 ) & |
|---|
| 480 | wspeed(ig,l+1)=sign(0.01,wspeed(ig,l+1)) |
|---|
| 481 | |
|---|
| 482 | else if (l .lt. nlayer) then |
|---|
| 483 | !case 3 |
|---|
| 484 | do ll=l,nlayer-1 |
|---|
| 485 | newzt(ig,ll)=envtemp(ll) |
|---|
| 486 | wspeed(ig,ll+1)=wadiabatic(ig,ll+1) |
|---|
| 487 | enddo !overshot |
|---|
| 488 | exit |
|---|
| 489 | |
|---|
| 490 | else |
|---|
| 491 | wspeed(ig,l+1)=0. |
|---|
| 492 | endif |
|---|
| 493 | |
|---|
| 494 | enddo |
|---|
| 495 | |
|---|
| 496 | endif !w0 |
|---|
| 497 | |
|---|
| 498 | ELSE |
|---|
| 499 | !! ********************************************************************** |
|---|
| 500 | !! 2.3 case 3: storm=true |
|---|
| 501 | if (storm(ig)) then |
|---|
| 502 | scheme(ig)=3 |
|---|
| 503 | alpha(ig) = totstormfract(ig) |
|---|
| 504 | do l=1,nlayer |
|---|
| 505 | wspeed(ig,l)=wadiabatic(ig,l) |
|---|
| 506 | enddo |
|---|
| 507 | endif ! storm=1 |
|---|
| 508 | |
|---|
| 509 | ENDIF ! rds or slp |
|---|
| 510 | |
|---|
| 511 | |
|---|
| 512 | !!!!!!!! estimate the amount of dust for diagnostics |
|---|
| 513 | DO l=1,nlayer |
|---|
| 514 | ! transfer background dust + storm dust (concentrated) |
|---|
| 515 | zqstormdustslp(ig,l) =zqi_mass(ig,l)+ & |
|---|
| 516 | zqstorm_mass(ig,l)/alpha(ig) |
|---|
| 517 | znstormdustslp(ig,l) =zqi_number(ig,l)+ & |
|---|
| 518 | zqstorm_number(ig,l)/alpha(ig) |
|---|
| 519 | zqdustslp(ig,l) = zqi_mass(ig,l) |
|---|
| 520 | zndustslp(ig,l) = zqi_number(ig,l) |
|---|
| 521 | |
|---|
| 522 | rdsdustqvl0(ig,l)=zqstormdustslp(ig,l) !for diagnosis |
|---|
| 523 | ENDDO |
|---|
| 524 | |
|---|
| 525 | ! ********************************************************************** |
|---|
| 526 | ! 3. Vertical transport |
|---|
| 527 | ! ********************************************************************** |
|---|
| 528 | do l=1,nlayer |
|---|
| 529 | masse(l)=(pplev(ig,l)-pplev(ig,l+1))/g |
|---|
| 530 | enddo |
|---|
| 531 | !Estimation of "dtmax" (s) to be used for vertical transport |
|---|
| 532 | dtmax=ptimestep |
|---|
| 533 | !secu is a margin on subtimstep to avoid dust crossing many layers |
|---|
| 534 | do l=2,nlayer |
|---|
| 535 | if (wspeed(ig,l).lt.0.) then ! case up |
|---|
| 536 | dtmax=min(dtmax,(pzlev(ig,l)-pzlev(ig,l-1))/ & |
|---|
| 537 | (secu*abs(wspeed(ig,l)))) |
|---|
| 538 | else if (wspeed(ig,l).gt.0.) then |
|---|
| 539 | dtmax=min(dtmax,(pzlev(ig,l+1)-pzlev(ig,l))/ & |
|---|
| 540 | (secu*abs(wspeed(ig,l)))) |
|---|
| 541 | endif |
|---|
| 542 | enddo |
|---|
| 543 | |
|---|
| 544 | nsubtimestep= int(ptimestep/dtmax) |
|---|
| 545 | subtimestep=ptimestep/float(nsubtimestep) |
|---|
| 546 | |
|---|
| 547 | do l=1,nlayer |
|---|
| 548 | wrds(l)=wspeed(ig,l)*pplev(ig,l)/(r*ztlev(l)) & |
|---|
| 549 | *subtimestep |
|---|
| 550 | enddo |
|---|
| 551 | |
|---|
| 552 | do l=1,nlayer |
|---|
| 553 | zqstorm_mass_col(l)= zqstormdustslp(ig,l) !zqstorm_mass(ig,l) |
|---|
| 554 | zqstorm_number_col(l)=znstormdustslp(ig,l) ! zqstorm_number(ig,l) |
|---|
| 555 | enddo |
|---|
| 556 | |
|---|
| 557 | do tsub=1,nsubtimestep |
|---|
| 558 | wqrdsmass(:)=0. |
|---|
| 559 | wqrdsnumber(:)=0. |
|---|
| 560 | CALL vl_storm(nlayer,zqstorm_mass_col,2., & |
|---|
| 561 | masse,wrds ,wqrdsmass) |
|---|
| 562 | CALL vl_storm(nlayer,zqstorm_number_col,2., & |
|---|
| 563 | masse,wrds ,wqrdsnumber) |
|---|
| 564 | enddo |
|---|
| 565 | |
|---|
| 566 | !!!!!generate the "extra" dust |
|---|
| 567 | do l=1,nlayer |
|---|
| 568 | rdsdustqvl1(ig,l)=zqstorm_mass_col(l) ! for diagnosis |
|---|
| 569 | |
|---|
| 570 | ! extra dust = storm dust |
|---|
| 571 | !zqdustslp(ig,l)=zqi_mass(ig,l) !(1.-alpha(ig))*zqi_mass(ig,l) |
|---|
| 572 | !zndustslp(ig,l)=zqi_number(ig,l) !(1.-alpha(ig))*zqi_number(ig,l) |
|---|
| 573 | !zqstorm_mass_col(l)=alpha(ig)*zqstorm_mass_col(l) |
|---|
| 574 | !zqstorm_number_col(l)=alpha(ig)*zqstorm_number_col(l) |
|---|
| 575 | |
|---|
| 576 | !with compensation |
|---|
| 577 | if (zqstorm_mass_col(l) .lt. zqi_mass(ig,l) ) then |
|---|
| 578 | ! the following two equations are easier to understand |
|---|
| 579 | zqdustslp(ig,l)=(1.-alpha(ig))*zqi_mass(ig,l)+alpha(ig)* & |
|---|
| 580 | zqstorm_mass_col(l) |
|---|
| 581 | zndustslp(ig,l)=(1.-alpha(ig))*zqi_number(ig,l)+alpha(ig)& |
|---|
| 582 | *zqstorm_number_col(l) |
|---|
| 583 | !with a bug, should be zqi+alpha**** |
|---|
| 584 | !zqdustslp(ig,l)=zqi_mass(ig,l)-alpha(ig)* & |
|---|
| 585 | ! (zqstorm_mass_col(l)-zqi_mass(ig,l)) |
|---|
| 586 | !zndustslp(ig,l)=zqi_number(ig,l)-alpha(ig)* & |
|---|
| 587 | ! (zqstorm_number_col(l)-zqi_number(ig,l)) |
|---|
| 588 | zqstorm_mass_col(l)=0. |
|---|
| 589 | zqstorm_number_col(l)=0. |
|---|
| 590 | else |
|---|
| 591 | zqstorm_mass_col(l)=alpha(ig)* & |
|---|
| 592 | (zqstorm_mass_col(l)-zqi_mass(ig,l)) |
|---|
| 593 | zqstorm_number_col(l)=alpha(ig)* & |
|---|
| 594 | (zqstorm_number_col(l)-zqi_number(ig,l)) |
|---|
| 595 | ! the mass mixing ratio of environmental dust doesn't change. |
|---|
| 596 | endif |
|---|
| 597 | !diagnose |
|---|
| 598 | stormdust_m1(ig,l)=zqstorm_mass_col(l) |
|---|
| 599 | enddo |
|---|
| 600 | |
|---|
| 601 | !======================================================================= |
|---|
| 602 | ! calculate the tendencies due to vertical transport |
|---|
| 603 | do l=1,nlayer |
|---|
| 604 | ! tendencies due to vertical transport |
|---|
| 605 | zdqvlstorm_mass(ig,l)= (zqstorm_mass_col(l)- & |
|---|
| 606 | zqstorm_mass(ig,l)) /ptimestep |
|---|
| 607 | zdqvlstorm_number(ig,l)= (zqstorm_number_col(l)- & |
|---|
| 608 | zqstorm_number(ig,l)) /ptimestep |
|---|
| 609 | |
|---|
| 610 | zdqenv_mass(ig,l)=(zqdustslp(ig,l)-zqi_mass(ig,l))/ptimestep |
|---|
| 611 | zdqenv_number(ig,l)=(zndustslp(ig,l)-zqi_number(ig,l)) & |
|---|
| 612 | /ptimestep |
|---|
| 613 | |
|---|
| 614 | ! chao for output only |
|---|
| 615 | !qstormdustvl1(ig,l)=zqstorm_mass_col(l) |
|---|
| 616 | !nstormdustvl1(ig,l)=zqstorm_number_col(l) |
|---|
| 617 | !stormdust_m_col1(ig)=stormdust_m_col1(ig)+zqstorm_mass_col(l) & |
|---|
| 618 | ! *(pplev(ig,l)-pplev(ig,l+1))/g |
|---|
| 619 | !rdsdustqvl1(ig,l)=zqstorm_mass_col(l) |
|---|
| 620 | enddo |
|---|
| 621 | |
|---|
| 622 | ! ********************************************************************** |
|---|
| 623 | ! 4. Detrainment: convert dust storm to background dust |
|---|
| 624 | ! ********************************************************************** |
|---|
| 625 | do l=1,nlayer |
|---|
| 626 | ! compute the coefficient of detrainment |
|---|
| 627 | if ((max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1))) .lt. & |
|---|
| 628 | detrainlim) .or. (zqdustslp(ig,l) .gt. & |
|---|
| 629 | 10000.*zqstorm_mass_col(l))) then |
|---|
| 630 | coefdetrain=1. |
|---|
| 631 | else if (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1))) & |
|---|
| 632 | .le. wlim) then |
|---|
| 633 | ! case where detrainment depends on vertical wind |
|---|
| 634 | ! coefdetrain=0.5*(((1-coefmin)/(detrainlim-3.)**2)* & |
|---|
| 635 | ! (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1)))-3.)**2 & |
|---|
| 636 | ! +coefmin) |
|---|
| 637 | coefdetrain=1.*(((1-coefmin)/(detrainlim-wlim)**2)* & |
|---|
| 638 | (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1)))-wlim)**2 & |
|---|
| 639 | +coefmin) |
|---|
| 640 | !coefdetrain=0.5 |
|---|
| 641 | else if (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1))).gt. 10. )& |
|---|
| 642 | then |
|---|
| 643 | coefdetrain=0.025 |
|---|
| 644 | else |
|---|
| 645 | coefdetrain=coefmin |
|---|
| 646 | endif |
|---|
| 647 | |
|---|
| 648 | detrainment(ig,l)=coefdetrain !for diagnose |
|---|
| 649 | |
|---|
| 650 | ! Calculate tendancies corresponding to pdq after detrainement |
|---|
| 651 | ! pdqdet = tendancy corresponding to detrainment only |
|---|
| 652 | zdqdetstorm_mass(ig,l)=-coefdetrain*zqstorm_mass_col(l) & |
|---|
| 653 | /ptimestep |
|---|
| 654 | zdqdetstorm_number(ig,l)=-coefdetrain*zqstorm_number_col(l) & |
|---|
| 655 | /ptimestep |
|---|
| 656 | |
|---|
| 657 | ! pdqrds ( tendancy corresponding to vertical transport and |
|---|
| 658 | ! detrainment) = zdqvlstorm + pdqdet |
|---|
| 659 | pdqrds(ig,l,igcm_stormdust_mass)=zdqdetstorm_mass(ig,l) & |
|---|
| 660 | +zdqvlstorm_mass(ig,l) |
|---|
| 661 | pdqrds(ig,l,igcm_stormdust_number)=zdqdetstorm_number(ig,l) & |
|---|
| 662 | +zdqvlstorm_number(ig,l) |
|---|
| 663 | pdqrds(ig,l,igcm_dust_mass)= zdqenv_mass(ig,l) & |
|---|
| 664 | -zdqdetstorm_mass(ig,l) |
|---|
| 665 | pdqrds(ig,l,igcm_dust_number)= zdqenv_number(ig,l) & |
|---|
| 666 | -zdqdetstorm_number(ig,l) |
|---|
| 667 | |
|---|
| 668 | !diagnose |
|---|
| 669 | stormdust_m2(ig,l)=zqstorm_mass_col(l)-coefdetrain*zqstorm_mass_col(l) |
|---|
| 670 | enddo ! nlayer |
|---|
| 671 | ! endif |
|---|
| 672 | !======================================================================= |
|---|
| 673 | enddo ! end do ig=1,ngrid |
|---|
| 674 | |
|---|
| 675 | ! !chao check conservation here |
|---|
| 676 | ! do l=1,nlayer |
|---|
| 677 | ! do ig=1,ngrid |
|---|
| 678 | ! totdust0(ig)=totdust0(ig)+ & |
|---|
| 679 | ! zq(ig,l,igcm_stormdust_mass)* & |
|---|
| 680 | ! ((pplev(ig,l) - pplev(ig,l+1)) / g) & |
|---|
| 681 | ! + zq(ig,l,igcm_dust_mass)* & |
|---|
| 682 | ! ((pplev(ig,l) - pplev(ig,l+1)) / g) |
|---|
| 683 | |
|---|
| 684 | ! totdust1(ig)=totdust1(ig)+ & |
|---|
| 685 | ! (zq(ig,l,igcm_stormdust_mass) + & |
|---|
| 686 | ! pdqrds(ig,l,igcm_stormdust_mass)*ptimestep)* & |
|---|
| 687 | ! ((pplev(ig,l) - pplev(ig,l+1)) / g) & |
|---|
| 688 | ! + ( zq(ig,l,igcm_dust_mass)+ & |
|---|
| 689 | ! pdqrds(ig,l,igcm_dust_mass)*ptimestep)* & |
|---|
| 690 | ! ((pplev(ig,l) - pplev(ig,l+1)) / g) |
|---|
| 691 | ! enddo |
|---|
| 692 | ! enddo |
|---|
| 693 | |
|---|
| 694 | ! call writediagfi(ngrid,'totdust0','total dust before rds', & |
|---|
| 695 | ! ' ',2,totdust0) |
|---|
| 696 | ! call writediagfi(ngrid,'totdust1','total dust after rds', & |
|---|
| 697 | ! ' ',2,totdust1) |
|---|
| 698 | !output for diagnosis |
|---|
| 699 | call WRITEDIAGFI(ngrid,'detrainment', & |
|---|
| 700 | 'coefficient of detrainment',' ',3,detrainment) |
|---|
| 701 | !call WRITEDIAGFI(ngrid,'qstormvl1','mmr of stormdust after rds_vl', & |
|---|
| 702 | ! & 'kg/kg',3,qstormdustvl1) |
|---|
| 703 | call WRITEDIAGFI(ngrid,'lapserate','lapse rate in the storm', & |
|---|
| 704 | & 'k/m',3,lapserate) |
|---|
| 705 | call WRITEDIAGFI(ngrid,'deltahr','extra heating rates', & |
|---|
| 706 | & 'K/s',3,deltahr) |
|---|
| 707 | call WRITEDIAGFI(ngrid,'wold', & |
|---|
| 708 | 'wind generated from adiabatic colling ', & |
|---|
| 709 | & 'm/s',3,wadiabatic) |
|---|
| 710 | call WRITEDIAGFI(ngrid,'newzt','perturbated temperature', & |
|---|
| 711 | & 'K/s',3,newzt) |
|---|
| 712 | call WRITEDIAGFI(ngrid,'zt','unperturbated temperature', & |
|---|
| 713 | & 'K/s',3,zt) |
|---|
| 714 | call WRITEDIAGFI(ngrid,'wtemp','under square root', & |
|---|
| 715 | & 'K/s',3,wtemp) |
|---|
| 716 | !call WRITEDIAGFI(ngrid,'stormdust_m_col1','mass of stormdust after rds_vl', & |
|---|
| 717 | ! & 'kg/kg',2,stormdust_m_col1) |
|---|
| 718 | !call WRITEDIAGFI(ngrid,'temprds','temp for calculating zdtvert', & |
|---|
| 719 | ! & 'k',3,temprds) |
|---|
| 720 | call WRITEDIAGFI(ngrid,'stormdust_m0','mass col of stormdust before rds_vl', & |
|---|
| 721 | & 'kg/kg',3,stormdust_m0) |
|---|
| 722 | call WRITEDIAGFI(ngrid,'stormdust_m1','mass col of stormdust after rds_vl', & |
|---|
| 723 | & 'kg/kg',3,stormdust_m1) |
|---|
| 724 | call WRITEDIAGFI(ngrid,'stormdust_m2','mass col of stormdust after rds_vl', & |
|---|
| 725 | & 'kg/kg',3,stormdust_m2) |
|---|
| 726 | |
|---|
| 727 | ! call writediagfi(ngrid,'wslp','estimated slope winds','m/s',3,wslp) |
|---|
| 728 | ! call writediagfi(ngrid,'wslp2','estimated slope winds2','m/s',3,wslp2) |
|---|
| 729 | ! call writediagfi(ngrid,'zhb','estimated slope winds2','m/s',3,zhb) |
|---|
| 730 | ! call writediagfi(ngrid,'bruntf','bouyancy frequency',' ',3,bruntf) |
|---|
| 731 | ! call writediagfi(ngrid,'slpdepth','slope depth','m',2,slpdepth) |
|---|
| 732 | ! call writediagfi(ngrid,'slpu','perbulation u','m/s',3,slpu) |
|---|
| 733 | ! call writediagfi(ngrid,'slpzh','perbulation zh',' ',3,slpzh) |
|---|
| 734 | ! call writediagfi(ngrid,'zqslp','zq in rocketduststorm','ikg/kg',3, & |
|---|
| 735 | ! zq(:,:,igcm_dust_mass)) |
|---|
| 736 | ! call writediagfi(ngrid,'zrdsqslp','zq in rocketduststorm','ikg/kg',3, & |
|---|
| 737 | ! zq(:,:,igcm_stormdust_mass)) |
|---|
| 738 | ! call writediagfi(ngrid,'wslplev','estimated slope winds','m/s',3,wslplev) |
|---|
| 739 | ! call writediagfi(ngrid,'slope','identified slope wind effect',' ',2,slope) |
|---|
| 740 | call writediagfi(ngrid,'w0','max of slope wind',' ',2,w0) |
|---|
| 741 | ! call writediagfi(ngrid,'w1','max of slope wind',' ',2,w1) |
|---|
| 742 | call writediagfi(ngrid,'mu0','cosine of solar incidence angle',& |
|---|
| 743 | ' ',2,mu0) |
|---|
| 744 | ! call writediagfi(ngrid,'storm','identified rocket dust storm',& |
|---|
| 745 | ! ' ',2,real(storm)) |
|---|
| 746 | call writediagfi(ngrid,'scheme','which scheme',& |
|---|
| 747 | ' ',2,real(scheme)) |
|---|
| 748 | call writediagfi(ngrid,'alpha','coefficient alpha',' ',2,alpha) |
|---|
| 749 | ! call writediagfi(ngrid,'q2rds','alpha zq',' ', & |
|---|
| 750 | ! 3,q2rds) |
|---|
| 751 | call writediagfi(ngrid,'rdsdustqvl0','not vl storm slp', & |
|---|
| 752 | 'kg/kg',3,zqstormdustslp) |
|---|
| 753 | call writediagfi(ngrid,'rdsdustqvl1','vled storm slp', & |
|---|
| 754 | 'kg/kg',3,rdsdustqvl1) |
|---|
| 755 | call writediagfi(ngrid,'dustqvl0','not vl slp', & |
|---|
| 756 | 'kg/kg',3,zqi_mass) |
|---|
| 757 | call writediagfi(ngrid,'dustqvl1','vled slp', & |
|---|
| 758 | 'kg/kg',3,zqdustslp) |
|---|
| 759 | ! call WRITEDIAGFI(ngrid,'lmax_th2', & |
|---|
| 760 | ! 'hauteur du thermique','point', & |
|---|
| 761 | ! 2,real(lmax_th(:))) |
|---|
| 762 | ! call WRITEDIAGFI(ngrid,'zmax_th2', & |
|---|
| 763 | ! 'hauteur du thermique','m', & |
|---|
| 764 | ! 2,zmax_th) |
|---|
| 765 | ! call writediagfi(ngrid,'lslope','lenght of slope',' ',2,lslope) |
|---|
| 766 | ! call writediagfi(ngrid,'hmon','identified slope wind effect',' ',2,hmon) |
|---|
| 767 | call writediagfi(ngrid,'envt','interpolated env. temp.', & |
|---|
| 768 | 'K',3,envt) |
|---|
| 769 | call writediagfi(ngrid,'hmons','identified slope wind effect', & |
|---|
| 770 | ' ',2,hmons) |
|---|
| 771 | ! call writediagfi(ngrid,'slpbg','temp. diff along slope', & |
|---|
| 772 | ! ' ',2,slpbg) |
|---|
| 773 | ! call writediagfi(ngrid,'zmea','identified slope wind effect', & |
|---|
| 774 | ! ' ',2,zmea) |
|---|
| 775 | ! call writediagfi(ngrid,'zsig','identified slope wind effect', & |
|---|
| 776 | ! ' ',2,zsig) |
|---|
| 777 | ! call writediagfi(ngrid,'zhslpenv','difference of zh above mons',' ',2,zhslpenv) |
|---|
| 778 | ! call writediagfi(ngrid,'lref','identified slope wind effect',' ',2,real(lref)) |
|---|
| 779 | END SUBROUTINE rocketduststorm |
|---|
| 780 | |
|---|
| 781 | !******************************************************************************** |
|---|
| 782 | !******************************************************************************** |
|---|
| 783 | SUBROUTINE vl_storm(nlay,q,pente_max,masse,w,wq) |
|---|
| 784 | ! |
|---|
| 785 | ! Auteurs: P.Le Van, F.Hourdin, F.Forget |
|---|
| 786 | ! |
|---|
| 787 | ! ******************************************************************** |
|---|
| 788 | ! Shema d'advection " pseudo amont " dans la verticale |
|---|
| 789 | ! pour appel dans la physique (sedimentation) |
|---|
| 790 | ! ******************************************************************** |
|---|
| 791 | ! q rapport de melange (kg/kg)... |
|---|
| 792 | ! masse : masse de la couche Dp/g |
|---|
| 793 | ! w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2) |
|---|
| 794 | ! pente_max = 2 conseillee |
|---|
| 795 | ! |
|---|
| 796 | ! |
|---|
| 797 | ! -------------------------------------------------------------------- |
|---|
| 798 | IMPLICIT NONE |
|---|
| 799 | ! |
|---|
| 800 | |
|---|
| 801 | ! Arguments: |
|---|
| 802 | ! ---------- |
|---|
| 803 | integer,intent(in) :: nlay ! number of atmospheric layers |
|---|
| 804 | real masse(nlay),pente_max |
|---|
| 805 | REAL q(nlay) |
|---|
| 806 | REAL w(nlay) |
|---|
| 807 | REAL wq(nlay+1) |
|---|
| 808 | ! |
|---|
| 809 | ! Local |
|---|
| 810 | ! --------- |
|---|
| 811 | ! |
|---|
| 812 | INTEGER i,l,j,ii |
|---|
| 813 | ! |
|---|
| 814 | |
|---|
| 815 | real dzq(nlay),dzqw(nlay),adzqw(nlay),dzqmax |
|---|
| 816 | real newmasse |
|---|
| 817 | real sigw, Mtot, MQtot |
|---|
| 818 | integer m |
|---|
| 819 | |
|---|
| 820 | REAL SSUM,CVMGP,CVMGT |
|---|
| 821 | integer ismax,ismin |
|---|
| 822 | |
|---|
| 823 | |
|---|
| 824 | ! On oriente tout dans le sens de la pression c'est a dire dans le |
|---|
| 825 | ! sens de W |
|---|
| 826 | |
|---|
| 827 | do l=2,nlay |
|---|
| 828 | dzqw(l)=q(l-1)-q(l) |
|---|
| 829 | adzqw(l)=abs(dzqw(l)) |
|---|
| 830 | enddo |
|---|
| 831 | |
|---|
| 832 | do l=2,nlay-1 |
|---|
| 833 | #ifdef CRAY |
|---|
| 834 | dzq(l)=0.5* |
|---|
| 835 | , cvmgp(dzqw(l)+dzqw(l+1),0.,dzqw(l)*dzqw(l+1)) |
|---|
| 836 | #else |
|---|
| 837 | if(dzqw(l)*dzqw(l+1).gt.0.) then |
|---|
| 838 | dzq(l)=0.5*(dzqw(l)+dzqw(l+1)) |
|---|
| 839 | else |
|---|
| 840 | dzq(l)=0. |
|---|
| 841 | endif |
|---|
| 842 | #endif |
|---|
| 843 | dzqmax=pente_max*min(adzqw(l),adzqw(l+1)) |
|---|
| 844 | dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l)) |
|---|
| 845 | enddo |
|---|
| 846 | |
|---|
| 847 | dzq(1)=0. |
|---|
| 848 | dzq(nlay)=0. |
|---|
| 849 | |
|---|
| 850 | ! write(*,*),'TB14 wq before up',wq(1,:) |
|---|
| 851 | ! write(*,*),'TB14 q before up',q(1,:) |
|---|
| 852 | ! --------------------------------------------------------------- |
|---|
| 853 | ! .... calcul des termes d'advection verticale ....... |
|---|
| 854 | ! --------------------------------------------------------------- |
|---|
| 855 | |
|---|
| 856 | ! calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq |
|---|
| 857 | ! |
|---|
| 858 | ! No flux at the model top: |
|---|
| 859 | wq(nlay+1)=0. |
|---|
| 860 | |
|---|
| 861 | ! 1) Compute wq where w < 0 (up) (NOT USEFUL FOR SEDIMENTATION) |
|---|
| 862 | ! =============================== |
|---|
| 863 | |
|---|
| 864 | ! Surface flux up: |
|---|
| 865 | if(w(1).lt.0.) wq(1)=0. ! warning : not always valid |
|---|
| 866 | |
|---|
| 867 | do l = 1,nlay-1 ! loop different than when w>0 |
|---|
| 868 | if(w(l+1).le.0)then |
|---|
| 869 | |
|---|
| 870 | ! Regular scheme (transfered mass < 1 layer) |
|---|
| 871 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 872 | if(-w(l+1).le.masse(l))then |
|---|
| 873 | sigw=w(l+1)/masse(l) |
|---|
| 874 | wq(l+1)=w(l+1)*(q(l)-0.5*(1.+sigw)*dzq(l)) |
|---|
| 875 | ! Extended scheme (transfered mass > 1 layer) |
|---|
| 876 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 877 | else |
|---|
| 878 | m = l-1 |
|---|
| 879 | Mtot = masse(m+1) |
|---|
| 880 | MQtot = masse(m+1)*q(m+1) |
|---|
| 881 | if (m.le.0)goto 77 |
|---|
| 882 | do while(-w(l+1).gt.(Mtot+masse(m))) |
|---|
| 883 | ! do while(-w(l+1).gt.Mtot) |
|---|
| 884 | m=m-1 |
|---|
| 885 | Mtot = Mtot + masse(m+1) |
|---|
| 886 | MQtot = MQtot + masse(m+1)*q(m+1) |
|---|
| 887 | if (m.le.0)goto 77 |
|---|
| 888 | end do |
|---|
| 889 | 77 continue |
|---|
| 890 | |
|---|
| 891 | if (m.gt.0) then |
|---|
| 892 | sigw=(w(l+1)+Mtot)/masse(m) |
|---|
| 893 | wq(l+1)= (MQtot + (-w(l+1)-Mtot)* & |
|---|
| 894 | (q(m)-0.5*(1.+sigw)*dzq(m)) ) |
|---|
| 895 | else |
|---|
| 896 | ! new |
|---|
| 897 | w(l+1) = -Mtot |
|---|
| 898 | wq(l+1) = -MQtot |
|---|
| 899 | ! write(*,*) 'TB14 MQtot = ',MQtot,l |
|---|
| 900 | |
|---|
| 901 | ! old |
|---|
| 902 | ! wq(l+1)= (MQtot + (-w(l+1)-Mtot)*qm(1)) |
|---|
| 903 | ! write(*,*) 'a rather weird situation in vlz_fi !' |
|---|
| 904 | ! stop |
|---|
| 905 | end if |
|---|
| 906 | endif |
|---|
| 907 | endif ! w<0 (up) |
|---|
| 908 | enddo |
|---|
| 909 | |
|---|
| 910 | do l = 1,nlay-1 ! loop different than when w>0 |
|---|
| 911 | |
|---|
| 912 | q(l)=q(l) + (wq(l+1)-wq(l))/masse(l) |
|---|
| 913 | |
|---|
| 914 | enddo |
|---|
| 915 | |
|---|
| 916 | ! write(*,*),'TB14 masse',masse(1,:) |
|---|
| 917 | ! write(*,*),'TB14 wq before down after up',wq(1,:) |
|---|
| 918 | ! write(*,*),'TB14 q before down',q(1,:) |
|---|
| 919 | |
|---|
| 920 | ! 2) Compute wq where w > 0 (down) (ALWAYS FOR SEDIMENTATION) |
|---|
| 921 | ! =============================== |
|---|
| 922 | |
|---|
| 923 | ! Initialisation wq = 0 to consider now only downward flux |
|---|
| 924 | wq(:)=0. ! |
|---|
| 925 | |
|---|
| 926 | do l = 1,nlay ! loop different than when w<0 |
|---|
| 927 | |
|---|
| 928 | if(w(l).gt.0.)then |
|---|
| 929 | |
|---|
| 930 | ! Regular scheme (transfered mass < 1 layer) |
|---|
| 931 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 932 | if(w(l).le.masse(l))then |
|---|
| 933 | sigw=w(l)/masse(l) |
|---|
| 934 | wq(l)=w(l)*(q(l)+0.5*(1.-sigw)*dzq(l)) |
|---|
| 935 | ! write(*,*),'TB14 wq after up',wq(1,:) |
|---|
| 936 | |
|---|
| 937 | |
|---|
| 938 | ! Extended scheme (transfered mass > 1 layer) |
|---|
| 939 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 940 | else |
|---|
| 941 | m=l |
|---|
| 942 | Mtot = masse(m) |
|---|
| 943 | MQtot = masse(m)*q(m) |
|---|
| 944 | if(m.ge.nlay)goto 88 |
|---|
| 945 | do while(w(l).gt.(Mtot+masse(m+1))) |
|---|
| 946 | m=m+1 |
|---|
| 947 | Mtot = Mtot + masse(m) |
|---|
| 948 | MQtot = MQtot + masse(m)*q(m) |
|---|
| 949 | if(m.ge.nlay)goto 88 |
|---|
| 950 | end do |
|---|
| 951 | 88 continue |
|---|
| 952 | if (m.lt.nlay) then |
|---|
| 953 | sigw=(w(l)-Mtot)/masse(m+1) |
|---|
| 954 | wq(l)=(MQtot + (w(l)-Mtot)* & |
|---|
| 955 | (q(m+1)+0.5*(1.-sigw)*dzq(m+1)) ) |
|---|
| 956 | else |
|---|
| 957 | w(l) = Mtot |
|---|
| 958 | wq(l) = MQtot |
|---|
| 959 | end if |
|---|
| 960 | end if |
|---|
| 961 | end if ! w>0 (down) |
|---|
| 962 | enddo |
|---|
| 963 | |
|---|
| 964 | do l = 1,nlay ! loop different than when w<0 |
|---|
| 965 | |
|---|
| 966 | q(l)=q(l) + (wq(l+1)-wq(l))/masse(l) |
|---|
| 967 | |
|---|
| 968 | enddo |
|---|
| 969 | end subroutine vl_storm |
|---|
| 970 | !******************************************************************************** |
|---|
| 971 | SUBROUTINE intep_vtemp(nlayer,hm,temp,zlay,envtemp,slpb) |
|---|
| 972 | |
|---|
| 973 | USE comcstfi_h, only: g,cpp |
|---|
| 974 | |
|---|
| 975 | implicit none |
|---|
| 976 | ! this subroutine is using for obtaining the environmental |
|---|
| 977 | ! temperature profile when a subgrid mountain exists. |
|---|
| 978 | |
|---|
| 979 | |
|---|
| 980 | ! input: |
|---|
| 981 | integer,intent(in) :: nlayer |
|---|
| 982 | real,intent(in) :: hm ! the height of mountain |
|---|
| 983 | real,intent(in) :: temp(nlayer) !large scale temp. profile of one mesh |
|---|
| 984 | real,intent(in) :: zlay(nlayer) ! height at the middle of each layer |
|---|
| 985 | |
|---|
| 986 | ! output: |
|---|
| 987 | real,intent(out) :: envtemp(nlayer) ! environment temperature |
|---|
| 988 | real,intent(out) :: slpb !the temperature difference between slope and |
|---|
| 989 | ! env. at the half height of mountain |
|---|
| 990 | |
|---|
| 991 | ! local variables |
|---|
| 992 | integer l,il,tmpl |
|---|
| 993 | integer lnew !the layer of atmosphere above the mountain |
|---|
| 994 | ! corresponding to the env. (for buoyancy |
|---|
| 995 | ! calc. ) |
|---|
| 996 | real newh(nlayer) !height at the middle of each layer |
|---|
| 997 | ! account for the exist of mountain |
|---|
| 998 | ! real g,cpp |
|---|
| 999 | real halfh ! half the height of a mountain |
|---|
| 1000 | |
|---|
| 1001 | !initilize the array |
|---|
| 1002 | lnew=0 |
|---|
| 1003 | newh(:)=0. |
|---|
| 1004 | envtemp(:)=0. |
|---|
| 1005 | tmpl=1 |
|---|
| 1006 | |
|---|
| 1007 | do l=1,nlayer |
|---|
| 1008 | newh(l)=hm+zlay(l) |
|---|
| 1009 | |
|---|
| 1010 | do il=tmpl,nlayer-1 !MV18: added the -1 |
|---|
| 1011 | if (newh(l) .ge. zlay(il) .and. newh(l) .lt. zlay(il+1))then |
|---|
| 1012 | ! find the corresponding layer |
|---|
| 1013 | lnew=il |
|---|
| 1014 | |
|---|
| 1015 | ! interpolate |
|---|
| 1016 | envtemp(l) = temp(il)+(newh(l)-zlay(lnew))*& |
|---|
| 1017 | (temp(il+1)-temp(il))/(zlay(il+1)-zlay(il)) |
|---|
| 1018 | |
|---|
| 1019 | exit !go to the next layer |
|---|
| 1020 | else if (newh(l) .ge. zlay(nlayer)) then |
|---|
| 1021 | ! higher than the highest layer |
|---|
| 1022 | lnew=nlayer |
|---|
| 1023 | envtemp(l)=temp(nlayer) |
|---|
| 1024 | |
|---|
| 1025 | endif |
|---|
| 1026 | enddo |
|---|
| 1027 | ! this can accelerate the loop |
|---|
| 1028 | tmpl=lnew |
|---|
| 1029 | |
|---|
| 1030 | enddo |
|---|
| 1031 | |
|---|
| 1032 | halfh=0.5*hm |
|---|
| 1033 | if (halfh .le. zlay(1) ) then |
|---|
| 1034 | slpb=0. |
|---|
| 1035 | else if (halfh .gt. zlay(nlayer)) then |
|---|
| 1036 | !normally, impossible for halfh gt zlay(l), anyway... |
|---|
| 1037 | tmpl=nlayer |
|---|
| 1038 | !difference between surface and atmosphere (env.) |
|---|
| 1039 | slpb=temp(1)-(temp(nlayer-1)+((halfh-zlay(nlayer-1))* & |
|---|
| 1040 | (temp(tmpl)-temp(tmpl-1))/(zlay(tmpl)-zlay(tmpl-1)))) |
|---|
| 1041 | else |
|---|
| 1042 | do l=1,nlayer-1 |
|---|
| 1043 | if ((halfh .gt. zlay(l)) .and. (halfh .le. zlay(l+1)))then |
|---|
| 1044 | tmpl= l |
|---|
| 1045 | slpb=temp(1)-(temp(tmpl)+(halfh-zlay(tmpl))* & |
|---|
| 1046 | (temp(tmpl+1)-temp(tmpl))/(zlay(tmpl+1)-zlay(tmpl))) |
|---|
| 1047 | endif |
|---|
| 1048 | enddo |
|---|
| 1049 | endif |
|---|
| 1050 | end subroutine intep_vtemp |
|---|
| 1051 | |
|---|
| 1052 | ! initialization module variables |
|---|
| 1053 | subroutine ini_rocketduststorm_mod(ngrid) |
|---|
| 1054 | |
|---|
| 1055 | implicit none |
|---|
| 1056 | |
|---|
| 1057 | integer, intent(in) :: ngrid |
|---|
| 1058 | |
|---|
| 1059 | allocate(dustliftday(ngrid)) |
|---|
| 1060 | allocate(alpha_hmons(ngrid)) |
|---|
| 1061 | |
|---|
| 1062 | end subroutine ini_rocketduststorm_mod |
|---|
| 1063 | |
|---|
| 1064 | subroutine end_rocketduststorm_mod |
|---|
| 1065 | |
|---|
| 1066 | implicit none |
|---|
| 1067 | |
|---|
| 1068 | if (allocated(dustliftday)) deallocate(dustliftday) |
|---|
| 1069 | if (allocated(alpha_hmons)) deallocate(alpha_hmons) |
|---|
| 1070 | |
|---|
| 1071 | end subroutine end_rocketduststorm_mod |
|---|
| 1072 | |
|---|
| 1073 | END MODULE rocketduststorm_mod |
|---|
