| 1 | MODULE topmons_mod |
|---|
| 2 | |
|---|
| 3 | IMPLICIT NONE |
|---|
| 4 | |
|---|
| 5 | CONTAINS |
|---|
| 6 | |
|---|
| 7 | !======================================================================= |
|---|
| 8 | ! ENTRAINMENT OF DUST ABOVE THE SUB-GRID SCALE TOPOGRAPHY |
|---|
| 9 | !======================================================================= |
|---|
| 10 | ! calculation of the vertical wind velocity of topdust tracers |
|---|
| 11 | ! transport scheme of topdust tracers |
|---|
| 12 | ! detrainement of topdust into background dust |
|---|
| 13 | ! ----------------------------------------------------------------------- |
|---|
| 14 | ! Authors: M. Vals; C. Wang; T. Bertrand; F. Forget |
|---|
| 15 | ! Institution: Laboratoire de Meteorologie Dynamique (LMD) Paris, France |
|---|
| 16 | ! ----------------------------------------------------------------------- |
|---|
| 17 | |
|---|
| 18 | SUBROUTINE topmons(ngrid,nlayer,nq,ptime,ptimestep, & |
|---|
| 19 | pq,pdq,pt,pdt,pplev,pplay,pzlev, & |
|---|
| 20 | pzlay,pdtsw,pdtlw, & |
|---|
| 21 | ! input for radiative transfer |
|---|
| 22 | icount,zday,zls,tsurf,igout,aerosol, & |
|---|
| 23 | tauscaling,dust_rad_adjust, & |
|---|
| 24 | IRtoVIScoef, & |
|---|
| 25 | ! input sub-grid scale rocket dust storm |
|---|
| 26 | totstormfract,clearatm, & |
|---|
| 27 | ! input sub-grid scale cloud |
|---|
| 28 | clearsky,totcloudfrac, & |
|---|
| 29 | ! input sub-grid scale mountain |
|---|
| 30 | nohmons, & |
|---|
| 31 | ! output |
|---|
| 32 | pdqtop,wfin,dsodust,dsords,dsotop, & |
|---|
| 33 | tau_pref_scenario,tau_pref_gcm) |
|---|
| 34 | |
|---|
| 35 | USE tracer_mod, only: igcm_topdust_mass,igcm_topdust_number & |
|---|
| 36 | ,igcm_dust_mass,igcm_dust_number & |
|---|
| 37 | ,rho_dust |
|---|
| 38 | USE comcstfi_h, only: r,g,cpp,rcp |
|---|
| 39 | USE dimradmars_mod, only: albedo,naerkind |
|---|
| 40 | USE comsaison_h, only: dist_sol,mu0,fract |
|---|
| 41 | USE surfdat_h, only: emis,co2ice,hmons,summit,alpha_hmons, & |
|---|
| 42 | hsummit,contains_mons |
|---|
| 43 | USE callradite_mod, only: callradite |
|---|
| 44 | |
|---|
| 45 | IMPLICIT NONE |
|---|
| 46 | |
|---|
| 47 | !-------------------------------------------------------- |
|---|
| 48 | ! Input Variables |
|---|
| 49 | !-------------------------------------------------------- |
|---|
| 50 | |
|---|
| 51 | INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points |
|---|
| 52 | INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points |
|---|
| 53 | INTEGER, INTENT(IN) :: nq ! number of tracer species |
|---|
| 54 | REAL, INTENT(IN) :: ptime |
|---|
| 55 | REAL, INTENT(IN) :: ptimestep |
|---|
| 56 | |
|---|
| 57 | REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) ! advected field nq |
|---|
| 58 | REAL, INTENT(IN) :: pdq(ngrid,nlayer,nq)! tendancy field pq |
|---|
| 59 | REAL, INTENT(IN) :: pt(ngrid,nlayer) ! temperature at mid-layer (K) |
|---|
| 60 | REAL, INTENT(IN) :: pdt(ngrid,nlayer) ! tendancy temperature at mid-layer (K) |
|---|
| 61 | |
|---|
| 62 | REAL, INTENT(IN) :: pplay(ngrid,nlayer) ! pressure at middle of the layers |
|---|
| 63 | REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) ! pressure at intermediate levels |
|---|
| 64 | REAL, INTENT(IN) :: pzlay(ngrid,nlayer) ! altitude at the middle of the layers |
|---|
| 65 | REAL, INTENT(IN) :: pzlev(ngrid,nlayer+1) ! altitude at layer boundaries |
|---|
| 66 | |
|---|
| 67 | REAL, INTENT(IN) :: pdtsw(ngrid,nlayer) ! (K/s) env |
|---|
| 68 | REAL, INTENT(IN) :: pdtlw(ngrid,nlayer) ! (K/s) env |
|---|
| 69 | |
|---|
| 70 | ! input for second radiative transfer |
|---|
| 71 | INTEGER, INTENT(INOUT) :: icount |
|---|
| 72 | REAL, INTENT(IN) :: zday |
|---|
| 73 | REAL, INTENT(IN) :: zls |
|---|
| 74 | REAL, INTENT(IN) :: tsurf(ngrid) |
|---|
| 75 | INTEGER, INTENT(IN) :: igout |
|---|
| 76 | REAL, INTENT(INOUT) :: tauscaling(ngrid) |
|---|
| 77 | REAL,INTENT(INOUT) :: dust_rad_adjust(ngrid) |
|---|
| 78 | REAL,INTENT(INOUT) :: IRtoVIScoef(ngrid) ! NB: not modified by this call to callradite, |
|---|
| 79 | ! the OUT is just here because callradite needs it |
|---|
| 80 | |
|---|
| 81 | ! input sub-grid scale rocket dust storm |
|---|
| 82 | LOGICAL, INTENT(IN) :: clearatm |
|---|
| 83 | REAL, INTENT(IN) :: totstormfract(ngrid) |
|---|
| 84 | ! input sub-grid scale water ice clouds |
|---|
| 85 | LOGICAL, INTENT(IN) :: clearsky |
|---|
| 86 | REAL, INTENT(IN) :: totcloudfrac(ngrid) |
|---|
| 87 | ! input sub-grid scale mountain |
|---|
| 88 | LOGICAL, INTENT(IN) :: nohmons |
|---|
| 89 | |
|---|
| 90 | !-------------------------------------------------------- |
|---|
| 91 | ! Output Variables |
|---|
| 92 | !-------------------------------------------------------- |
|---|
| 93 | |
|---|
| 94 | REAL, INTENT(OUT) :: pdqtop(ngrid,nlayer,nq) ! tendancy field for dust |
|---|
| 95 | REAL, INTENT(OUT) :: wfin(ngrid,nlayer+1) ! final vertical wind velocity: combination of both wup and wrad |
|---|
| 96 | ! wfin < 0 means upward wind (in pressure levels/s) |
|---|
| 97 | |
|---|
| 98 | ! output for second radiative transfer |
|---|
| 99 | REAL, INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) |
|---|
| 100 | REAL, INTENT(OUT) :: dsodust(ngrid,nlayer) |
|---|
| 101 | REAL, INTENT(OUT) :: dsords(ngrid,nlayer) |
|---|
| 102 | REAL, INTENT(OUT) :: dsotop(ngrid,nlayer) |
|---|
| 103 | REAL,INTENT(OUT) :: tau_pref_scenario(ngrid) ! prescribed dust column |
|---|
| 104 | ! visible opacity at odpref |
|---|
| 105 | REAL,INTENT(OUT) :: tau_pref_gcm(ngrid) ! dust column visible opacity at |
|---|
| 106 | ! odpref in the GCM |
|---|
| 107 | |
|---|
| 108 | !-------------------------------------------------------- |
|---|
| 109 | ! Local variables |
|---|
| 110 | !-------------------------------------------------------- |
|---|
| 111 | LOGICAL,SAVE :: firstcall=.true. |
|---|
| 112 | |
|---|
| 113 | !$OMP THREADPRIVATE(firstcall) |
|---|
| 114 | |
|---|
| 115 | INTEGER l,ig,tsub,iq,ll |
|---|
| 116 | REAL zq0(ngrid,nlayer,nq) ! initial tracers |
|---|
| 117 | REAL zq(ngrid,nlayer,nq) ! updated tracers |
|---|
| 118 | REAL,PARAMETER:: mu0lim=0.001 ! solar zenith angle limit to determine the daytime |
|---|
| 119 | |
|---|
| 120 | ! Variables for radiative transfer |
|---|
| 121 | REAL fluxsurf_lw1(ngrid) |
|---|
| 122 | REAL fluxsurf_dn_sw1(ngrid,2),fluxsurf_up_sw1(ngrid,2) |
|---|
| 123 | REAL fluxtop_lw1(ngrid) |
|---|
| 124 | REAL fluxtop_dn_sw1(ngrid,2),fluxtop_up_sw1(ngrid,2) |
|---|
| 125 | REAL tau(ngrid,naerkind) |
|---|
| 126 | REAL taucloudtes(ngrid) |
|---|
| 127 | REAL rdust(ngrid,nlayer) |
|---|
| 128 | REAL rstormdust(ngrid,nlayer) |
|---|
| 129 | REAL rtopdust(ngrid,nlayer) |
|---|
| 130 | REAL rice(ngrid,nlayer) |
|---|
| 131 | REAL nuice(ngrid,nlayer) |
|---|
| 132 | DOUBLE PRECISION riceco2(ngrid,nlayer) |
|---|
| 133 | REAL nuiceco2(ngrid,nlayer) |
|---|
| 134 | |
|---|
| 135 | ! Temperature profile |
|---|
| 136 | REAL zt(ngrid,nlayer) ! actual temperature at mid-layer (K) |
|---|
| 137 | REAL ztlev(ngrid,nlayer) ! temperature at intermediate levels l+1/2 without last level |
|---|
| 138 | REAL zdtlw1(ngrid,nlayer) ! (K/s) topdust |
|---|
| 139 | REAL zdtsw1(ngrid,nlayer) ! (K/s) topdust |
|---|
| 140 | REAL zdtlw1_lev(ngrid,nlayer),zdtsw1_lev(ngrid,nlayer) ! rad. heating rate at intermediate levels l+1/2 |
|---|
| 141 | REAL zdtlw_lev(ngrid,nlayer),zdtsw_lev(ngrid,nlayer) ! rad. heating rate at intermediate levels l+1/2 |
|---|
| 142 | REAL zdtvert(ngrid,nlayer) ! dT/dz , lapse rate |
|---|
| 143 | REAL deltahr(ngrid,nlayer+1) ! extra heating between the concentrated and non-concentrated dust |
|---|
| 144 | |
|---|
| 145 | REAL newzt(ngrid,nlayer) ! recalculated temperature with extra radiative heating |
|---|
| 146 | REAL t_top(ngrid,nlayer) ! estimated temperature above the mountain |
|---|
| 147 | REAL t_top_lev(ngrid,nlayer)! estimated temperature above the mountain at the levels |
|---|
| 148 | REAL dt_top(ngrid) ! temperature difference between the summit and the environment |
|---|
| 149 | INTEGER lmons(ngrid) ! layer of the mountain's slope |
|---|
| 150 | INTEGER lsummit(ngrid) ! layer of the mountain's summit |
|---|
| 151 | REAL t_env(ngrid,nlayer) ! recalculated temperature of the environment air "next" to the mountain |
|---|
| 152 | |
|---|
| 153 | ! Vertical velocity profile |
|---|
| 154 | REAL wrad(ngrid,nlayer+1) ! vertical wind velocity resulting from radiative heating |
|---|
| 155 | REAL wup(ngrid,nlayer+1) ! vertical wind velocity calculated above the sub-grid scale mountain |
|---|
| 156 | REAL wup2(ngrid,nlayer+1) ! square of the vertical velocity calculated above the sub-grid scale mountain |
|---|
| 157 | REAL w0(ngrid) ! vertical wind velocity convergence at the top of the sub-grid scale mountain |
|---|
| 158 | INTEGER lwmax(ngrid) ! level of the maximum vertical velocity |
|---|
| 159 | |
|---|
| 160 | REAL,PARAMETER:: wmax=10. ! maximum vertical velocity resulting from radiative heating (m/s) |
|---|
| 161 | REAL,PARAMETER:: secu=3. ! coefficient on wfin to avoid dust crossing many layers during subtimestep |
|---|
| 162 | REAL,PARAMETER:: k0=0.25 !max 10m/s: k0=1. !max 3m/s: k0=0.25 |
|---|
| 163 | REAL,PARAMETER:: k1=5.e-4 !max 10m/s: k1=5.e-4 !max 3m/s: k1=5.e-4 |
|---|
| 164 | REAL,PARAMETER:: k2=5.e-3 !max 10m/s: k2=30.e-3 !max 3m/s: k2=5.e-3 |
|---|
| 165 | |
|---|
| 166 | ! Entrainment |
|---|
| 167 | REAL entr(ngrid) ! entrainment flux |
|---|
| 168 | REAL rho(ngrid,nlayer),rhobarz(ngrid,nlayer+1) ! density, density at the levels |
|---|
| 169 | REAL masse(ngrid,nlayer) ! air mass |
|---|
| 170 | REAL masse_pbl(ngrid) ! total air mass within the PBL |
|---|
| 171 | REAL masse_pbl_dust_mass(ngrid),masse_pbl_dust_number(ngrid) ! total air mass + dust mass/number within the PBL (kg/m2) |
|---|
| 172 | REAL dqm_mass_pbl(ngrid),dqm_number_pbl(ngrid) ! total mass of the entrained dust mass/number from the PBL (kg/m2) |
|---|
| 173 | |
|---|
| 174 | ! Van Leer |
|---|
| 175 | REAL zq_dust_mass(ngrid,nlayer),zq_dust_number(ngrid,nlayer) ! mixing ratio background dust (kg/kg) |
|---|
| 176 | REAL zq_topdust_mass(ngrid,nlayer),zq_topdust_number(ngrid,nlayer) ! mixing ratio topdust (kg/kg) |
|---|
| 177 | REAL mr_dust_mass(ngrid,nlayer),mr_dust_number(ngrid,nlayer) ! mixing ratio background dust (kg/kg) |
|---|
| 178 | REAL mr_topdust_mass(ngrid,nlayer),mr_topdust_number(ngrid,nlayer) ! mixing ratio background dust + concentrated topdust (kg/kg) |
|---|
| 179 | |
|---|
| 180 | REAL w(ngrid,nlayer) ! vertical flux above the sub-grid scale mountain |
|---|
| 181 | REAL wqmass(ngrid,nlayer+1) ! flux mass |
|---|
| 182 | REAL wqnumber(ngrid,nlayer+1) ! flux number |
|---|
| 183 | REAL qbar_mass_pbl(ngrid) ! total dust mass mixing ratio within the PBL for Van Leer |
|---|
| 184 | REAL qbar_number_pbl(ngrid) ! total dust number mixing ratio within the PBL for Van Leer |
|---|
| 185 | REAL masse_col(nlayer) ! masse of the atmosphere in one column |
|---|
| 186 | |
|---|
| 187 | REAL dqvl_dust_mass(ngrid,nlayer) ! tendancy pdq dust mass after vertical transport |
|---|
| 188 | REAL dqvl_dust_number(ngrid,nlayer) ! tendancy pdq dust number after vertical transport |
|---|
| 189 | REAL dqvl_topdust_mass(ngrid,nlayer) ! tendancy pdq topdust mass after vertical transport |
|---|
| 190 | REAL dqvl_topdust_number(ngrid,nlayer) ! tendancy pdq topdust number after vertical transport |
|---|
| 191 | |
|---|
| 192 | INTEGER nsubtimestep(ngrid) ! number of subtimestep when calling Van Leer scheme |
|---|
| 193 | REAL subtimestep(ngrid) ! ptimestep/nsubtimestep |
|---|
| 194 | REAL dtmax ! considered minimum time interval needed for the dust to cross one vertical layer |
|---|
| 195 | |
|---|
| 196 | ! Detrainment |
|---|
| 197 | REAL coefdetrain(ngrid,nlayer) ! coefficient for detrainment : % of topdust detrained |
|---|
| 198 | REAL dqdet_topdust_mass(ngrid,nlayer) ! tendancy pdq topdust mass after detrainment only |
|---|
| 199 | REAL dqdet_topdust_number(ngrid,nlayer) ! tendancy pdq topdust number after detrainment only |
|---|
| 200 | |
|---|
| 201 | ! Diagnostics |
|---|
| 202 | REAL zlaywmax(ngrid) |
|---|
| 203 | REAL detrain_rate(ngrid,nlayer) |
|---|
| 204 | |
|---|
| 205 | ! ********************************************************************** |
|---|
| 206 | ! ********************************************************************** |
|---|
| 207 | ! Parametrization of the entrainment of dust by slope winds above the |
|---|
| 208 | ! converging sub-grid scale topography ("mountain top dust flows") |
|---|
| 209 | ! ********************************************************************** |
|---|
| 210 | ! ********************************************************************** |
|---|
| 211 | !! 1. Radiative transfer in topdust |
|---|
| 212 | !! 2. Compute vertical velocity for topdust |
|---|
| 213 | !! 3. Vertical transport |
|---|
| 214 | !! 4. Detrainment |
|---|
| 215 | ! ********************************************************************** |
|---|
| 216 | ! ********************************************************************** |
|---|
| 217 | |
|---|
| 218 | ! ********************************************************************** |
|---|
| 219 | ! 0. Initializations |
|---|
| 220 | ! ********************************************************************** |
|---|
| 221 | aerosol(:,:,:)=0. |
|---|
| 222 | pdqtop(:,:,:) = 0. |
|---|
| 223 | dqvl_topdust_mass(:,:)=0. |
|---|
| 224 | dqvl_topdust_number(:,:)=0. |
|---|
| 225 | dqvl_dust_mass(:,:)=0. |
|---|
| 226 | dqvl_dust_number(:,:)=0. |
|---|
| 227 | dqdet_topdust_mass(:,:)=0. |
|---|
| 228 | dqdet_topdust_number(:,:)=0. |
|---|
| 229 | wup(:,:)=0. |
|---|
| 230 | wup2(:,:)=0. |
|---|
| 231 | wrad(:,:)=0. |
|---|
| 232 | w0(:)=0. |
|---|
| 233 | wfin(:,:)=0. |
|---|
| 234 | w(:,:)=0. |
|---|
| 235 | zq(:,:,:) = 0. |
|---|
| 236 | zq0(:,:,:) = 0. |
|---|
| 237 | entr(:) = 0. |
|---|
| 238 | mr_dust_mass(:,:) = 0. |
|---|
| 239 | mr_dust_number(:,:) = 0. |
|---|
| 240 | mr_topdust_mass(:,:) = 0. |
|---|
| 241 | mr_topdust_number(:,:) = 0. |
|---|
| 242 | t_top(:,:) = 0. |
|---|
| 243 | dt_top(:) = 0. |
|---|
| 244 | newzt(:,:) = 0. |
|---|
| 245 | lmons(:) = 1 |
|---|
| 246 | lsummit(:) =1 |
|---|
| 247 | lwmax(:) = 1 |
|---|
| 248 | deltahr(:,:) = 0. |
|---|
| 249 | zdtvert(:,:) = 0. |
|---|
| 250 | wqmass(:,:) = 0. |
|---|
| 251 | wqnumber(:,:) = 0. |
|---|
| 252 | coefdetrain(:,:) = 1. |
|---|
| 253 | dqm_mass_pbl(:)=0. |
|---|
| 254 | dqm_number_pbl(:)=0. |
|---|
| 255 | nsubtimestep(:)=1 |
|---|
| 256 | masse_pbl(:)=0. |
|---|
| 257 | detrain_rate(:,:) = 0. |
|---|
| 258 | |
|---|
| 259 | !! Update the initial temperature |
|---|
| 260 | zt(1:ngrid,1:nlayer)=pt(1:ngrid,1:nlayer)+pdt(1:ngrid,1:nlayer)*ptimestep |
|---|
| 261 | newzt(1:ngrid,1:nlayer)=zt(1:ngrid,1:nlayer) |
|---|
| 262 | t_env(1:ngrid,1:nlayer)=zt(1:ngrid,1:nlayer) |
|---|
| 263 | t_top(1:ngrid,1:nlayer)=zt(1:ngrid,1:nlayer) |
|---|
| 264 | |
|---|
| 265 | !! Update the initial mixing ratios |
|---|
| 266 | zq0(1:ngrid,1:nlayer,1:nq)=pq(1:ngrid,1:nlayer,1:nq)+pdq(1:ngrid,1:nlayer,1:nq)*ptimestep ! update of the background dust after rocket dust storm scheme |
|---|
| 267 | zq(1:ngrid,1:nlayer,1:nq)=zq0(1:ngrid,1:nlayer,1:nq) |
|---|
| 268 | zq_dust_mass(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_dust_mass) |
|---|
| 269 | zq_dust_number(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_dust_number) |
|---|
| 270 | zq_topdust_mass(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_topdust_mass) |
|---|
| 271 | zq_topdust_number(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_topdust_number) |
|---|
| 272 | |
|---|
| 273 | ! ---------------------------------------------------------------------- |
|---|
| 274 | ! 1. Radiative heating |
|---|
| 275 | ! ---------------------------------------------------------------------- |
|---|
| 276 | |
|---|
| 277 | ! ********************************************************************* |
|---|
| 278 | ! 1.1. Call the second radiative transfer for topdust, obtain the extra heating |
|---|
| 279 | ! ********************************************************************* |
|---|
| 280 | |
|---|
| 281 | ! NB: since the only grid meshes that matter for the radiative transfer |
|---|
| 282 | ! are the ones that contain a mount (contains_mons(ig)=.true.), |
|---|
| 283 | ! it could be relevant to optimize the code by calling the RT |
|---|
| 284 | ! only for those grid meshes. |
|---|
| 285 | ! See Ticket #92 on trac.lmd.jussieu.fr/Planeto |
|---|
| 286 | CALL callradite(icount,ngrid,nlayer,nq,zday,zls,zq,albedo, & |
|---|
| 287 | emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, & |
|---|
| 288 | zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_dn_sw1,fluxsurf_up_sw1, & |
|---|
| 289 | fluxtop_lw1,fluxtop_dn_sw1,fluxtop_up_sw1, & |
|---|
| 290 | tau_pref_scenario,tau_pref_gcm, & |
|---|
| 291 | tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef, & |
|---|
| 292 | taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice,rstormdust,rtopdust, & |
|---|
| 293 | totstormfract,clearatm,dsords,dsotop,nohmons,& |
|---|
| 294 | clearsky,totcloudfrac) |
|---|
| 295 | ! ********************************************************************** |
|---|
| 296 | ! 1.2. Compute radiative vertical velocity for entrained topdust |
|---|
| 297 | ! ********************************************************************** |
|---|
| 298 | DO ig=1,ngrid |
|---|
| 299 | IF ( (mu0(ig) .gt. mu0lim) .and. (contains_mons(ig)) ) THEN |
|---|
| 300 | !! mu0(ig)>mu0lim ensures that it is daytime |
|---|
| 301 | !! contains_mons=True ensures that there is a mount in the mesh and alpha_hmons>0 |
|---|
| 302 | |
|---|
| 303 | !! ********************************************************************** |
|---|
| 304 | !! Temperature profile above the mountain and in the close environment |
|---|
| 305 | call t_topmons(nlayer,summit(ig),hsummit(ig),hmons(ig),zt(ig,:),pzlay(ig,:), & |
|---|
| 306 | t_top(ig,:),dt_top(ig),t_env(ig,:),lsummit(ig),lmons(ig)) |
|---|
| 307 | !! ********************************************************************** |
|---|
| 308 | !! Calculation of the extra heating : computing heating rates |
|---|
| 309 | !! gradient at boundaries of each layer |
|---|
| 310 | zdtlw1_lev(ig,1)=0. |
|---|
| 311 | zdtsw1_lev(ig,1)=0. |
|---|
| 312 | zdtlw_lev(ig,1)=0. |
|---|
| 313 | zdtsw_lev(ig,1)=0. |
|---|
| 314 | ztlev(ig,1)=zt(ig,1) |
|---|
| 315 | t_top_lev(ig,1)=tsurf(ig) |
|---|
| 316 | |
|---|
| 317 | DO l=1,nlayer-1 |
|---|
| 318 | !! Calculation for the background dust AND the concentrated topdust |
|---|
| 319 | zdtlw1_lev(ig,l+1)=(zdtlw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ & |
|---|
| 320 | zdtlw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / & |
|---|
| 321 | (pzlay(ig,l+1)-pzlay(ig,l)) |
|---|
| 322 | |
|---|
| 323 | zdtsw1_lev(ig,l+1)=(zdtsw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ & |
|---|
| 324 | zdtsw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / & |
|---|
| 325 | (pzlay(ig,l+1)-pzlay(ig,l)) |
|---|
| 326 | !! Calculation for the background dust only |
|---|
| 327 | zdtlw_lev(ig,l+1)=(pdtlw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ & |
|---|
| 328 | pdtlw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / & |
|---|
| 329 | (pzlay(ig,l+1)-pzlay(ig,l)) |
|---|
| 330 | |
|---|
| 331 | zdtsw_lev(ig,l+1)=(pdtsw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ & |
|---|
| 332 | pdtsw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / & |
|---|
| 333 | (pzlay(ig,l+1)-pzlay(ig,l)) |
|---|
| 334 | !! Temperature at the levels |
|---|
| 335 | ztlev(ig,l+1)=(zt(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ & |
|---|
| 336 | zt(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / & |
|---|
| 337 | (pzlay(ig,l+1)-pzlay(ig,l)) |
|---|
| 338 | !! Temperature above the mountain at the levels |
|---|
| 339 | t_top_lev(ig,l+1)=(t_top(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ & |
|---|
| 340 | t_top(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / & |
|---|
| 341 | (pzlay(ig,l+1)-pzlay(ig,l)) |
|---|
| 342 | ENDDO ! DO l=1,nlayer-1 |
|---|
| 343 | !! ********************************************************************** |
|---|
| 344 | !! Vertical velocity profile for extra heating balanced by adiabatic cooling |
|---|
| 345 | !! Gradient at boundaries of each layer, start from surface |
|---|
| 346 | zdtvert(ig,1)=0. |
|---|
| 347 | DO l=1,nlayer-1 |
|---|
| 348 | zdtvert(ig,l+1)=(t_top_lev(ig,l+1)-t_top_lev(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l)) |
|---|
| 349 | ENDDO |
|---|
| 350 | !! Extra heating balanced by adiabatic cooling |
|---|
| 351 | DO l=1,nlayer |
|---|
| 352 | deltahr(ig,l)=(zdtlw1_lev(ig,l)+zdtsw1_lev(ig,l)) & |
|---|
| 353 | -(zdtlw_lev(ig,l)+zdtsw_lev(ig,l)) |
|---|
| 354 | wrad(ig,l)=-deltahr(ig,l)/(g/cpp+ & |
|---|
| 355 | max(zdtvert(ig,l),-0.99*g/cpp)) |
|---|
| 356 | !! Limit of the vertical wind velocity in case of lapse rate close to adiabatic |
|---|
| 357 | wrad(ig,l)=max(wrad(ig,l),-wmax) |
|---|
| 358 | wrad(ig,l)=min(wrad(ig,l),wmax) |
|---|
| 359 | !! ATTENTION !! to simplify here, no downward velocity calculated |
|---|
| 360 | if (wrad(ig,l).gt.0.) then |
|---|
| 361 | wrad(ig,l)=0. |
|---|
| 362 | endif |
|---|
| 363 | ENDDO |
|---|
| 364 | ENDIF ! IF ((mu0(ig) .gt. mu0lim) .and. contains_mons(ig)) |
|---|
| 365 | ENDDO ! DO ig=1,ngrid |
|---|
| 366 | |
|---|
| 367 | ! ---------------------------------------------------------------------- |
|---|
| 368 | ! 2. Vertical velocity profile |
|---|
| 369 | ! ---------------------------------------------------------------------- |
|---|
| 370 | |
|---|
| 371 | ! ********************************************************************** |
|---|
| 372 | ! 2.1 Compute total vertical velocity for entrained topdust: buoyancy + radiative |
|---|
| 373 | ! ********************************************************************** |
|---|
| 374 | DO ig=1,ngrid |
|---|
| 375 | IF ( (mu0(ig) .gt. mu0lim) .and. (contains_mons(ig)) ) THEN |
|---|
| 376 | !! Positive buoyancy: negative vertical velocity entrains UP |
|---|
| 377 | IF (dt_top(ig) .gt. 0.) THEN |
|---|
| 378 | !! Convergence of vertical winds at the top of the mountain |
|---|
| 379 | w0(ig)=-k0*g*sqrt((dt_top(ig)/t_env(ig,lsummit(ig)))) |
|---|
| 380 | !! Case the vertical velocity at the top of the mountain > wrad |
|---|
| 381 | IF (w0(ig).lt.wrad(ig,lsummit(ig)+1)) then |
|---|
| 382 | wup(ig,lsummit(ig)+1)=w0(ig) |
|---|
| 383 | wfin(ig,lsummit(ig)+1)=w0(ig) |
|---|
| 384 | newzt(ig,lsummit(ig))= t_top(ig,lsummit(ig)) |
|---|
| 385 | !! Temperature profile from the top of the summit to the top of the atmosphere |
|---|
| 386 | DO l=lsummit(ig)+1,nlayer-1 |
|---|
| 387 | !ATTENTION: test without wrad, T is only decreasing because of the adiabatic gradient |
|---|
| 388 | !newzt(ig,l)=newzt(ig,l-1)-(g/cpp)* & |
|---|
| 389 | ! (pzlay(ig,l)-pzlay(ig,l-1)) |
|---|
| 390 | !! Case of superadiabatic layer |
|---|
| 391 | IF (zdtvert(ig,l) .lt. -g/cpp) then |
|---|
| 392 | newzt(ig,l)=t_top(ig,l) |
|---|
| 393 | !! New temperature due to radiative heating balanced by adiabatic cooling |
|---|
| 394 | ELSE |
|---|
| 395 | newzt(ig,l)=newzt(ig,l-1)+ & |
|---|
| 396 | ( (deltahr(ig,l)/(-wfin(ig,l)))-g/cpp )* & |
|---|
| 397 | ( pzlay(ig,l)-pzlay(ig,l-1) ) |
|---|
| 398 | ENDIF ! (zdtvert(ig,l) .lt. -g/cpp) |
|---|
| 399 | !! Vertical velocity profile due to the presence of the mountain with the new temperature |
|---|
| 400 | !! Square of the vertical velocity profile with the new temperature |
|---|
| 401 | wup2(ig,l+1)=(1.-2*k1*(pzlev(ig,l+1)-pzlev(ig,l)))*& |
|---|
| 402 | (wup(ig,l)**2)+2.*k2*g*(pzlev(ig,l+1) & |
|---|
| 403 | -pzlev(ig,l))*((newzt(ig,l) & |
|---|
| 404 | -t_env(ig,l))/t_env(ig,l)) |
|---|
| 405 | !! Vertical velocity profile if W2 > 0 |
|---|
| 406 | IF (wup2(ig,l+1) .gt. 0.) THEN |
|---|
| 407 | wup(ig,l+1)=-sqrt(wup2(ig,l+1)) |
|---|
| 408 | wfin(ig,l+1)=wup(ig,l+1) |
|---|
| 409 | IF (wup(ig,l+1) .gt. wrad(ig,l+1)) then |
|---|
| 410 | DO ll=l,nlayer-1 |
|---|
| 411 | newzt(ig,ll)=zt(ig,ll) |
|---|
| 412 | wfin(ig,ll+1)=wrad(ig,ll+1) |
|---|
| 413 | ENDDO |
|---|
| 414 | EXIT |
|---|
| 415 | ENDIF |
|---|
| 416 | !! Vertical velocity profile if W2 < 0 |
|---|
| 417 | ELSE |
|---|
| 418 | DO ll=l,nlayer-1 |
|---|
| 419 | newzt(ig,ll)=zt(ig,ll) |
|---|
| 420 | wfin(ig,ll+1)=wrad(ig,ll+1) |
|---|
| 421 | ENDDO |
|---|
| 422 | EXIT |
|---|
| 423 | ENDIF ! IF (wup2(ig,l+1) .gt. 0.) |
|---|
| 424 | ENDDO ! DO l=lsummit(ig)+1,nlayer-1 |
|---|
| 425 | !! Case the vertical velocity at the top of the mountain < wrad |
|---|
| 426 | ELSE |
|---|
| 427 | DO ll=lsummit(ig),nlayer-1 |
|---|
| 428 | newzt(ig,ll)=zt(ig,ll) |
|---|
| 429 | wfin(ig,ll+1)=wrad(ig,ll+1) |
|---|
| 430 | ENDDO |
|---|
| 431 | ENDIF !(w0(ig).lt.wrad(ig,lsummit(ig)+1)) |
|---|
| 432 | !! Positive buoyancy: positive vertical velocity entrains DOWN |
|---|
| 433 | ELSE |
|---|
| 434 | DO l=lsummit(ig)+1,nlayer |
|---|
| 435 | wfin(ig,l)=wrad(ig,l) |
|---|
| 436 | ENDDO |
|---|
| 437 | ENDIF ! IF (dt_top(ig) .gt. 0.) |
|---|
| 438 | ! ********************************************************************** |
|---|
| 439 | ! 2.2 Find the level lwmax of the maximum vertical velocity caused by the |
|---|
| 440 | ! mountain presence (wup), where to entrain the dust from the PBL |
|---|
| 441 | ! ********************************************************************** |
|---|
| 442 | IF (minloc(wup(ig,:),dim=1).lt.nlayer) THEN |
|---|
| 443 | lwmax(ig)=minloc(wup(ig,:),dim=1) |
|---|
| 444 | ENDIF |
|---|
| 445 | zlaywmax(ig)=pzlay(ig,lwmax(ig)) ! wup(lwmax) acts on level lwmax, going to layer lwmax |
|---|
| 446 | ! ********************************************************************** |
|---|
| 447 | ! 2.3 Compute the subtimestep to conserve the mass in the Van Leer transport |
|---|
| 448 | ! ********************************************************************** |
|---|
| 449 | !! Calculation of the subtimesteps |
|---|
| 450 | dtmax=ptimestep |
|---|
| 451 | DO l=lwmax(ig),nlayer |
|---|
| 452 | IF (wfin(ig,l).lt.0.) THEN |
|---|
| 453 | dtmax=min(dtmax,(pzlev(ig,l)-pzlev(ig,l-1))/ & |
|---|
| 454 | (secu*abs(wfin(ig,l)))) |
|---|
| 455 | ELSE IF (wfin(ig,l).gt.0.) then |
|---|
| 456 | dtmax=min(dtmax,(pzlev(ig,l+1)-pzlev(ig,l))/ & |
|---|
| 457 | (secu*abs(wfin(ig,l)))) |
|---|
| 458 | ENDIF |
|---|
| 459 | ENDDO |
|---|
| 460 | nsubtimestep(ig)= int(ptimestep/dtmax) |
|---|
| 461 | subtimestep(ig)=ptimestep/float(nsubtimestep(ig)) |
|---|
| 462 | !! Mass flux generated by wup in kg/m2 |
|---|
| 463 | DO l=1,nlayer!+1 |
|---|
| 464 | w(ig,l)=wfin(ig,l)*pplev(ig,l)/(r*ztlev(ig,l)) & |
|---|
| 465 | *subtimestep(ig) |
|---|
| 466 | ENDDO ! l=1,nlayer |
|---|
| 467 | |
|---|
| 468 | ENDIF ! ((mu0(ig) .gt. mu0lim)... |
|---|
| 469 | ENDDO ! DO ig=1,ngrid |
|---|
| 470 | |
|---|
| 471 | ! ---------------------------------------------------------------------- |
|---|
| 472 | ! 3. Dust entrainment by slope winds above the mountain |
|---|
| 473 | ! ---------------------------------------------------------------------- |
|---|
| 474 | !! rho on the layer |
|---|
| 475 | rho(:,:)=pplay(:,:)/(r*pt(:,:)) |
|---|
| 476 | !! rhobarz(l) = rho(l+1/2): rho on the levels |
|---|
| 477 | rhobarz(:,1)=rho(:,1) |
|---|
| 478 | DO l=2,nlayer |
|---|
| 479 | rhobarz(:,l)=pplev(:,l)/(r*0.5*(pt(:,l)+pt(:,l-1))) |
|---|
| 480 | ENDDO |
|---|
| 481 | rhobarz(:,nlayer+1)=0. !top of the model |
|---|
| 482 | !! Mass computation |
|---|
| 483 | DO l=1,nlayer |
|---|
| 484 | masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g |
|---|
| 485 | ENDDO |
|---|
| 486 | |
|---|
| 487 | DO ig=1,ngrid |
|---|
| 488 | IF ( (mu0(ig) .gt. mu0lim) .and. (contains_mons(ig)) ) THEN |
|---|
| 489 | !! Total air mass within the PBL before entrainment (=> by PBL we mean between the surface and the layer where the vertical wind is maximum) |
|---|
| 490 | masse_pbl(ig)=0. |
|---|
| 491 | DO l=1,lwmax(ig)-1 |
|---|
| 492 | masse_pbl(ig)=masse_pbl(ig)+(pplev(ig,l)-pplev(ig,l+1))/g |
|---|
| 493 | ENDDO |
|---|
| 494 | ENDIF ! ((mu0(ig) .gt. mu0lim)... |
|---|
| 495 | ENDDO ! ig=1,ngrid |
|---|
| 496 | ! ********************************************************************** |
|---|
| 497 | ! 3.1. Entrainment |
|---|
| 498 | ! ********************************************************************** |
|---|
| 499 | DO ig=1,ngrid |
|---|
| 500 | IF ( (mu0(ig) .gt. mu0lim) .and. (contains_mons(ig)) .and. (masse_pbl(ig) .gt. 0.) ) THEN |
|---|
| 501 | !! Transport of background dust + concentrated topdust above lwmax |
|---|
| 502 | DO l=lwmax(ig),nlayer |
|---|
| 503 | mr_dust_mass(ig,l) = zq_dust_mass(ig,l) |
|---|
| 504 | mr_dust_number(ig,l) = zq_dust_number(ig,l) |
|---|
| 505 | mr_topdust_mass(ig,l) = zq_dust_mass(ig,l) & |
|---|
| 506 | +zq_topdust_mass(ig,l)/alpha_hmons(ig) |
|---|
| 507 | mr_topdust_number(ig,l) = zq_dust_number(ig,l) & |
|---|
| 508 | +zq_topdust_number(ig,l)/alpha_hmons(ig) |
|---|
| 509 | ENDDO ! l=lwmax(ig),nlayer |
|---|
| 510 | !! Entrainment flux from the PBL |
|---|
| 511 | entr(ig) = alpha_hmons(ig)*wup(ig,lwmax(ig))*rhobarz(ig,lwmax(ig))/masse_pbl(ig) |
|---|
| 512 | !! If entrains more than available masse_pbl (abs(entr) x mass_pbl x ptimestep > masse_pbl) |
|---|
| 513 | if (abs(entr(ig)) .gt. 1./ptimestep) then |
|---|
| 514 | entr(ig) = -(1./ptimestep) |
|---|
| 515 | end if |
|---|
| 516 | ! ********************************************************************** |
|---|
| 517 | ! 3.2. Vertical transport: Van Leer |
|---|
| 518 | ! ********************************************************************** |
|---|
| 519 | !! Loop over the subtimesteps |
|---|
| 520 | DO tsub=1,nsubtimestep(ig) |
|---|
| 521 | !! qbar_pbl: total mixing ratio within the PBL |
|---|
| 522 | masse_pbl_dust_mass(ig)=0. |
|---|
| 523 | masse_pbl_dust_number(ig)=0. |
|---|
| 524 | DO l=1,lwmax(ig)-1 |
|---|
| 525 | masse_pbl_dust_mass(ig)=masse_pbl_dust_mass(ig)+zq_dust_mass(ig,l)*(pplev(ig,l)-pplev(ig,l+1))/g |
|---|
| 526 | masse_pbl_dust_number(ig)=masse_pbl_dust_number(ig)+zq_dust_number(ig,l)*(pplev(ig,l)-pplev(ig,l+1))/g |
|---|
| 527 | ENDDO |
|---|
| 528 | qbar_mass_pbl(ig)=masse_pbl_dust_mass(ig)/masse_pbl(ig) |
|---|
| 529 | qbar_number_pbl(ig)=masse_pbl_dust_mass(ig)/masse_pbl(ig) |
|---|
| 530 | !! integration of the equation dq/dt = -(entr)q, with entr constant --> w0 < 0 when up so the equation is dq/dt = (entr)q |
|---|
| 531 | dqm_mass_pbl(:)=0. |
|---|
| 532 | dqm_number_pbl(:)=0. |
|---|
| 533 | DO l=1,lwmax(ig)-1 ! this time we take the dust from the surface up to the level, where the velocity is maximum |
|---|
| 534 | !! Total dust entrained from the PBL: |
|---|
| 535 | dqm_mass_pbl(ig)=dqm_mass_pbl(ig)+(1.-exp(entr(ig)*subtimestep(ig)))*zq_dust_mass(ig,l)*(pplev(ig,l)-pplev(ig,l+1))/g |
|---|
| 536 | dqm_number_pbl(ig)=dqm_number_pbl(ig)+(1.-exp(entr(ig)*subtimestep(ig)))*zq_dust_number(ig,l)*(pplev(ig,l)-pplev(ig,l+1))/g |
|---|
| 537 | !! Background dust after entrainment (explicit: qdust(ig,l)=qdust(ig,l)-entr(ig)*qdust(ig,l)*ptimestep) |
|---|
| 538 | zq_dust_mass(ig,l)=zq_dust_mass(ig,l)-(1.-exp(entr(ig)*subtimestep(ig)))*zq_dust_mass(ig,l) |
|---|
| 539 | zq_dust_number(ig,l)=zq_dust_number(ig,l)-(1.-exp(entr(ig)*subtimestep(ig)))*zq_dust_number(ig,l) |
|---|
| 540 | ENDDO |
|---|
| 541 | !! Van Leer scheme (from lwmax to nlayer) |
|---|
| 542 | wqmass(ig,:)=0. |
|---|
| 543 | wqnumber(ig,:)=0. |
|---|
| 544 | CALL van_leer(nlayer,mr_topdust_mass(ig,:),2.,lwmax(ig), & |
|---|
| 545 | masse(ig,:),w(ig,:),wqmass(ig,:),masse_pbl(ig),dqm_mass_pbl(ig),alpha_hmons(ig),qbar_mass_pbl(ig)) |
|---|
| 546 | CALL van_leer(nlayer,mr_topdust_number(ig,:),2.,lwmax(ig), & |
|---|
| 547 | masse(ig,:),w(ig,:),wqnumber(ig,:),masse_pbl(ig),dqm_number_pbl(ig),alpha_hmons(ig),qbar_number_pbl(ig)) |
|---|
| 548 | ENDDO !tsub=... |
|---|
| 549 | ! ********************************************************************** |
|---|
| 550 | ! 3.3. Re-calculation of the mixing ratios after vertical transport |
|---|
| 551 | ! ********************************************************************** |
|---|
| 552 | !! Redistribution from lwmax to nlayer |
|---|
| 553 | DO l=lwmax(ig),nlayer |
|---|
| 554 | !! General and "healthy" case |
|---|
| 555 | IF (mr_topdust_mass(ig,l).ge.mr_dust_mass(ig,l)) THEN |
|---|
| 556 | zq_dust_mass(ig,l) = mr_dust_mass(ig,l) |
|---|
| 557 | zq_dust_number(ig,l) = mr_dust_number(ig,l) |
|---|
| 558 | zq_topdust_mass(ig,l) = alpha_hmons(ig)*(mr_topdust_mass(ig,l)-mr_dust_mass(ig,l)) |
|---|
| 559 | zq_topdust_number(ig,l) = alpha_hmons(ig)*(mr_topdust_number(ig,l)-mr_dust_number(ig,l)) |
|---|
| 560 | !! Particular case: there is less than initial dust mixing ratio after the vertical transport |
|---|
| 561 | ELSE |
|---|
| 562 | zq_dust_mass(ig,l) = (1.-alpha_hmons(ig))*mr_dust_mass(ig,l)+alpha_hmons(ig)*mr_topdust_mass(ig,l) |
|---|
| 563 | zq_dust_number(ig,l) = (1.-alpha_hmons(ig))*mr_dust_number(ig,l)+alpha_hmons(ig)*mr_topdust_number(ig,l) |
|---|
| 564 | zq_topdust_mass(ig,l) = 0. |
|---|
| 565 | zq_topdust_number(ig,l) = 0. |
|---|
| 566 | ENDIF |
|---|
| 567 | ENDDO ! DO l=lwmax(ig),nlayer |
|---|
| 568 | ! ********************************************************************** |
|---|
| 569 | ! 3.4. Calculation of the tendencies of the vertical transport from the surface up to the top of the atmosphere |
|---|
| 570 | ! ********************************************************************** |
|---|
| 571 | DO l=1,nlayer |
|---|
| 572 | dqvl_topdust_mass(ig,l) = (zq_topdust_mass(ig,l)- & |
|---|
| 573 | zq0(ig,l,igcm_topdust_mass)) /ptimestep |
|---|
| 574 | dqvl_topdust_number(ig,l) = (zq_topdust_number(ig,l)- & |
|---|
| 575 | zq0(ig,l,igcm_topdust_number)) /ptimestep |
|---|
| 576 | dqvl_dust_mass(ig,l) = (zq_dust_mass(ig,l)-zq0(ig,l,igcm_dust_mass)) /ptimestep |
|---|
| 577 | dqvl_dust_number(ig,l) = (zq_dust_number(ig,l)-zq0(ig,l,igcm_dust_number)) /ptimestep |
|---|
| 578 | ENDDO |
|---|
| 579 | |
|---|
| 580 | ENDIF ! ((mu0(ig) .gt. mu0lim)... |
|---|
| 581 | ENDDO ! ig=1,ngrid |
|---|
| 582 | |
|---|
| 583 | ! ---------------------------------------------------------------------- |
|---|
| 584 | ! 4. Detrainment: topdust is converted to background dust |
|---|
| 585 | ! ---------------------------------------------------------------------- |
|---|
| 586 | |
|---|
| 587 | !! ********************************************************************** |
|---|
| 588 | !! 4.1 Compute the coefficient of detrainment |
|---|
| 589 | !! ********************************************************************** |
|---|
| 590 | DO ig=1,ngrid |
|---|
| 591 | !DO l=lwmax(ig),nlayer-1 |
|---|
| 592 | DO l=1,nlayer!-1 |
|---|
| 593 | !! Detrainment during the day |
|---|
| 594 | IF ( (mu0(ig) .gt. mu0lim) .and. (zq_topdust_mass(ig,l) .gt. zq_dust_mass(ig,l)*0.01) .and. (contains_mons(ig)) ) THEN |
|---|
| 595 | coefdetrain(ig,l)=1.*( rhobarz(ig,l+1)*abs(wfin(ig,l+1)) - rhobarz(ig,l)*abs(wfin(ig,l)) ) / masse(ig,l) |
|---|
| 596 | !! Detrainment when abs(w(l)) > abs(w(l+1)), i.e. coefdetrain < 0 |
|---|
| 597 | if ( coefdetrain(ig,l).lt.0.) then |
|---|
| 598 | dqdet_topdust_mass(ig,l)=-(1.-exp(coefdetrain(ig,l)*ptimestep))*zq_topdust_mass(ig,l)/ptimestep |
|---|
| 599 | dqdet_topdust_number(ig,l)=-(1.-exp(coefdetrain(ig,l)*ptimestep))*zq_topdust_number(ig,l)/ptimestep |
|---|
| 600 | ! for diagnostics |
|---|
| 601 | detrain_rate(ig,l)=(1.-exp(coefdetrain(ig,l)*ptimestep)) |
|---|
| 602 | !! One cannot detrain more than available topdust |
|---|
| 603 | if ( (abs(dqdet_topdust_mass(ig,l)).gt.zq_topdust_mass(ig,l)/ptimestep) .or. (abs(dqdet_topdust_number(ig,l)).gt.zq_topdust_number(ig,l)/ptimestep) ) then |
|---|
| 604 | dqdet_topdust_mass(ig,l)=-zq_topdust_mass(ig,l)/ptimestep |
|---|
| 605 | dqdet_topdust_number(ig,l)=-zq_topdust_number(ig,l)/ptimestep |
|---|
| 606 | detrain_rate(ig,l)=1. |
|---|
| 607 | endif |
|---|
| 608 | !else!entrainment |
|---|
| 609 | ! dqdet_topdust_mass(ig,l)=-(1.-exp(coefdetrain(ig,l)*ptimestep))*zq_dust_mass(ig,l)/ptimestep |
|---|
| 610 | ! dqdet_topdust_number(ig,l)=-(1.-exp(coefdetrain(ig,l)*ptimestep))*zq_dust_number(ig,l)/ptimestep |
|---|
| 611 | endif |
|---|
| 612 | !! Full detrainment imposed during the night or when topdust leaves its origin column (contains_mons=False) |
|---|
| 613 | ELSE |
|---|
| 614 | dqdet_topdust_mass(ig,l)=-zq_topdust_mass(ig,l)/ptimestep |
|---|
| 615 | dqdet_topdust_number(ig,l)=-zq_topdust_number(ig,l)/ptimestep |
|---|
| 616 | detrain_rate(ig,l)=1. |
|---|
| 617 | ENDIF |
|---|
| 618 | ENDDO ! DO l=1,nlayer |
|---|
| 619 | ENDDO ! DO ig=1,ngrid |
|---|
| 620 | |
|---|
| 621 | ! ---------------------------------------------------------------------- |
|---|
| 622 | ! 5. Final tendencies: vertical transport + detrainment |
|---|
| 623 | ! ---------------------------------------------------------------------- |
|---|
| 624 | |
|---|
| 625 | DO ig=1,ngrid |
|---|
| 626 | DO l=1,nlayer |
|---|
| 627 | pdqtop(ig,l,igcm_topdust_mass)=dqdet_topdust_mass(ig,l) & |
|---|
| 628 | +dqvl_topdust_mass(ig,l) |
|---|
| 629 | pdqtop(ig,l,igcm_topdust_number)=dqdet_topdust_number(ig,l) & |
|---|
| 630 | +dqvl_topdust_number(ig,l) |
|---|
| 631 | pdqtop(ig,l,igcm_dust_mass)= -dqdet_topdust_mass(ig,l) & |
|---|
| 632 | +dqvl_dust_mass(ig,l) |
|---|
| 633 | pdqtop(ig,l,igcm_dust_number)= -dqdet_topdust_number(ig,l) & |
|---|
| 634 | +dqvl_dust_number(ig,l) |
|---|
| 635 | ENDDO ! DO l=1,nlayer |
|---|
| 636 | ENDDO ! DO ig=1,ngrid |
|---|
| 637 | |
|---|
| 638 | |
|---|
| 639 | ! ********************************************************************** |
|---|
| 640 | ! WRITEDIAGFI |
|---|
| 641 | ! ********************************************************************** |
|---|
| 642 | IF (ngrid.eq.1) THEN |
|---|
| 643 | CALL WRITEDIAGFI(ngrid,'wup_top', & |
|---|
| 644 | 'wup_top','',1,wup(:,:)) |
|---|
| 645 | CALL WRITEDIAGFI(ngrid,'wrad_top', & |
|---|
| 646 | 'wrad_top','',1,wrad(:,:)) |
|---|
| 647 | CALL WRITEDIAGFI(ngrid,'wfin_top', & |
|---|
| 648 | 'wfin_top','',1,wfin(:,:)) |
|---|
| 649 | CALL WRITEDIAGFI(ngrid,'wlwmax_top', & |
|---|
| 650 | 'wlwmax_top','',0,wup(:,lwmax(1))) |
|---|
| 651 | CALL WRITEDIAGFI(ngrid,'w0_top', & |
|---|
| 652 | 'w0_top','',0,wup(:,lsummit(1)+1)) |
|---|
| 653 | CALL WRITEDIAGFI(ngrid,'alpha_hmons', & |
|---|
| 654 | 'alpha_hmons','',0,alpha_hmons) |
|---|
| 655 | CALL WRITEDIAGFI(ngrid,'zt_top', & |
|---|
| 656 | 'zt_top','',1,t_top(:,:)) |
|---|
| 657 | CALL WRITEDIAGFI(ngrid,'dt_top', & |
|---|
| 658 | 'dt_top','',0,dt_top(:)) |
|---|
| 659 | CALL WRITEDIAGFI(ngrid,'envtemp', & |
|---|
| 660 | 'envtemp','',1,zt(:,:)) |
|---|
| 661 | CALL WRITEDIAGFI(ngrid,'t_env', & |
|---|
| 662 | 't_env','',1,t_env(:,:)) |
|---|
| 663 | CALL WRITEDIAGFI(ngrid,'newzt_top', & |
|---|
| 664 | 'newzt_top','',1,newzt(:,:)) |
|---|
| 665 | CALL WRITEDIAGFI(ngrid,'deltahr_top', & |
|---|
| 666 | 'deltahr_top','',1,deltahr(:,:)) |
|---|
| 667 | CALL WRITEDIAGFI(ngrid,'rholev', & |
|---|
| 668 | 'rholev','',1,rho(:,:)) |
|---|
| 669 | CALL WRITEDIAGFI(ngrid,'rhobarz', & |
|---|
| 670 | 'rhobarz','',1,rhobarz(:,:)) |
|---|
| 671 | CALL WRITEDIAGFI(ngrid,'zlsummit', & |
|---|
| 672 | 'zlsummit','',0,pzlay(:,lsummit(1))) |
|---|
| 673 | CALL WRITEDIAGFI(ngrid,'zlwmax', & |
|---|
| 674 | 'zlwmax','',0,pzlay(:,lwmax(1))) |
|---|
| 675 | CALL WRITEDIAGFI(ngrid,'pzlev_top', & |
|---|
| 676 | 'pzlev_top','',1,pzlev(:,:)) |
|---|
| 677 | CALL WRITEDIAGFI(ngrid,'coefdetrain', & |
|---|
| 678 | 'coefdetrain','',1,coefdetrain(:,:)) |
|---|
| 679 | CALL WRITEDIAGFI(ngrid,'zdtvert', & |
|---|
| 680 | 'zdtvert','',1,zdtvert(:,:)) |
|---|
| 681 | CALL WRITEDIAGFI(ngrid,'entr', & |
|---|
| 682 | 'entr','',0,entr(:)) |
|---|
| 683 | ELSE |
|---|
| 684 | CALL WRITEDIAGFI(ngrid,'wup_top', & |
|---|
| 685 | 'wup_top','',3,wup(:,:)) |
|---|
| 686 | CALL WRITEDIAGFI(ngrid,'wup2_top', & |
|---|
| 687 | 'wup2_top','',3,wup2(:,:)) |
|---|
| 688 | CALL WRITEDIAGFI(ngrid,'wrad_top', & |
|---|
| 689 | 'wrad_top','',3,wrad(:,:)) |
|---|
| 690 | CALL WRITEDIAGFI(ngrid,'wfin_top', & |
|---|
| 691 | 'wfin_top','',3,wfin(:,:)) |
|---|
| 692 | CALL WRITEDIAGFI(ngrid,'hmons', & |
|---|
| 693 | 'hmons','',2,hmons) |
|---|
| 694 | CALL WRITEDIAGFI(ngrid,'hsummit', & |
|---|
| 695 | 'hsummit','',2,hsummit) |
|---|
| 696 | CALL WRITEDIAGFI(ngrid,'alpha_hmons', & |
|---|
| 697 | 'alpha_hmons','',2,alpha_hmons) |
|---|
| 698 | CALL WRITEDIAGFI(ngrid,'zt_top', & |
|---|
| 699 | 'zt_top','',3,t_top(:,:)) |
|---|
| 700 | CALL WRITEDIAGFI(ngrid,'ztemp', & |
|---|
| 701 | 'envtemp','',3,zt(:,:)) |
|---|
| 702 | CALL WRITEDIAGFI(ngrid,'zt_env', & |
|---|
| 703 | 't_env','',3,t_env(:,:)) |
|---|
| 704 | CALL WRITEDIAGFI(ngrid,'zdtvert_top', & |
|---|
| 705 | 'zdtvert_top','',3,zdtvert(:,:)) |
|---|
| 706 | CALL WRITEDIAGFI(ngrid,'newzt_top', & |
|---|
| 707 | 'newzt_top','',3,newzt(:,:)) |
|---|
| 708 | CALL WRITEDIAGFI(ngrid,'deltahr_top', & |
|---|
| 709 | 'deltahr_top','',3,deltahr(:,:)) |
|---|
| 710 | CALL WRITEDIAGFI(ngrid,'rhobarz', & |
|---|
| 711 | 'rhobarz','',3,rhobarz(:,:)) |
|---|
| 712 | CALL WRITEDIAGFI(ngrid,'zlwmax', & |
|---|
| 713 | 'zlwmax','',2,zlaywmax(:)) |
|---|
| 714 | CALL WRITEDIAGFI(ngrid,'coefdetrain', & |
|---|
| 715 | 'coefdetrain','',3,coefdetrain(:,:)) |
|---|
| 716 | CALL WRITEDIAGFI(ngrid,'detrain_rate', & |
|---|
| 717 | 'detrain_rate', & |
|---|
| 718 | '',3,detrain_rate(:,:)) |
|---|
| 719 | |
|---|
| 720 | ENDIF |
|---|
| 721 | |
|---|
| 722 | END SUBROUTINE topmons |
|---|
| 723 | |
|---|
| 724 | !======================================================================= |
|---|
| 725 | |
|---|
| 726 | ! ********************************************************************** |
|---|
| 727 | ! Subroutine to determine both temperature profiles above and near |
|---|
| 728 | ! a sub-grid scale mountain |
|---|
| 729 | !*********************************************************************** |
|---|
| 730 | SUBROUTINE t_topmons(nlayer,summit,hsummit,hmons,zt,zlay,t_top,dt_top,t_env,lsummit,lmons) |
|---|
| 731 | |
|---|
| 732 | USE comcstfi_h, only: g,cpp |
|---|
| 733 | |
|---|
| 734 | IMPLICIT NONE |
|---|
| 735 | |
|---|
| 736 | !-------------------------------------------------------- |
|---|
| 737 | ! Input Variables |
|---|
| 738 | !-------------------------------------------------------- |
|---|
| 739 | integer,intent(in) :: nlayer |
|---|
| 740 | real,intent(in) :: summit ! Height of the mountain summit |
|---|
| 741 | real,intent(in) :: hmons ! Height of the mountain slope |
|---|
| 742 | real,intent(in) :: hsummit ! Height of the mountain summit from the GCM surface |
|---|
| 743 | real,intent(in) :: zt(nlayer) ! large scale temperature profile |
|---|
| 744 | real,intent(in) :: zlay(nlayer) ! height at the middle of each layer |
|---|
| 745 | !-------------------------------------------------------- |
|---|
| 746 | ! Output Variables |
|---|
| 747 | !-------------------------------------------------------- |
|---|
| 748 | real,intent(out) :: t_top(nlayer) ! temperature above the mountain |
|---|
| 749 | real,intent(out) :: dt_top ! temperature difference between the slope and |
|---|
| 750 | ! the environment |
|---|
| 751 | real,intent(out) :: t_env(nlayer) ! temperature of the environment next to the mountain |
|---|
| 752 | integer,intent(out) :: lsummit ! layer reached by the summit height in the GCM |
|---|
| 753 | integer,intent(out) :: lmons ! layer of the real temperature to be considered near to the mountain top |
|---|
| 754 | |
|---|
| 755 | !-------------------------------------------------------- |
|---|
| 756 | ! Local Variables |
|---|
| 757 | !-------------------------------------------------------- |
|---|
| 758 | integer l,ll |
|---|
| 759 | real zlmons |
|---|
| 760 | |
|---|
| 761 | ! ********************************************************************** |
|---|
| 762 | t_env(:)=zt(:) |
|---|
| 763 | t_top(:)=zt(:) |
|---|
| 764 | dt_top=0. |
|---|
| 765 | lsummit=1 |
|---|
| 766 | lmons=1 |
|---|
| 767 | |
|---|
| 768 | !! Layer where the dust is injected |
|---|
| 769 | do while (hsummit .ge. zlay(lsummit)) |
|---|
| 770 | !! highest layer reached by the mountain (we don't interpolate and define zlay(lsummit) as the top of the mountain) |
|---|
| 771 | lsummit=lsummit+1 |
|---|
| 772 | enddo |
|---|
| 773 | !! temperature above the mountain is set to the surface temperature |
|---|
| 774 | t_top(lsummit)=zt(1) |
|---|
| 775 | do l=lsummit,nlayer-1 |
|---|
| 776 | !! temperature above the mountain follows the adiabatic |
|---|
| 777 | t_top(l+1)=t_top(l)-(zlay(l+1)-zlay(l))*g/cpp |
|---|
| 778 | enddo |
|---|
| 779 | |
|---|
| 780 | !! Layer where to take the temperature above the mountain |
|---|
| 781 | do while (hmons .ge. zlay(lmons)) |
|---|
| 782 | lmons=lmons+1 |
|---|
| 783 | enddo |
|---|
| 784 | !! Real temperature of the environment near to the mountain is t_env(lsummit)=zt(lmons) |
|---|
| 785 | !! (More simple and makes no real difference is to impose t_env(:)=zt(:)) |
|---|
| 786 | if (lmons.gt.lsummit) then |
|---|
| 787 | t_env(lsummit)=zt(lmons) |
|---|
| 788 | zlmons=zlay(lmons) |
|---|
| 789 | do l=0,nlayer-lsummit-1!2 |
|---|
| 790 | zlmons=zlmons+(zlay(lsummit+l+1)-zlay(lsummit+l)) |
|---|
| 791 | ll=lmons |
|---|
| 792 | do while (zlmons .ge. zlay(ll) .and. ll.lt.nlayer ) |
|---|
| 793 | ll=ll+1 |
|---|
| 794 | enddo |
|---|
| 795 | ll=ll-1 |
|---|
| 796 | !! Interpolation above lsummit |
|---|
| 797 | t_env(lsummit+l+1)= zt(ll) + ((zlmons-zlay(ll))*(zt(ll+1)-zt(ll))/(zlay(ll+1)-zlay(ll))) |
|---|
| 798 | enddo |
|---|
| 799 | t_env(nlayer)=zt(nlayer) |
|---|
| 800 | endif ! (lmons.gt.lsummit) |
|---|
| 801 | do l=1,nlayer |
|---|
| 802 | !! temperature above the mountain follows the adiabatic up to reach the environment temperature |
|---|
| 803 | if (t_top(l).le.t_env(l)) then |
|---|
| 804 | t_top(l)=t_env(l) |
|---|
| 805 | endif |
|---|
| 806 | enddo |
|---|
| 807 | !! Temperature delta created at the top of the mountain |
|---|
| 808 | dt_top = t_top(lsummit)-t_env(lsummit) |
|---|
| 809 | |
|---|
| 810 | END SUBROUTINE t_topmons |
|---|
| 811 | |
|---|
| 812 | !======================================================================= |
|---|
| 813 | ! ********************************************************************** |
|---|
| 814 | ! Subroutine to determine the vertical transport with |
|---|
| 815 | ! a Van Leer advection scheme (copied from the sedimentation scheme --> see vlz_fi.F) |
|---|
| 816 | !*********************************************************************** |
|---|
| 817 | SUBROUTINE van_leer(nlay,q,pente_max,lwmax,masse,w,wq,masse_pbl,dqm_pbl,alpha_hmons,qbar_pbl) |
|---|
| 818 | |
|---|
| 819 | IMPLICIT NONE |
|---|
| 820 | |
|---|
| 821 | !-------------------------------------------------------- |
|---|
| 822 | ! Input/Output Variables |
|---|
| 823 | !-------------------------------------------------------- |
|---|
| 824 | INTEGER,INTENT(IN) :: nlay ! number of atmospheric layers |
|---|
| 825 | REAL,INTENT(IN) :: masse(nlay) ! mass of the layer Dp/g |
|---|
| 826 | REAL,INTENT(IN) :: pente_max != 2 conseillee |
|---|
| 827 | INTEGER,INTENT(IN) :: lwmax ! layer of dust injection above the mountain, where the vertical velocity is maximum |
|---|
| 828 | REAL,INTENT(INOUT) :: w(nlay) ! atmospheric mass "transferred" at each timestep (kg.m-2) |
|---|
| 829 | REAL,INTENT(INOUT) :: wq(nlay+1) |
|---|
| 830 | REAL,INTENT(INOUT) :: q(nlay) ! mixing ratio (kg/kg) |
|---|
| 831 | REAL,INTENT(IN) :: masse_pbl |
|---|
| 832 | REAL,INTENT(IN) :: dqm_pbl |
|---|
| 833 | REAL,INTENT(IN) :: alpha_hmons |
|---|
| 834 | REAL,INTENT(IN) :: qbar_pbl |
|---|
| 835 | |
|---|
| 836 | !-------------------------------------------------------- |
|---|
| 837 | ! Local Variables |
|---|
| 838 | !-------------------------------------------------------- |
|---|
| 839 | |
|---|
| 840 | INTEGER i,l,j,ii |
|---|
| 841 | REAL dzq(nlay),dzqw(nlay),adzqw(nlay),dzqmax |
|---|
| 842 | REAL sigw, Mtot, MQtot |
|---|
| 843 | INTEGER m |
|---|
| 844 | |
|---|
| 845 | ! ********************************************************************** |
|---|
| 846 | ! Mixing ratio vertical gradient at the levels |
|---|
| 847 | ! ********************************************************************** |
|---|
| 848 | do l=lwmax+1,nlay !l=2,nlay |
|---|
| 849 | dzqw(l)=q(l-1)-q(l) |
|---|
| 850 | adzqw(l)=abs(dzqw(l)) |
|---|
| 851 | enddo |
|---|
| 852 | |
|---|
| 853 | do l=lwmax+1,nlay-1 !l=2,nlay-1 |
|---|
| 854 | if(dzqw(l)*dzqw(l+1).gt.0.) then |
|---|
| 855 | dzq(l)=0.5*(dzqw(l)+dzqw(l+1)) |
|---|
| 856 | else |
|---|
| 857 | dzq(l)=0. |
|---|
| 858 | endif |
|---|
| 859 | dzqmax=pente_max*min(adzqw(l),adzqw(l+1)) |
|---|
| 860 | dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l)) |
|---|
| 861 | enddo |
|---|
| 862 | |
|---|
| 863 | dzq(lwmax)=0. !! dzq(1)=0. => Attention: in the case of van leer only above the mountain it is very important to initialize the transport from lwmax !! |
|---|
| 864 | dzq(nlay)=0. |
|---|
| 865 | ! ********************************************************************** |
|---|
| 866 | ! Vertical advection |
|---|
| 867 | ! ********************************************************************** |
|---|
| 868 | |
|---|
| 869 | !! No flux at the model top: |
|---|
| 870 | wq(nlay+1)=0. |
|---|
| 871 | |
|---|
| 872 | !! Mass flux injected at lwmax |
|---|
| 873 | wq(lwmax) = -dqm_pbl/alpha_hmons ! dqm_pbl = alpha_hmons x wup(wmax) x rho(lwmax) x ptimestep x qbar_pbl |
|---|
| 874 | ! /alpha_hmons because the variable considered is mr_topdust |
|---|
| 875 | do l = lwmax,nlay-1 |
|---|
| 876 | |
|---|
| 877 | ! 1) Compute wq where w < 0 (up) (UPWARD TRANSPORT) |
|---|
| 878 | ! =============================== |
|---|
| 879 | |
|---|
| 880 | if (w(l+1).le.0) then |
|---|
| 881 | ! Regular scheme (transfered mass < 1 layer) |
|---|
| 882 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 883 | if(-w(l+1).le.masse(l))then |
|---|
| 884 | sigw=w(l+1)/masse(l) |
|---|
| 885 | wq(l+1)=w(l+1)*(q(l)-0.5*(1.+sigw)*dzq(l)) |
|---|
| 886 | !!------------------------------------------------------- |
|---|
| 887 | ! The following part should not be needed in the |
|---|
| 888 | ! case of an integration over subtimesteps |
|---|
| 889 | !!------------------------------------------------------- |
|---|
| 890 | !! Extended scheme (transfered mass > 1 layer) |
|---|
| 891 | !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 892 | !! else |
|---|
| 893 | !! m = l-1 |
|---|
| 894 | !! Mtot = masse(m+1) |
|---|
| 895 | !! MQtot = masse(m+1)*q(m+1) |
|---|
| 896 | !! if (m.lt.lwmax)goto 77!if (m.le.0)goto 77 |
|---|
| 897 | !! do while(-w(l+1).gt.(Mtot+masse(m))) |
|---|
| 898 | !! m=m-1 |
|---|
| 899 | !! Mtot = Mtot + masse(m+1) |
|---|
| 900 | !! MQtot = MQtot + masse(m+1)*q(m+1) |
|---|
| 901 | !! if (m.lt.lwmax)goto 77!if (m.le.0)goto 77 |
|---|
| 902 | !! end do |
|---|
| 903 | !! 77 continue |
|---|
| 904 | !! |
|---|
| 905 | !! if (m.gt.lwmax) then!if (m.gt.0) then |
|---|
| 906 | !! sigw=(w(l+1)+Mtot)/masse(m) |
|---|
| 907 | !! wq(l+1)= -(MQtot + (-w(l+1)-Mtot)* & |
|---|
| 908 | !! (q(m)-0.5*(1.+sigw)*dzq(m))) |
|---|
| 909 | !! else if ((-w(l+1).gt.(Mtot))) then |
|---|
| 910 | !! Mtot=Mtot+masse_pbl!*alpha_hmons |
|---|
| 911 | !! MQtot=MQtot+masse_pbl*qbar_pbl!*alpha_hmons |
|---|
| 912 | !! sigw=(w(l+1)+Mtot)/masse(m) |
|---|
| 913 | !! wq(l+1)= -(MQtot + (-w(l+1)-Mtot)*qbar_pbl) |
|---|
| 914 | !! else |
|---|
| 915 | !! w(l+1) = -Mtot |
|---|
| 916 | !! wq(l+1) = -MQtot |
|---|
| 917 | !! end if |
|---|
| 918 | !! |
|---|
| 919 | endif |
|---|
| 920 | |
|---|
| 921 | ! 2) Compute wq where w > 0 (down) (DOWNWARD TRANSPORT) |
|---|
| 922 | ! =============================== |
|---|
| 923 | |
|---|
| 924 | else if (w(l).gt.0.) then |
|---|
| 925 | ! Regular scheme (transfered mass < 1 layer) |
|---|
| 926 | ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 927 | if(w(l).le.masse(l))then |
|---|
| 928 | sigw=w(l)/masse(l) |
|---|
| 929 | wq(l)=w(l)*(q(l)+0.5*(1.-sigw)*dzq(l)) |
|---|
| 930 | !!------------------------------------------------------- |
|---|
| 931 | ! The following part should not be needed in the |
|---|
| 932 | ! case of an integration over subtimesteps |
|---|
| 933 | !!------------------------------------------------------- |
|---|
| 934 | !! Extended scheme (transfered mass > 1 layer) |
|---|
| 935 | !! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 936 | !! else |
|---|
| 937 | !! m=l |
|---|
| 938 | !! Mtot = masse(m) |
|---|
| 939 | !! MQtot = masse(m)*q(m) |
|---|
| 940 | !! if(m.ge.nlay)goto 88 |
|---|
| 941 | !! do while(w(l).gt.(Mtot+masse(m+1))) |
|---|
| 942 | !! m=m+1 |
|---|
| 943 | !! Mtot = Mtot + masse(m) |
|---|
| 944 | !! MQtot = MQtot + masse(m)*q(m) |
|---|
| 945 | !! if(m.ge.nlay)goto 88 |
|---|
| 946 | !! end do |
|---|
| 947 | !! 88 continue |
|---|
| 948 | !! if (m.lt.nlay) then |
|---|
| 949 | !! sigw=(w(l)-Mtot)/masse(m+1) |
|---|
| 950 | !! wq(l)=(MQtot + (w(l)-Mtot)* & |
|---|
| 951 | !! (q(m+1)+0.5*(1.-sigw)*dzq(m+1)) ) |
|---|
| 952 | !! else |
|---|
| 953 | !! w(l) = Mtot |
|---|
| 954 | !! wq(l) = MQtot |
|---|
| 955 | !! end if |
|---|
| 956 | end if ! (w(l).le.masse(l)) |
|---|
| 957 | |
|---|
| 958 | end if ! (w(l+1).le.0) |
|---|
| 959 | |
|---|
| 960 | enddo ! l = lwmax+1,nlay-1 |
|---|
| 961 | |
|---|
| 962 | ! ********************************************************************** |
|---|
| 963 | ! Mixing ratio update after the vertical transport |
|---|
| 964 | ! ********************************************************************** |
|---|
| 965 | do l = lwmax,nlay-1 !l = 1,nlay-1 ! loop different than when w>0 |
|---|
| 966 | |
|---|
| 967 | if ( (wq(l+1)-wq(l)) .lt. -(masse(l)*q(l)) .and. (abs(wq(l+1)).gt. 0.) .and. (abs(wq(l)).gt. 0.) ) then |
|---|
| 968 | wq(l+1) = wq(l)-masse(l)*q(l) |
|---|
| 969 | end if |
|---|
| 970 | |
|---|
| 971 | q(l)=q(l) + (wq(l+1)-wq(l))/masse(l) |
|---|
| 972 | |
|---|
| 973 | enddo |
|---|
| 974 | |
|---|
| 975 | |
|---|
| 976 | end subroutine van_leer |
|---|
| 977 | |
|---|
| 978 | !******************************************************************************** |
|---|
| 979 | subroutine topmons_setup(ngrid) |
|---|
| 980 | ! Purpose: |
|---|
| 981 | ! 1) Fill the logical array contains_mons(:), |
|---|
| 982 | ! with contains_mons(ig)=True if there is |
|---|
| 983 | ! a mountain in the mesh ig |
|---|
| 984 | ! 2) Compute alpha_hmons(:) and hsummit(:) |
|---|
| 985 | use surfdat_h,only: phisfi,hmons,summit,base,& |
|---|
| 986 | alpha_hmons,hsummit,contains_mons |
|---|
| 987 | use comcstfi_h,only: pi,g |
|---|
| 988 | use planetwide_mod,only: planetwide_maxval, planetwide_minval |
|---|
| 989 | use geometry_mod,only: longitude_deg,latitude_deg,& |
|---|
| 990 | boundslon,boundslat !boundslon/lat(ngrid,4) : |
|---|
| 991 | ! |------------------------------| |
|---|
| 992 | ! |north_west=2 north_east=1| |
|---|
| 993 | ! | | |
|---|
| 994 | ! | (ig) | |
|---|
| 995 | ! | | |
|---|
| 996 | ! |south_west=3 south_east=4| |
|---|
| 997 | ! |------------------------------| |
|---|
| 998 | |
|---|
| 999 | implicit none |
|---|
| 1000 | integer,intent(in) :: ngrid |
|---|
| 1001 | |
|---|
| 1002 | ! Local variables |
|---|
| 1003 | integer,parameter :: ntop_max = 19 ! total number of mounts written in the hmons list |
|---|
| 1004 | real lon_top(ntop_max),lat_top(ntop_max) ! coordinates of the mounts (in deg) |
|---|
| 1005 | ! Mountains list : |
|---|
| 1006 | ! Olympus Mons,Ascraeus Mons,Elysium Mons,Arsia Mons,Pavonis Mons, |
|---|
| 1007 | ! Hecates Tholus,Tharsis Tholus,Ceraunius Tholus,Alba Mons,Apollinaris Mons, |
|---|
| 1008 | ! Albor Tholus,Biblis Tholus,Anseris Mons,Ulysses Tholus,Aeolis Mons, |
|---|
| 1009 | ! Euripus Mons,Hadriacus Mons,Tyrrhenus Mons,Uranius Mons |
|---|
| 1010 | ! |
|---|
| 1011 | ! NB: in 64x48 horiz. resolution, Biblis Tholus & Ulysses Tholus fall in the |
|---|
| 1012 | ! same mesh, hence only Biblis Tholus is kept by the alpha_hmons computation |
|---|
| 1013 | data lon_top/-134.,-104.5,146.9,-121.1,-113.4,& |
|---|
| 1014 | 150.2,-90.8,-97.4,-109.6,174.4,& |
|---|
| 1015 | 150.4,-124.6,86.6,-121.6,137.8,& |
|---|
| 1016 | 105.,91.8,106.5,-92.2/ |
|---|
| 1017 | data lat_top/18.4,11.8,24.8,-8.4,-0.8,& |
|---|
| 1018 | 31.8,13.4,24.,40.4,-9.3,& |
|---|
| 1019 | 18.8,2.6,-29.8,2.9,-5.4,& |
|---|
| 1020 | -44.8,-32.1,-21.1,26.8/ |
|---|
| 1021 | integer,parameter :: ntop = 19 ! the topmons scheme is limited to the first ntop mounts |
|---|
| 1022 | real :: boundslon_deg(4),boundslat_deg(4) |
|---|
| 1023 | real :: hmax,hmin |
|---|
| 1024 | integer :: ig,itop |
|---|
| 1025 | |
|---|
| 1026 | |
|---|
| 1027 | |
|---|
| 1028 | IF (ngrid.gt.1) THEN |
|---|
| 1029 | ! Sanity check |
|---|
| 1030 | if (ntop.gt.ntop_max) then |
|---|
| 1031 | call abort_physic("topmons_setup","Number of mountains ntop greater than ntop_max",1) |
|---|
| 1032 | endif |
|---|
| 1033 | |
|---|
| 1034 | ! Determine contains_mons |
|---|
| 1035 | contains_mons(:)=.false. |
|---|
| 1036 | |
|---|
| 1037 | do ig=1,ngrid |
|---|
| 1038 | boundslon_deg(:)=boundslon(ig,:)/pi*180. |
|---|
| 1039 | boundslat_deg(:)=boundslat(ig,:)/pi*180. |
|---|
| 1040 | |
|---|
| 1041 | do itop=1,ntop |
|---|
| 1042 | if ( (lon_top(itop).gt.boundslon_deg(2)).and.(lon_top(itop).lt.boundslon_deg(1)) & |
|---|
| 1043 | .and.(lat_top(itop).gt.boundslat_deg(3)).and.(lat_top(itop).lt.boundslat_deg(2)) ) then |
|---|
| 1044 | contains_mons(ig)=.true. |
|---|
| 1045 | write(*,*) "topmons_setup: Found a mount at:" |
|---|
| 1046 | write(*,*) "(",boundslon_deg(2),",",boundslat_deg(2),") (",boundslon_deg(1),",",boundslat_deg(1),")" |
|---|
| 1047 | write(*,*) " ((",lon_top(itop),",",lat_top(itop),"))" |
|---|
| 1048 | write(*,*) "(",boundslon_deg(3),",",boundslat_deg(3),") (",boundslon_deg(4),",",boundslat_deg(4),")" |
|---|
| 1049 | endif |
|---|
| 1050 | enddo |
|---|
| 1051 | enddo |
|---|
| 1052 | |
|---|
| 1053 | ! Compute alpha_hmons |
|---|
| 1054 | call planetwide_maxval(hmons,hmax) |
|---|
| 1055 | call planetwide_minval(hmons,hmin) |
|---|
| 1056 | do ig=1,ngrid |
|---|
| 1057 | if (contains_mons(ig)) then |
|---|
| 1058 | ! the mesh ig contains a mountain |
|---|
| 1059 | alpha_hmons(ig)= 0.5*(hmons(ig)-hmin)/(hmax-hmin) |
|---|
| 1060 | ! Sanity check |
|---|
| 1061 | if (alpha_hmons(ig).le.0) then |
|---|
| 1062 | call abort_physic("topmons_setup","ERROR: alpha_hmons cannot be <0 "//& |
|---|
| 1063 | "if the mesh contains a mountain. Please check your "//& |
|---|
| 1064 | "formula or your start files.",1) |
|---|
| 1065 | endif |
|---|
| 1066 | else |
|---|
| 1067 | ! the mesh ig doesn't contain a mountain |
|---|
| 1068 | alpha_hmons(ig)= 0 |
|---|
| 1069 | endif |
|---|
| 1070 | enddo |
|---|
| 1071 | |
|---|
| 1072 | ! Compute hsummit |
|---|
| 1073 | hsummit(:)=summit(:)-phisfi(:)/g |
|---|
| 1074 | |
|---|
| 1075 | ELSE ! 1D case |
|---|
| 1076 | hmin=0. |
|---|
| 1077 | hmax=23162.1 !set here the height of the sub-grid scale topography |
|---|
| 1078 | do ig=1,ngrid |
|---|
| 1079 | alpha_hmons(ig)= (hmons(ig)-hmin)/(hmax-hmin) !0.1*(hmons(ig)-hmin)/(hmax-hmin) |
|---|
| 1080 | print*,"1D, hmons=",hmons(ig),"alpha=",alpha_hmons(ig) |
|---|
| 1081 | enddo |
|---|
| 1082 | |
|---|
| 1083 | hsummit(:)=14000. |
|---|
| 1084 | ENDIF ! (ngrid.gt.1) |
|---|
| 1085 | |
|---|
| 1086 | end subroutine topmons_setup |
|---|
| 1087 | !******************************************************************************** |
|---|
| 1088 | |
|---|
| 1089 | END MODULE topmons_mod |
|---|
