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