Changeset 2201


Ignore:
Timestamp:
Dec 18, 2019, 4:00:47 PM (5 years ago)
Author:
mvals
Message:

Mars GCM
Update of the rocketduststorm parametrization: clean-up of the code + back to the sub-timesteps method for the Van Leer transport, homogeneous with the one used in topmons_mod.F.
MV

Location:
trunk/LMDZ.MARS
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r2199 r2201  
    28002800== 13/12/2019 == MV
    28012801Implementation of a new parametrization of the dust entrainment by slope winds above the sub-grid scale topography. The parametrization is activated with the flag slpwind=.true. (set to "false" by default) in callphys.def. The new parametrization involves the new tracers topdust_mass and topdust_number.
     2802
     2803== 18/12/2019 == MV
     2804Update of the rocketduststorm parametrization: clean-up of the code + back to the sub-timesteps method for the Van Leer transport, homogeneous with the one used in topmons_mod.F.

  • trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90

    r2199 r2201  
    1111!=======================================================================
    1212! calculation of the vertical flux
    13 ! call of vl_storm : Van Leer transport scheme of the dust tracers
     13! call of van_leer : Van Leer transport scheme of the dust tracers
    1414! detrainement of stormdust in the background dust
    1515! -----------------------------------------------------------------------
    … …  
    9797! Local variables
    9898!--------------------------------------------------------
    99       INTEGER l,ig,tsub,iq,ll
     99      INTEGER l,ig,iq,ll
    100100!     local variables from callradite.F       
    101101      REAL zdtlw1(ngrid,nlayer)    ! (K/s) storm
    … …  
    117117      REAL mr_dust_mass(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (dust mass)
    118118      REAL mr_dust_number(ngrid,nlayer) ! intermediate mixing ratio to calculate van leer transport with the "real" concentration (sdust number)
    119 
    120       REAL zq_vl_col(nlayer) ! column intermediate tracer used by Van Leer (number)
    121       REAL zn_vl_col(nlayer) ! column intermediate tracer used by Van Leer (mass)
    122119                   
    123120      REAL dqvl_stormdust_mass(ngrid,nlayer)    ! tendancy of vertical transport (stormdust mass)
    … …  
    140137           
    141138      REAL,PARAMETER:: coefmin =0.025 ! 0<coefmin<1 Minimum fraction of stormdust detrained 
    142       REAL,PARAMETER:: wmin =0.1!0.25 ! stormdust detrainment if wrad < wmin 
     139      REAL,PARAMETER:: wmin =0.25 ! stormdust detrainment if wrad < wmin 
    143140      REAL,PARAMETER:: wmax =10.   ! maximum vertical velocity of the rocket dust storms (m/s)
    144      
     141
     142!     subtimestep
     143      INTEGER tsub
     144      INTEGER nsubtimestep    !number of subtimestep when calling van_leer
     145      REAL subtimestep        !ptimestep/nsubtimestep
     146      REAL dtmax              !considered time needed for dust to cross one layer
     147      REAL,PARAMETER:: secu=3.!3.      !coefficient on wspeed to avoid dust crossing many layers during subtimestep
     148
    145149!     diagnostics
    146150      REAL lapserate(ngrid,nlayer)
    … …  
    265269
    266270            scheme(ig)=2
    267             !! This is the env. lapse rate
    268             zdtvert(ig,1)=0.
    269             DO l=1,nlayer-1
    270               zdtvert(ig,l+1)=(zt(ig,l+1)-zt(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l))
    271               lapserate(ig,l+1)=zdtvert(ig,l+1)
    272             ENDDO
    273271     
    274272            !! computing heating rates gradient at boundraies of each layer
    … …  
    301299                           zt(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /      &
    302300                              (pzlay(ig,l+1)-pzlay(ig,l))
    303             ENDDO ! DO l=1,nlayer-1   
     301            ENDDO ! DO l=1,nlayer-1
     302
     303            !! This is the env. lapse rate
     304            zdtvert(ig,1)=0.
     305            DO l=1,nlayer-1
     306              zdtvert(ig,l+1)=(ztlev(ig,l+1)-ztlev(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l))
     307              lapserate(ig,l+1)=zdtvert(ig,l+1)
     308            ENDDO
    304309
    305310            !! **********************************************************************
    … …  
    338343        ENDDO ! DO ig=1,ngrid
    339344
    340         !! **********************************************************************
    341         !! 3.2 Vertical transport by a Van Leer scheme
    342345        DO ig=1,ngrid
    343346          IF (storm(ig)) THEN
    344            
     347        !! **********************************************************************
     348        !! 3.2 Compute the subtimestep to conserve the mass in the Van Leer transport
     349            dtmax=ptimestep
     350            DO l=2,nlayer
     351              IF (wrad(ig,l).lt.0.) THEN
     352                 dtmax=min(dtmax,(pzlev(ig,l)-pzlev(ig,l-1))/  &
     353                                   (secu*abs(wrad(ig,l))))
     354              ELSE IF (wrad(ig,l).gt.0.) then
     355                 dtmax=min(dtmax,(pzlev(ig,l+1)-pzlev(ig,l))/  &
     356                                   (secu*abs(wrad(ig,l))))
     357              ENDIF
     358            ENDDO
     359            nsubtimestep= int(ptimestep/dtmax)
     360            subtimestep=ptimestep/float(nsubtimestep)
     361            !! Mass flux generated by wup in kg/m2
     362            DO l=1,nlayer
     363               w(ig,l)=wrad(ig,l)*pplev(ig,l)/(r*ztlev(ig,l)) &
     364                        *subtimestep
     365            ENDDO ! l=1,nlayer
     366
     367        !! **********************************************************************
     368        !! 3.3 Vertical transport by a Van Leer scheme
    345369            !! Mass of atmosphere in the layer
    346370            DO l=1,nlayer
    347371               masse_col(l)=(pplev(ig,l)-pplev(ig,l+1))/g
    348372            ENDDO
    349            
    350             !! Mass flux in kg/m2
    351             DO l=1,nlayer
    352                w(ig,l)=wrad(ig,l)*(pplev(ig,l)/(r*ztlev(ig,l)))*ptimestep
    353             ENDDO
    354            
    355             !! Vector column
    356             DO l=1,nlayer
    357                zq_vl_col(l)= mr_stormdust_mass(ig,l)
    358                zn_vl_col(l)= mr_stormdust_number(ig,l)
    359             ENDDO
    360            
    361             !! Van Leer scheme
    362             CALL vl_storm(nlayer,zq_vl_col,2.,   &
    363                   masse_col,w(ig,:),wqmass(ig,:))
    364             CALL vl_storm(nlayer,zn_vl_col,2.,  &
    365                   masse_col,w(ig,:),wqnumber(ig,:))
    366             !! Mass mixing ratio after transport     
    367             mr_stormdust_mass(ig,:) = zq_vl_col(:)
    368             mr_stormdust_number(ig,:) = zn_vl_col(:)
    369                        
     373            !! Mass flux in kg/m2 if you are not using the subtimestep
     374            !DO l=1,nlayer
     375            !   w(ig,l)=wrad(ig,l)*(pplev(ig,l)/(r*ztlev(ig,l)))*ptimestep
     376            !ENDDO
     377            !! Loop over the subtimestep
     378            DO tsub=1,nsubtimestep
     379              !! Van Leer scheme
     380              wqmass(ig,:)=0.
     381              wqnumber(ig,:)=0.
     382              CALL van_leer(nlayer,mr_stormdust_mass(ig,:),2.,   &
     383                    masse_col,w(ig,:),wqmass(ig,:))
     384              CALL van_leer(nlayer,mr_stormdust_number(ig,:),2.,  &
     385                    masse_col,w(ig,:),wqnumber(ig,:))
     386            ENDDO !tsub=...           
     387             
    370388          ENDIF ! IF storm(ig)
    371389        ENDDO ! DO ig=1,ngrid 
    372390
    373391        !! **********************************************************************
    374         !! 3.3 Re-calculation of the mixing ratios after vertical transport
     392        !! 3.4 Re-calculation of the mixing ratios after vertical transport
    375393        DO ig=1,ngrid
    376394         IF (storm(ig)) THEN
    … …  
    396414
    397415        !! **********************************************************************
    398         !! 3.4 Calculation of the tendencies of the vertical transport
     416        !! 3.5 Calculation of the tendencies of the vertical transport
    399417        DO ig=1,ngrid
    400418         IF (storm(ig)) THEN
    … …  
    462480        ENDDO ! DO ig=1,ngrid
    463481
    464       ! **********************************************************************
    465       ! 6. To prevent from negative values
    466       ! **********************************************************************
    467         DO ig=1,ngrid
    468           DO l=1,nlayer
    469             IF ((pq(ig,l,igcm_stormdust_mass) &
    470                +pdqrds(ig,l,igcm_stormdust_mass)*ptimestep .le. 0.) .or. &
    471               (pq(ig,l,igcm_stormdust_number) &
    472                +pdqrds(ig,l,igcm_stormdust_number)*ptimestep .le. 0.)) THEN
    473                pdqrds(ig,l,igcm_stormdust_mass)=-pq(ig,l,igcm_stormdust_mass)/ptimestep
    474                pdqrds(ig,l,igcm_stormdust_number)=-pq(ig,l,igcm_stormdust_number)/ptimestep
    475             ENDIF
    476           ENDDO ! nlayer
    477         ENDDO ! DO ig=1,ngrid
    478 
    479         DO ig=1,ngrid
    480           DO l=1,nlayer           
    481             IF ((pq(ig,l,igcm_dust_mass) &
    482                +pdqrds(ig,l,igcm_dust_mass)*ptimestep .le. 0.) .or. &
    483               (pq(ig,l,igcm_dust_number) &
    484                +pdqrds(ig,l,igcm_dust_number)*ptimestep .le. 0.)) THEN
    485                pdqrds(ig,l,igcm_dust_mass)=-pq(ig,l,igcm_dust_mass)/ptimestep
    486                pdqrds(ig,l,igcm_dust_number)=-pq(ig,l,igcm_dust_number)/ptimestep
    487             ENDIF
    488           ENDDO ! nlayer
    489         ENDDO ! DO ig=1,ngrid
     482!      ! **********************************************************************
     483!      ! 6. To prevent from negative values
     484!      ! **********************************************************************
     485!        DO ig=1,ngrid
     486!          DO l=1,nlayer
     487!            IF ((pq(ig,l,igcm_stormdust_mass) &
     488!               +pdqrds(ig,l,igcm_stormdust_mass)*ptimestep .le. 0.) .or. &
     489!              (pq(ig,l,igcm_stormdust_number) &
     490!               +pdqrds(ig,l,igcm_stormdust_number)*ptimestep .le. 0.)) THEN
     491!               pdqrds(ig,l,igcm_stormdust_mass)=-pq(ig,l,igcm_stormdust_mass)/ptimestep
     492!               pdqrds(ig,l,igcm_stormdust_number)=-pq(ig,l,igcm_stormdust_number)/ptimestep
     493!            ENDIF
     494!          ENDDO ! nlayer
     495!        ENDDO ! DO ig=1,ngrid
     496!
     497!        DO ig=1,ngrid
     498!          DO l=1,nlayer           
     499!            IF ((pq(ig,l,igcm_dust_mass) &
     500!               +pdqrds(ig,l,igcm_dust_mass)*ptimestep .le. 0.) .or. &
     501!              (pq(ig,l,igcm_dust_number) &
     502!               +pdqrds(ig,l,igcm_dust_number)*ptimestep .le. 0.)) THEN
     503!               pdqrds(ig,l,igcm_dust_mass)=-pq(ig,l,igcm_dust_mass)/ptimestep
     504!               pdqrds(ig,l,igcm_dust_number)=-pq(ig,l,igcm_dust_number)/ptimestep
     505!            ENDIF
     506!          ENDDO ! nlayer
     507!        ENDDO ! DO ig=1,ngrid
    490508       
    491509!=======================================================================
    … …  
    501519        END SUBROUTINE rocketduststorm
    502520
    503 !=======================================================================
    504 ! VAN LEER
    505 
    506       SUBROUTINE vl_storm(nlay,q,pente_max,masse,w,wq)
    507 !
    508 !     Auteurs:   P.Le Van, F.Hourdin, F.Forget
    509 !
    510 !    ********************************************************************
    511 !     Shema  d'advection " pseudo amont " dans la verticale
    512 !    pour appel dans la physique (sedimentation)
    513 !    ********************************************************************
    514 !    q rapport de melange (kg/kg)...
    515 !    masse : masse de la couche Dp/g
    516 !    w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
    517 !    pente_max = 2 conseillee
    518 !
    519 !
    520 !   --------------------------------------------------------------------
     521!======================================================================= 
     522! **********************************************************************
     523! Subroutine to determine the vertical transport with
     524! a Van Leer advection scheme (copied from the sedimentation scheme --> see vlz_fi.F)
     525!***********************************************************************
     526      SUBROUTINE van_leer(nlay,q,pente_max,masse,w,wq)
     527
    521528      IMPLICIT NONE
    522 !
    523 
    524 !   Arguments:
    525 !   ----------
    526       INTEGER,INTENT(IN) :: nlay ! number of atmospheric layers
    527       REAL masse(nlay),pente_max
    528       REAL q(nlay)
    529       REAL w(nlay)
    530       REAL wq(nlay+1)
    531 !
    532 !      Local
    533 !   ---------
    534 !
     529
     530!--------------------------------------------------------
     531! Input/Output Variables
     532!--------------------------------------------------------
     533      INTEGER,INTENT(IN) :: nlay       ! number of atmospheric layers
     534      REAL,INTENT(IN) ::  masse(nlay)  ! mass of the layer Dp/g
     535      REAL,INTENT(IN) :: pente_max     != 2 conseillee
     536      REAL,INTENT(INOUT) :: q(nlay)    ! mixing ratio (kg/kg)
     537      REAL,INTENT(INOUT) :: w(nlay)    ! atmospheric mass "transferred" at each timestep (kg.m-2)
     538      REAL,INTENT(INOUT) :: wq(nlay+1)
     539
     540!--------------------------------------------------------
     541! Local Variables
     542!--------------------------------------------------------
     543
    535544      INTEGER i,l,j,ii
    536 !
    537 
    538545      REAL dzq(nlay),dzqw(nlay),adzqw(nlay),dzqmax
    539546      REAL newmasse
    … …  
    541548      INTEGER m
    542549
    543       REAL      SSUM,CVMGP,CVMGT
    544       INTEGER ismax,ismin
    545 
    546 
    547 !    On oriente tout dans le sens de la pression c'est a dire dans le
    548 !    sens de W
    549 
     550! **********************************************************************
     551!  Mixing ratio vertical gradient at the levels
     552! **********************************************************************
    550553      do l=2,nlay
    551554            dzqw(l)=q(l-1)-q(l)
    … …  
    554557
    555558      do l=2,nlay-1
    556 #ifdef CRAY
    557             dzq(l)=0.5*
    558      ,      cvmgp(dzqw(l)+dzqw(l+1),0.,dzqw(l)*dzqw(l+1))
    559 #else
    560559            if(dzqw(l)*dzqw(l+1).gt.0.) then
    561560                dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
    … …  
    563562                dzq(l)=0.
    564563            endif
    565 #endif
    566564            dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
    567565            dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
    … …  
    571569      dzq(nlay)=0.
    572570
    573 ! ---------------------------------------------------------------
    574 !   .... calcul des termes d'advection verticale  .......
    575 ! ---------------------------------------------------------------
    576 
    577 ! calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
    578 !
    579 !      No flux at the model top:
     571! **********************************************************************
     572!  Vertical advection
     573! **********************************************************************
     574
     575       !! No flux at the model top:
    580576       wq(nlay+1)=0.
    581577
    582 !      Surface flux up:
     578       !! Surface flux up:
    583579       if(w(1).lt.0.) wq(1)=0. ! warning : not always valid
    584580
    585        do l = 1,nlay-1  ! loop different than when w>0
    586 
    587 !      1) Compute wq where w < 0 (up)
     581       do l = 1,nlay-1
     582
     583!      1) Compute wq where w < 0 (up) (UPWARD TRANSPORT)
    588584!      ===============================
    589585
    590586         if(w(l+1).le.0)then
    591 
    592587!         Regular scheme (transfered mass < 1 layer)
    593588!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    595590            sigw=w(l+1)/masse(l)
    596591            wq(l+1)=w(l+1)*(q(l)-0.5*(1.+sigw)*dzq(l))
     592!!-------------------------------------------------------
     593!          The following part should not be needed in the
     594!          case of an integration over subtimesteps
     595!!-------------------------------------------------------
    597596!         Extended scheme (transfered mass > 1 layer)
    598597!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    621620          endif ! (-w(l+1).le.masse(l))
    622621     
    623 !      2) Compute wq where w > 0 (down)    
     622!      2) Compute wq where w > 0 (down) (DOWNWARD TRANSPORT)     
    624623!      ===============================
    625624
    … …  
    630629          if(w(l).le.masse(l))then
    631630            sigw=w(l)/masse(l)
    632             wq(l)=w(l)*(q(l)+0.5*(1.-sigw)*dzq(l))
    633            
    634 
     631            wq(l)=w(l)*(q(l)+0.5*(1.-sigw)*dzq(l))           
     632!!-------------------------------------------------------
     633!          The following part should not be needed in the
     634!          case of an integration over subtimesteps
     635!!-------------------------------------------------------
    635636!         Extended scheme (transfered mass > 1 layer)
    636637!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    661662       enddo ! l = 1,nlay-1
    662663       
    663        do l = 1,nlay          ! loop different than when w<0
     664       do l = 1,nlay
    664665
    665666!         it cannot entrain more than available mass !
    … …  
    672673       enddo
    673674       
    674       END SUBROUTINE vl_storm
     675      END SUBROUTINE van_leer
    675676
    676677!=======================================================================
Note: See TracChangeset for help on using the changeset viewer.