source: trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90 @ 2673

Last change on this file since 2673 was 2643, checked in by abierjon, 3 years ago

Mars GCM:
Implementation of the IRabs-to-VISext dust scenario conversion, depending on the GCM effective radius :

  • new flag 'reff_driven_IRtoVIS_scenario', false by default, must be set to true to use this dynamic conversion (otherwise, the coefficient takes the constant value of 2.6, eqv to reff=1.5um, as before) A specific line is added in deftank/callphys.def.GCM6
  • this flag requires the 'dustiropacity' to be set to 'tes' to match the IR scenario's wavelength. 'mcs'-like dso diagnostics can only be produced by the use of the post-processing tool util/aeroptical.F90
  • the variable IRtoVIScoef is computed in aeropacity via the GCM CDOD ratio : tau_pref_gcm_VIS/tau_pref_gcm_IR (only at the first call to callradite during each physical timestep)
  • change read_dust_scenario into module read_dust_scenario_mod

AB

File size: 31.3 KB
RevLine 
[1974]1      MODULE rocketduststorm_mod
2
3      IMPLICIT NONE
4
5      REAL, SAVE, ALLOCATABLE :: dustliftday(:) ! dust lifting rate (s-1)
[2584]6
7!$OMP THREADPRIVATE(dustliftday)
[1974]8     
9      CONTAINS
10
11!=======================================================================
12! ROCKET DUST STORM - vertical transport and detrainment
13!=======================================================================
[2079]14! calculation of the vertical flux
[2201]15! call of van_leer : Van Leer transport scheme of the dust tracers
[2079]16! detrainement of stormdust in the background dust
[1974]17! -----------------------------------------------------------------------
[2079]18! Authors: M. Vals; C. Wang; F. Forget; T. Bertrand
[1974]19! Institution: Laboratoire de Meteorologie Dynamique (LMD) Paris, France
20! -----------------------------------------------------------------------
21
22      SUBROUTINE rocketduststorm(ngrid,nlayer,nq,ptime,ptimestep,      &
23                                 pq,pdqfi,pt,pdtfi,pplev,pplay,pzlev,  &
24                                 pzlay,pdtsw,pdtlw,                    &
25!             input for radiative transfer
26                                 clearatm,icount,zday,zls,             &
27                                 tsurf,igout,totstormfract,            &
[2417]28                                 tauscaling,dust_rad_adjust,           &
[2643]29                                 IRtoVIScoef,                          &
[1974]30!             input sub-grid scale cloud
31                                 clearsky,totcloudfrac,                &
[2199]32!             input sub-grid scale topography
[2628]33                                 nohmons,                              &
[2079]34!             output
[2246]35                                 pdqrds,wrad,dsodust,dsords,dsotop,    &
[2415]36                                 tau_pref_scenario,tau_pref_gcm)
[1974]37
[2079]38      USE tracer_mod, only: igcm_stormdust_mass,igcm_stormdust_number &
[1974]39                            ,igcm_dust_mass,igcm_dust_number          &
40                            ,rho_dust
41      USE comcstfi_h, only: r,g,cpp,rcp
[2079]42      USE dimradmars_mod, only: albedo,naerkind
43      USE comsaison_h, only: dist_sol,mu0,fract
44      USE surfdat_h, only: emis,co2ice,zmea, zstd, zsig, hmons
[2226]45      USE callradite_mod, only: callradite
[2079]46      IMPLICIT NONE
[1974]47
[2160]48      include "callkeys.h"
49
[1974]50!--------------------------------------------------------
51! Input Variables
52!--------------------------------------------------------
53
54      INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points
55      INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points
56      INTEGER, INTENT(IN) :: nq ! number of tracer species
57      REAL, INTENT(IN) :: ptime
58      REAL, INTENT(IN) :: ptimestep
59
60      REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) ! advected field nq
61      REAL, INTENT(IN) :: pdqfi(ngrid,nlayer,nq)! tendancy field pq
62      REAL, INTENT(IN) :: pt(ngrid,nlayer)      ! temperature at mid-layer (K)
63      REAL, INTENT(IN) :: pdtfi(ngrid,nlayer)   ! tendancy temperature at mid-layer (K)
64
65      REAL, INTENT(IN) :: pplay(ngrid,nlayer)     ! pressure at middle of the layers
66      REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) ! pressure at intermediate levels
67      REAL, INTENT(IN) :: pzlay(ngrid,nlayer)     ! altitude at the middle of the layers
68      REAL, INTENT(IN) :: pzlev(ngrid,nlayer+1)   ! altitude at layer boundaries
69
70      REAL, INTENT(IN) :: pdtsw(ngrid,nlayer)     ! (K/s) env
71      REAL, INTENT(IN) :: pdtlw(ngrid,nlayer)     ! (K/s) env
72
73!     input for second radiative transfer
[2079]74      LOGICAL, INTENT(IN) :: clearatm
[1974]75      INTEGER, INTENT(INOUT) :: icount
[2079]76      REAL, INTENT(IN) :: zday
77      REAL, INTENT(IN) :: zls
78      REAL, INTENT(IN) :: tsurf(ngrid)
79      INTEGER, INTENT(IN) :: igout
80      REAL, INTENT(IN) :: totstormfract(ngrid)
[2226]81      REAL, INTENT(INOUT) :: tauscaling(ngrid)
[2643]82      REAL,INTENT(INOUT) :: dust_rad_adjust(ngrid)                             
83      REAL,INTENT(INOUT) :: IRtoVIScoef(ngrid) ! NB: not modified by this call to callradite,
84                                               ! the OUT is just here because callradite needs it
[2079]85     
[2628]86!     subgrid scale water ice clouds
[1974]87      logical, intent(in) :: clearsky
[2199]88      real, intent(in) :: totcloudfrac(ngrid)
89
[2628]90!     subgrid scale topography
[2199]91      LOGICAL, INTENT(IN) :: nohmons
[1974]92 
93!--------------------------------------------------------
94! Output Variables
95!--------------------------------------------------------
96
97      REAL, INTENT(OUT) :: pdqrds(ngrid,nlayer,nq) ! tendancy field for dust when detraining
[2079]98      REAL, INTENT(OUT) :: wrad(ngrid,nlayer+1)   ! vertical speed within the rocket dust storm
[2415]99      REAL, INTENT(OUT) :: dsodust(ngrid,nlayer) ! density scaled opacity of env. dust
100      REAL, INTENT(OUT) :: dsords(ngrid,nlayer) ! density scaled opacity of storm dust
101      REAL, INTENT(OUT) :: dsotop(ngrid,nlayer) ! density scaled opacity of topmons dust
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
[1974]106!--------------------------------------------------------
107! Local variables
108!--------------------------------------------------------
[2201]109      INTEGER l,ig,iq,ll
[2079]110!     local variables from callradite.F       
[1974]111      REAL zdtlw1(ngrid,nlayer)    ! (K/s) storm
112      REAL zdtsw1(ngrid,nlayer)    ! (K/s) storm
113      REAL zt(ngrid,nlayer)        ! actual temperature at mid-layer (K)
[2091]114      REAL zdtvert(ngrid,nlayer)   ! dT/dz , lapse rate
115      REAL ztlev(ngrid,nlayer)     ! temperature at intermediate levels l+1/2 without last level
[1974]116
[2079]117      REAL zdtlw1_lev(nlayer),zdtsw1_lev(nlayer) ! rad. heating rate at intermediate levels l+1/2 for stormdust
118      REAL zdtlw_lev(nlayer),zdtsw_lev(nlayer)   ! rad. heating rate at intermediate levels l+1/2 for background dust
[1974]119
[2079]120      REAL zq_stormdust_mass(ngrid,nlayer) ! intermediate tracer stormdust mass
121      REAL zq_stormdust_number(ngrid,nlayer) ! intermediate tracer stormdust number
122      REAL zq_dust_mass(ngrid,nlayer)           ! intermediate tracer dust mass
123      REAL zq_dust_number(ngrid,nlayer)         ! intermediate tracer dust number
[1974]124
[2079]125      REAL mr_stormdust_mass(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (stormdust mass)
126      REAL mr_stormdust_number(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (stormdust number)
127      REAL mr_dust_mass(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (dust mass)
128      REAL mr_dust_number(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (sdust number)
129                   
130      REAL dqvl_stormdust_mass(ngrid,nlayer)    ! tendancy of vertical transport (stormdust mass)
131      REAL dqvl_stormdust_number(ngrid,nlayer)  ! tendancy of vertical transport (stormdust number)
132      REAL dqvl_dust_mass(ngrid,nlayer)    ! tendancy of vertical transport (dust mass)
133      REAL dqvl_dust_number(ngrid,nlayer)  ! tendancy of vertical transport (dust number)
134      REAL dqdet_stormdust_mass(ngrid,nlayer)   ! tendancy of detrainement (stormdust mass)
135      REAL dqdet_stormdust_number(ngrid,nlayer) ! tendancy of detrainement (stormdust number)
[1974]136
[2090]137      REAL masse_col(nlayer)     ! mass of atmosphere (kg/m2)
[1974]138      REAL zq(ngrid,nlayer,nq)   ! updated tracers
139     
[2090]140      REAL w(ngrid,nlayer)          ! air mass flux (calculated with the vertical wind velocity profile) used as input in Van Leer (kgair/m2)
141      REAL wqmass(ngrid,nlayer+1)   ! tracer (dust_mass) mass flux in Van Leer (kg/m2)
142      REAL wqnumber(ngrid,nlayer+1) ! tracer (dust_number) mass flux in Van Leer (kg/m2)
[1974]143
[2079]144      LOGICAL storm(ngrid)    ! true when there is a dust storm (if the opacity is high): trigger the rocket dust storm scheme
[2160]145      REAL detrain(ngrid,nlayer)  ! coefficient for detrainment : % of stormdust detrained
[2079]146      INTEGER scheme(ngrid)   ! triggered scheme
147           
148      REAL,PARAMETER:: coefmin =0.025 ! 0<coefmin<1 Minimum fraction of stormdust detrained 
[2201]149      REAL,PARAMETER:: wmin =0.25 ! stormdust detrainment if wrad < wmin 
[2079]150      REAL,PARAMETER:: wmax =10.   ! maximum vertical velocity of the rocket dust storms (m/s)
[2201]151
152!     subtimestep
153      INTEGER tsub
154      INTEGER nsubtimestep    !number of subtimestep when calling van_leer
155      REAL subtimestep        !ptimestep/nsubtimestep
156      REAL dtmax              !considered time needed for dust to cross one layer
157      REAL,PARAMETER:: secu=3.!3.      !coefficient on wspeed to avoid dust crossing many layers during subtimestep
158
[2079]159!     diagnostics
160      REAL lapserate(ngrid,nlayer)
161      REAL deltahr(ngrid,nlayer+1)
[1974]162   
[2079]163      LOGICAL,SAVE :: firstcall=.true.
[2616]164     
165!$OMP THREADPRIVATE(firstcall)
[2079]166
167!     variables for the radiative transfer
[1974]168      REAL  fluxsurf_lw1(ngrid)
169      REAL  fluxsurf_sw1(ngrid,2)
170      REAL  fluxtop_lw1(ngrid)
171      REAL  fluxtop_sw1(ngrid,2)
172      REAL  tau(ngrid,naerkind)
173      REAL  aerosol(ngrid,nlayer,naerkind)
174      REAL  taucloudtes(ngrid)
175      REAL  rdust(ngrid,nlayer)
176      REAL  rstormdust(ngrid,nlayer)
[2199]177      REAL  rtopdust(ngrid,nlayer)
[1974]178      REAL  rice(ngrid,nlayer)
179      REAL  nuice(ngrid,nlayer)
[2459]180      DOUBLE PRECISION  riceco2(ngrid,nlayer)
[2448]181      REAL  nuiceco2(ngrid,nlayer)
[1974]182
183      ! **********************************************************************
184      ! **********************************************************************
[2079]185      ! Rocket dust storm parametrization to reproduce the detached dust layers
186      ! during the dust storm season:
[1974]187      !     The radiative warming due to the presence of storm dust is
188      !     balanced by the adiabatic cooling. The tracer "storm dust" 
189      !     is transported by the upward/downward flow.
190      ! **********************************************************************
191      ! **********************************************************************     
192      !!    1. Radiative transfer in storm dust
193      !!    2. Compute vertical velocity for storm dust
[2079]194      !!      case 1 storm = false: nothing to do
195      !!      case 2 rocket dust storm (storm=true)
196      !!    3. Vertical transport (Van Leer)
[1974]197      !!    4. Detrainment
198      ! **********************************************************************
199      ! **********************************************************************
200
201     
202      ! **********************************************************************
[2079]203      ! Initializations
[1974]204      ! **********************************************************************
205      storm(:)=.false.
206      pdqrds(:,:,:) = 0.
[2079]207      mr_dust_mass(:,:)=0.
208      mr_dust_number(:,:)=0.
209      mr_stormdust_mass(:,:)=0.
210      mr_stormdust_number(:,:)=0.
211      dqvl_dust_mass(:,:)=0.
212      dqvl_dust_number(:,:)=0.
213      dqvl_stormdust_mass(:,:)=0.
214      dqvl_stormdust_number(:,:)=0.
215      dqdet_stormdust_mass(:,:)=0.
216      dqdet_stormdust_number(:,:)=0.
217      wrad(:,:)=0.
[2090]218      w(:,:)=0.
219      wqmass(:,:)=0.
220      wqnumber(:,:)=0.
[2091]221      zdtvert(:,:)=0.
[1974]222      lapserate(:,:)=0.
223      deltahr(:,:)=0.
224      scheme(:)=0
[2160]225      detrain(:,:)=1.
[2083]226
[1974]227      !! no update for the stormdust tracer and temperature fields
228      !! because previous callradite was for background dust
229      zq(1:ngrid,1:nlayer,1:nq)=pq(1:ngrid,1:nlayer,1:nq)
230      zt(1:ngrid,1:nlayer)=pt(1:ngrid,1:nlayer)
231
232
[2079]233      zq_dust_mass(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_dust_mass)
234      zq_dust_number(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_dust_number)
235      zq_stormdust_mass(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_stormdust_mass)
236      zq_stormdust_number(1:ngrid,1:nlayer)=zq(1:ngrid,1:nlayer,igcm_stormdust_number)
237
238      ! *********************************************************************
239      ! 0. Check if there is a storm
240      ! *********************************************************************
241      DO ig=1,ngrid
[1974]242        storm(ig)=.false.
[2079]243        DO l=1,nlayer
244          IF (zq(ig,l,igcm_stormdust_mass) &
245          .gt. zq(ig,l,igcm_dust_mass)*(1.E-4)) THEN
[1974]246            storm(ig)=.true.
[2079]247            EXIT
248          ENDIF
249        ENDDO     
250      ENDDO
[1974]251     
252      ! *********************************************************************
253      ! 1. Call the second radiative transfer for stormdust, obtain the extra heating
254      ! *********************************************************************
[2199]255      CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,          &
256                 emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,      &
257                 zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1,     &
[2415]258                 fluxtop_sw1,tau_pref_scenario,tau_pref_gcm, &
[2643]259                 tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef, &
[2448]260                 taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice,rstormdust,rtopdust, &
[2628]261                 totstormfract,clearatm,dsords,dsotop,nohmons,&
[1974]262                 clearsky,totcloudfrac)
263
264      ! **********************************************************************
265      ! 2. Compute vertical velocity for storm dust
[2079]266      ! **********************************************************************
[1974]267        !! **********************************************************************
[2079]268        !! 2.1 Nothing to do when no storm
269        !!             no storm
270        DO ig=1,ngrid       
271          IF (.NOT.(storm(ig))) then
272            scheme(ig)=1
273            cycle
274          ENDIF ! IF (storm(ig))
275        ENDDO ! DO ig=1,ngrid                 
[1974]276       
[2079]277        !! **********************************************************************
278        !! 2.2 Calculation of the extra heating : computing heating rates
279        !! gradient at boundaries of each layer, start from surface
280        DO ig=1,ngrid
[2091]281          IF (storm(ig)) THEN
282
283            scheme(ig)=2
[1974]284     
[2091]285            !! computing heating rates gradient at boundraies of each layer
286            !! start from surface
287            zdtlw1_lev(1)=0.
288            zdtsw1_lev(1)=0.
289            zdtlw_lev(1)=0.
290            zdtsw_lev(1)=0.
291            ztlev(ig,1)=zt(ig,1)
[1974]292
[2091]293            DO l=1,nlayer-1
294              !! Calculation for the dust storm fraction
295              zdtlw1_lev(l+1)=(zdtlw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
296                           zdtlw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /  &
297                              (pzlay(ig,l+1)-pzlay(ig,l))
298           
299              zdtsw1_lev(l+1)=(zdtsw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
300                           zdtsw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /  &
301                              (pzlay(ig,l+1)-pzlay(ig,l))
302              !! Calculation for the background dust fraction
303              zdtlw_lev(l+1)=(pdtlw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+   &
304                           pdtlw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /   &
305                              (pzlay(ig,l+1)-pzlay(ig,l))
[2079]306           
[2091]307              zdtsw_lev(l+1)=(pdtsw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+   &
308                           pdtsw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /   &
309                              (pzlay(ig,l+1)-pzlay(ig,l))
[2079]310           
[2091]311              ztlev(ig,l+1)=(zt(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+          &
312                           zt(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /      &
313                              (pzlay(ig,l+1)-pzlay(ig,l))
[2201]314            ENDDO ! DO l=1,nlayer-1
[1974]315
[2201]316            !! This is the env. lapse rate
317            zdtvert(ig,1)=0.
318            DO l=1,nlayer-1
319              zdtvert(ig,l+1)=(ztlev(ig,l+1)-ztlev(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l))
320              lapserate(ig,l+1)=zdtvert(ig,l+1)
321            ENDDO
322
[2091]323            !! **********************************************************************
324            !! 2.3 Calculation of the vertical velocity : extra heating
325            !! balanced by adiabatic cooling
326           
[2079]327            DO l=1,nlayer
328              deltahr(ig,l)=(zdtlw1_lev(l)+zdtsw1_lev(l))  &
329                                          -(zdtlw_lev(l)+zdtsw_lev(l))
330              wrad(ig,l)=-deltahr(ig,l)/(g/cpp+   &
[2091]331                                         max(zdtvert(ig,l),-0.99*g/cpp))       
[2079]332              !! Limit vertical wind in case of lapse rate close to adiabatic
333              wrad(ig,l)=max(wrad(ig,l),-wmax)
334              wrad(ig,l)=min(wrad(ig,l),wmax)
335            ENDDO
[2091]336
[2079]337          ENDIF ! IF (storm(ig))
338        ENDDO ! DO ig=1,ngrid
[1974]339
340      ! **********************************************************************
341      ! 3. Vertical transport
342      ! **********************************************************************
[2079]343        !! **********************************************************************
344        !! 3.1 Transport of background dust + storm dust (concentrated)
345        DO ig=1,ngrid
346          IF (storm(ig)) THEN
347            DO l=1,nlayer
348              mr_dust_mass(ig,l) = zq_dust_mass(ig,l)
349              mr_dust_number(ig,l) = zq_dust_number(ig,l)
350              mr_stormdust_mass(ig,l) = zq_dust_mass(ig,l)+ &
351                                      zq_stormdust_mass(ig,l)/totstormfract(ig)
352              mr_stormdust_number(ig,l) = zq_dust_number(ig,l)+ &
353                                      zq_stormdust_number(ig,l)/totstormfract(ig)
354            ENDDO ! DO l=1,nlayer
355          ENDIF ! IF (storm(ig))
356        ENDDO ! DO ig=1,ngrid
[1974]357
[2079]358        DO ig=1,ngrid
359          IF (storm(ig)) THEN
[2201]360        !! **********************************************************************
361        !! 3.2 Compute the subtimestep to conserve the mass in the Van Leer transport
362            dtmax=ptimestep
363            DO l=2,nlayer
364              IF (wrad(ig,l).lt.0.) THEN
365                 dtmax=min(dtmax,(pzlev(ig,l)-pzlev(ig,l-1))/  &
366                                   (secu*abs(wrad(ig,l))))
367              ELSE IF (wrad(ig,l).gt.0.) then
368                 dtmax=min(dtmax,(pzlev(ig,l+1)-pzlev(ig,l))/  &
369                                   (secu*abs(wrad(ig,l))))
370              ENDIF
371            ENDDO
372            nsubtimestep= int(ptimestep/dtmax)
373            subtimestep=ptimestep/float(nsubtimestep)
374            !! Mass flux generated by wup in kg/m2
375            DO l=1,nlayer
376               w(ig,l)=wrad(ig,l)*pplev(ig,l)/(r*ztlev(ig,l)) &
377                        *subtimestep
378            ENDDO ! l=1,nlayer
379
380        !! **********************************************************************
381        !! 3.3 Vertical transport by a Van Leer scheme
[2079]382            !! Mass of atmosphere in the layer
383            DO l=1,nlayer
[2090]384               masse_col(l)=(pplev(ig,l)-pplev(ig,l+1))/g
[2079]385            ENDDO
[2201]386            !! Mass flux in kg/m2 if you are not using the subtimestep
387            !DO l=1,nlayer
388            !   w(ig,l)=wrad(ig,l)*(pplev(ig,l)/(r*ztlev(ig,l)))*ptimestep
389            !ENDDO
390            !! Loop over the subtimestep
391            DO tsub=1,nsubtimestep
392              !! Van Leer scheme
393              wqmass(ig,:)=0.
394              wqnumber(ig,:)=0.
395              CALL van_leer(nlayer,mr_stormdust_mass(ig,:),2.,   &
396                    masse_col,w(ig,:),wqmass(ig,:))
397              CALL van_leer(nlayer,mr_stormdust_number(ig,:),2.,  &
398                    masse_col,w(ig,:),wqnumber(ig,:))
399            ENDDO !tsub=...           
400             
[2079]401          ENDIF ! IF storm(ig)
402        ENDDO ! DO ig=1,ngrid 
[1974]403
[2079]404        !! **********************************************************************
[2201]405        !! 3.4 Re-calculation of the mixing ratios after vertical transport
[2079]406        DO ig=1,ngrid
407         IF (storm(ig)) THEN
408           DO l=1,nlayer
409           
410             !! General and "healthy" case
411             IF (mr_stormdust_mass(ig,l).ge.mr_dust_mass(ig,l)) THEN
412               zq_dust_mass(ig,l) = mr_dust_mass(ig,l)
413               zq_dust_number(ig,l) = mr_dust_number(ig,l)
414               zq_stormdust_mass(ig,l) = totstormfract(ig)*(mr_stormdust_mass(ig,l)-mr_dust_mass(ig,l))
415               zq_stormdust_number(ig,l) = totstormfract(ig)*(mr_stormdust_number(ig,l)-mr_dust_number(ig,l))
416             !! Particular case: there is less than initial dust mixing ratio after the vertical transport
417             ELSE
418               zq_dust_mass(ig,l) = (1.-totstormfract(ig))*mr_dust_mass(ig,l)+totstormfract(ig)*mr_stormdust_mass(ig,l)
419               zq_dust_number(ig,l) = (1.-totstormfract(ig))*mr_dust_number(ig,l)+totstormfract(ig)*mr_stormdust_number(ig,l)
420               zq_stormdust_mass(ig,l) = 0.
421               zq_stormdust_number(ig,l) = 0.
422             ENDIF
423             
424           ENDDO ! DO l=1,nlayer           
425         ENDIF ! IF storm(ig)
426        ENDDO ! DO ig=1,ngrid
[1974]427
[2079]428        !! **********************************************************************
[2201]429        !! 3.5 Calculation of the tendencies of the vertical transport
[2079]430        DO ig=1,ngrid
431         IF (storm(ig)) THEN
432           DO l=1,nlayer
433             dqvl_stormdust_mass(ig,l) = (zq_stormdust_mass(ig,l)- &
434                                   zq(ig,l,igcm_stormdust_mass)) /ptimestep
435             dqvl_stormdust_number(ig,l) = (zq_stormdust_number(ig,l)- &
436                                     zq(ig,l,igcm_stormdust_number)) /ptimestep
437             dqvl_dust_mass(ig,l) = (zq_dust_mass(ig,l)-zq(ig,l,igcm_dust_mass)) /ptimestep
438             dqvl_dust_number(ig,l) = (zq_dust_number(ig,l)-zq(ig,l,igcm_dust_number)) /ptimestep
439           ENDDO
440         ENDIF ! IF storm(ig)
441        ENDDO ! DO ig=1,ngrid           
[1974]442
443      ! **********************************************************************
[2079]444      ! 4. Detrainment: stormdust is converted to background dust
[1974]445      ! **********************************************************************
[2079]446        !! **********************************************************************
447        !! 4.1 Compute the coefficient of detrainmen
448        DO ig=1,ngrid
449          DO l=1,nlayer
450            IF ((max(abs(wrad(ig,l)),abs(wrad(ig,l+1))) .lt.  &
451                          wmin) .or. (zq_dust_mass(ig,l) .gt.  &
452                                     10000.*zq_stormdust_mass(ig,l))) THEN
[2160]453               detrain(ig,l)=1.
[2079]454            ELSE IF (max(abs(wrad(ig,l)),abs(wrad(ig,l+1)))   &
455                                                       .le. wmax) THEN
[2160]456               detrain(ig,l)=coeff_detrainment*                 &
457                             (((1-coefmin)/(wmin-wmax)**2)*     &
458                             (max(abs(wrad(ig,l)),abs(wrad(ig,l+1)))-wmax)**2 &
459                              +coefmin)
[2079]460            ELSE IF (max(abs(wrad(ig,l)),abs(wrad(ig,l+1))).gt. wmax ) THEN
[2160]461               detrain(ig,l)=coefmin
[2079]462            ELSE
[2160]463               detrain(ig,l)=coefmin
[2079]464            ENDIF
465          ENDDO ! DO l=1,nlayer
466        ENDDO ! DO ig=1,ngrid
467       
468        !! **********************************************************************
469        !! 4.2 Calculation of the tendencies of the detrainment
470        DO ig=1,ngrid
471          DO l=1,nlayer
[2160]472            dqdet_stormdust_mass(ig,l)=-detrain(ig,l)*zq_stormdust_mass(ig,l) &
[1974]473                                                        /ptimestep
[2160]474            dqdet_stormdust_number(ig,l)=-detrain(ig,l)*zq_stormdust_number(ig,l) &
[1974]475                                                        /ptimestep
[2079]476          ENDDO ! DO l=1,nlayer
477        ENDDO ! DO ig=1,ngrid
478           
479      ! **********************************************************************
480      ! 5. Final tendencies: vertical transport + detrainment
481      ! **********************************************************************
482        DO ig=1,ngrid
483          DO l=1,nlayer     
484          pdqrds(ig,l,igcm_stormdust_mass)=dqdet_stormdust_mass(ig,l) &
485                                                 +dqvl_stormdust_mass(ig,l)
486          pdqrds(ig,l,igcm_stormdust_number)=dqdet_stormdust_number(ig,l) &
487                                                 +dqvl_stormdust_number(ig,l)
488          pdqrds(ig,l,igcm_dust_mass)= -dqdet_stormdust_mass(ig,l) &
489                                       +dqvl_dust_mass(ig,l)
490          pdqrds(ig,l,igcm_dust_number)= -dqdet_stormdust_number(ig,l) &
491                                       +dqvl_dust_number(ig,l)
492          ENDDO ! DO l=1,nlayer
493        ENDDO ! DO ig=1,ngrid
[1974]494
[2201]495!      ! **********************************************************************
496!      ! 6. To prevent from negative values
497!      ! **********************************************************************
498!        DO ig=1,ngrid
499!          DO l=1,nlayer
500!            IF ((pq(ig,l,igcm_stormdust_mass) &
501!               +pdqrds(ig,l,igcm_stormdust_mass)*ptimestep .le. 0.) .or. &
502!              (pq(ig,l,igcm_stormdust_number) &
503!               +pdqrds(ig,l,igcm_stormdust_number)*ptimestep .le. 0.)) THEN
504!               pdqrds(ig,l,igcm_stormdust_mass)=-pq(ig,l,igcm_stormdust_mass)/ptimestep
505!               pdqrds(ig,l,igcm_stormdust_number)=-pq(ig,l,igcm_stormdust_number)/ptimestep
506!            ENDIF
507!          ENDDO ! nlayer
508!        ENDDO ! DO ig=1,ngrid
509!
510!        DO ig=1,ngrid
511!          DO l=1,nlayer           
512!            IF ((pq(ig,l,igcm_dust_mass) &
513!               +pdqrds(ig,l,igcm_dust_mass)*ptimestep .le. 0.) .or. &
514!              (pq(ig,l,igcm_dust_number) &
515!               +pdqrds(ig,l,igcm_dust_number)*ptimestep .le. 0.)) THEN
516!               pdqrds(ig,l,igcm_dust_mass)=-pq(ig,l,igcm_dust_mass)/ptimestep
517!               pdqrds(ig,l,igcm_dust_number)=-pq(ig,l,igcm_dust_number)/ptimestep
518!            ENDIF
519!          ENDDO ! nlayer
520!        ENDDO ! DO ig=1,ngrid
[2079]521       
522!=======================================================================
523! WRITEDIAGFI
524
[1974]525        call WRITEDIAGFI(ngrid,'lapserate','lapse rate in the storm', &
526           &                        'k/m',3,lapserate)
527        call WRITEDIAGFI(ngrid,'deltahr','extra heating rates', &
528           &                        'K/s',3,deltahr)
529        call writediagfi(ngrid,'scheme','which scheme',&
530                                                   ' ',2,real(scheme))
[2079]531
[1974]532        END SUBROUTINE rocketduststorm
533
[2201]534!======================================================================= 
535! **********************************************************************
536! Subroutine to determine the vertical transport with
537! a Van Leer advection scheme (copied from the sedimentation scheme --> see vlz_fi.F)
538!***********************************************************************
539      SUBROUTINE van_leer(nlay,q,pente_max,masse,w,wq)
[2079]540
[1974]541      IMPLICIT NONE
542
[2201]543!--------------------------------------------------------
544! Input/Output Variables
545!--------------------------------------------------------
546      INTEGER,INTENT(IN) :: nlay       ! number of atmospheric layers
547      REAL,INTENT(IN) ::  masse(nlay)  ! mass of the layer Dp/g
548      REAL,INTENT(IN) :: pente_max     != 2 conseillee
549      REAL,INTENT(INOUT) :: q(nlay)    ! mixing ratio (kg/kg)
550      REAL,INTENT(INOUT) :: w(nlay)    ! atmospheric mass "transferred" at each timestep (kg.m-2)
551      REAL,INTENT(INOUT) :: wq(nlay+1)
552
553!--------------------------------------------------------
554! Local Variables
555!--------------------------------------------------------
556
[1974]557      INTEGER i,l,j,ii
[2079]558      REAL dzq(nlay),dzqw(nlay),adzqw(nlay),dzqmax
559      REAL newmasse
560      REAL sigw, Mtot, MQtot
561      INTEGER m
[1974]562
[2201]563! **********************************************************************
564!  Mixing ratio vertical gradient at the levels
565! **********************************************************************
[1974]566      do l=2,nlay
567            dzqw(l)=q(l-1)-q(l)
568            adzqw(l)=abs(dzqw(l))
569      enddo
570
571      do l=2,nlay-1
572            if(dzqw(l)*dzqw(l+1).gt.0.) then
573                dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
574            else
575                dzq(l)=0.
576            endif
577            dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
578            dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
579      enddo
580
581      dzq(1)=0.
582      dzq(nlay)=0.
[2079]583
[2201]584! **********************************************************************
585!  Vertical advection
586! **********************************************************************
[1974]587
[2201]588       !! No flux at the model top:
[1974]589       wq(nlay+1)=0.
590
[2201]591       !! Surface flux up:
[2119]592       if(w(1).lt.0.) wq(1)=0. ! warning : not always valid
[1974]593
[2201]594       do l = 1,nlay-1
[2119]595
[2201]596!      1) Compute wq where w < 0 (up) (UPWARD TRANSPORT)
[2119]597!      ===============================
598
[1974]599         if(w(l+1).le.0)then
600!         Regular scheme (transfered mass < 1 layer)
601!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
602          if(-w(l+1).le.masse(l))then
603            sigw=w(l+1)/masse(l)
604            wq(l+1)=w(l+1)*(q(l)-0.5*(1.+sigw)*dzq(l))
[2201]605!!-------------------------------------------------------
606!          The following part should not be needed in the
607!          case of an integration over subtimesteps
608!!-------------------------------------------------------
[1974]609!         Extended scheme (transfered mass > 1 layer)
610!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
611          else
612             m = l-1
613             Mtot = masse(m+1)
614             MQtot = masse(m+1)*q(m+1)
615             if (m.le.0)goto 77
616             do while(-w(l+1).gt.(Mtot+masse(m)))
617!             do while(-w(l+1).gt.Mtot)
618                m=m-1
619                Mtot = Mtot + masse(m+1)
620                MQtot = MQtot + masse(m+1)*q(m+1)
621                if (m.le.0)goto 77
622             end do
623 77          continue
624
625             if (m.gt.0) then
626                sigw=(w(l+1)+Mtot)/masse(m)
[2079]627                wq(l+1)= -(MQtot + (-w(l+1)-Mtot)*         &
[1974]628                (q(m)-0.5*(1.+sigw)*dzq(m))  )
629             else
630                w(l+1) = -Mtot
631                wq(l+1) = -MQtot
632             end if
[2119]633          endif ! (-w(l+1).le.masse(l))
[1974]634     
[2201]635!      2) Compute wq where w > 0 (down) (DOWNWARD TRANSPORT)     
[1974]636!      ===============================
637
[2119]638         else if(w(l).gt.0.)then
[1974]639
640!         Regular scheme (transfered mass < 1 layer)
641!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
642          if(w(l).le.masse(l))then
643            sigw=w(l)/masse(l)
[2201]644            wq(l)=w(l)*(q(l)+0.5*(1.-sigw)*dzq(l))           
645!!-------------------------------------------------------
646!          The following part should not be needed in the
647!          case of an integration over subtimesteps
648!!-------------------------------------------------------
[1974]649!         Extended scheme (transfered mass > 1 layer)
650!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
651          else
652            m=l
653            Mtot = masse(m)
654            MQtot = masse(m)*q(m)
655            if(m.ge.nlay)goto 88
656            do while(w(l).gt.(Mtot+masse(m+1)))
657                m=m+1
658                Mtot = Mtot + masse(m)
659                MQtot = MQtot + masse(m)*q(m)
660                if(m.ge.nlay)goto 88
661            end do
662 88         continue
663            if (m.lt.nlay) then
664                sigw=(w(l)-Mtot)/masse(m+1)
665                wq(l)=(MQtot + (w(l)-Mtot)* &
666                          (q(m+1)+0.5*(1.-sigw)*dzq(m+1)) )
667            else
668                w(l) = Mtot
669                wq(l) = MQtot
670            end if
[2119]671          end if
[2100]672
[2119]673         end if ! w<0 (up)
[2100]674
[2119]675       enddo ! l = 1,nlay-1
[1974]676       
[2201]677       do l = 1,nlay
[1974]678
[2119]679!         it cannot entrain more than available mass !
680          if ( (wq(l+1)-wq(l)) .lt. -(masse(l)*q(l)) ) then
681            wq(l+1) = wq(l)-masse(l)*q(l)
682          end if
683
[1974]684          q(l)=q(l) +  (wq(l+1)-wq(l))/masse(l)
685
686       enddo
[2079]687       
[2201]688      END SUBROUTINE van_leer
[1974]689
[2079]690!=======================================================================
691! Initialization of the module variables
[1974]692
693       subroutine ini_rocketduststorm_mod(ngrid)
694       
695       implicit none
696       
697       integer, intent(in) :: ngrid
698       
699       allocate(dustliftday(ngrid))
700       
701       end subroutine ini_rocketduststorm_mod
702       
703       subroutine end_rocketduststorm_mod
704       
705       implicit none
706       
707       if (allocated(dustliftday)) deallocate(dustliftday)
708
709       end subroutine end_rocketduststorm_mod       
710     
711      END MODULE rocketduststorm_mod
Note: See TracBrowser for help on using the repository browser.