Changeset 1974 for trunk/LMDZ.MARS/libf


Ignore:
Timestamp:
Jul 18, 2018, 4:48:34 PM (6 years ago)
Author:
mvals
Message:

Mars GCM:
Integration of the detached dust layer parametrizations (rocket dust storm, slope wind lifting, CW, and dust injection scheme, DB).
Still experimental, default behaviour (rdstorm=.false., dustinjection=0) identical to previous revision.
NB: Updated newstart requires an updated "surface.nc" containing the "hmons" field.
EM+MV

Location:
trunk/LMDZ.MARS/libf
Files:
5 added
20 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/datareadnc.F

    r1918 r1974  
    11c=======================================================================
    22      SUBROUTINE datareadnc(relief,phisinit,alb,ith,z0,
    3      &                    zmea,zstd,zsig,zgam,zthe)
     3     &                    zmea,zstd,zsig,zgam,zthe,
     4     &                    hmons,summit,zavg)
    45c=======================================================================
    56c
    … …  
    4243c=======================================================================
    4344
    44       use ioipsl_getincom, only: getin
    45       USE comconst_mod, ONLY: g,pi
     45! to use  'getin'
     46      use ioipsl_getincom, only: getin
     47      use comconst_mod, only: g,pi
    4648      use datafile_mod, only: datadir
     49      use avg_horiz_mod, only: avg_horiz
     50      use mvc_horiz_mod, only: mvc_horiz
    4751
    4852      implicit none
    … …  
    7478      REAL,intent(out) :: zgam(imdp1*jmdp1)
    7579      REAL,intent(out) :: zthe(imdp1*jmdp1)
    76 
     80      REAL,intent(out) :: hmons(imdp1*jmdp1) !CW17,361*180 hmons
     81      REAL,intent(out) :: summit(imdp1*jmdp1) !CW17,361*180 hmons
     82      REAL,intent(out) :: zavg(imdp1*jmdp1) !CW17,361*180 hmons
     83     
    7784! Local variables:
    7885      REAL        zdata(imd*jmdp1)
    … …  
    95102
    96103      CHARACTER*20 string
    97       DIMENSION string(0:4)
     104      DIMENSION string(0:6)
    98105
    99106
    … …  
    233240          call grid_noro1(360, 180, longitude, latitude, zdata,
    234241     .         iim, jjp1, rlonv, rlatu, zmea,zstd,zsig,zgam,zthe)
     242     
     243          !CW17
     244          call avg_horiz(imd,jmdp1,iim,jjm,longitude,
     245     .                latitude,rlonv,rlatu,zdata,pfield)
     246
     247          do j=1,jjp1 !CW17 49, iimp1=65, the last column = first column
     248             pfield(iimp1*j) =  pfield(1+iimp1*(j-1))
     249          enddo
     250          do i=1,iimp1*jjp1
     251c             if (pfield(i) .ne. -999999.) then
     252                zavg(i) = pfield(i)
     253c             else
     254c               zavg(i)=-999999.
     255c             endif
     256          enddo
    235257
    236258      endif
    … …  
    297319c    FIN
    298320c-----------------------------------------------------------------------
    299 
     321 
     322c=======================================================================
     323c<<<<<<<<<<<<<<<<<<<<<<<add new dataset>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
     324      ! name of dataset
     325      string(5) = 'hmons' !subgrid hmons
     326c     the following data could be useful in future, but currently,
     327c     we don't need to put them into restartfi.nc
     328      string(6) = 'summit' !subgrid summit
     329c     string(7) = 'base'   !base of summit (not subgrid)
     330c     string(8) = 'bottom' !subgrid base
     331      do k=5,6
     332          write(*,*) 'string',k,string(k) 
     333c-----------------------------------------------------------------------
     334c    Lecture NetCDF 
     335c-----------------------------------------------------------------------
     336
     337          ierr = NF_INQ_VARID (unit, string(k), nvarid)
     338          if (ierr.ne.nf_noerr) then
     339            write(*,*) 'datareadnc error, cannot find ',trim(string(k))
     340            write(*,*) ' in file ',trim(datadir),'/surface.nc'
     341            stop
     342          endif
     343
     344c-----------------------------------------------------------------------
     345c    initialisation
     346c-----------------------------------------------------------------------
     347          call initial0(iimp1*jjp1,pfield)
     348          call initial0(imd*jmdp1,zdata)
     349          call initial0(imdp1*jmdp1,zdataS)
     350
     351#ifdef NC_DOUBLE
     352          ierr = NF_GET_VAR_DOUBLE(unit, nvarid, zdata)
     353#else
     354          ierr = NF_GET_VAR_REAL(unit, nvarid, zdata)
     355#endif
     356          if (ierr.ne.nf_noerr) then
     357            write(*,*) 'datareadnc: error failed loading ',
     358     .                  trim(string(k))
     359            stop
     360          endif
     361
     362c-----------------------------------------------------------------------
     363c   Passage de zdata en grille (imdp1 jmdp1)
     364c-----------------------------------------------------------------------
     365          do j=1,jmdp1     !  180
     366              do i=1,imd   !  360
     367                  !copy zdata to zdataS, line by line
     368                  !      i+ 361 *(j-1)          i+  360*(j-1)
     369                  zdataS(i+imdp1*(j-1)) = zdata(i+ngridmxgdat*(j-1))
     370              enddo
     371              ! the last column = the first column of zdata
     372              zdataS(imdp1+imdp1*(j-1)) = zdata(1+ngridmxgdat*(j-1))
     373          enddo
     374
     375          if (k .eq. 5 .or. k .eq. 6) then
     376              ! hmons, summit, keep the maximum of subgrids
     377              call mvc_horiz(imd,jmdp1,iim,jjm,longitude,
     378     .                latitude,rlonv,rlatu,zdata,pfield)
     379          endif
     380
     381
     382c-----------------------------------------------------------------------
     383c    Periodicite   
     384c-----------------------------------------------------------------------
     385
     386          do j=1,jjp1 ! 49, iimp1=65, the last column = first column
     387             pfield(iimp1*j) =  pfield(1+iimp1*(j-1))
     388          enddo
     389 
     390c-----------------------------------------------------------------------
     391c    Sauvegarde des champs   
     392c-----------------------------------------------------------------------
     393
     394          if (k.eq.5) then ! hmons
     395             do i=1,iimp1*jjp1
     396                  if (pfield(i) .ne. -999999.) then
     397                    hmons(i) = pfield(i)
     398                  else
     399                    hmons(i)=0.
     400                  endif
     401             enddo
     402          endif
     403         
     404          if (k.eq.6) then  ! summit
     405             do i=1,iimp1*jjp1
     406                    summit(i) = pfield(i)
     407             enddo
     408          endif
     409
     410      enddo
     411c<<<<<<<<<<<<<<<<<<<<<<<<<done add new dataset>>>>>>>>>>>>>>>>>>>>>>>>>>
     412c=======================================================================
    300413      END

  • trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F

    r1944 r1974  
    2727      use surfdat_h, only: phisfi, z0, zmea, zstd, zsig, zgam, zthe,
    2828     &                     albedodat, z0_default, qsurf, tsurf,
    29      &                     co2ice, emis
     29     &                     co2ice, emis, hmons
    3030      use comsoil_h, only: inertiedat, layer, mlayer, nsoilmx, tsoil
    3131      use control_mod, only: day_step, iphysiq, anneeref, planet_type
    … …  
    9797      real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
    9898      real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
     99      real hmonsS(iip1,jjp1)
     100      real summitS(iip1,jjp1)
     101      real zavgS(iip1,jjp1)
    99102      real z0S(iip1,jjp1)
    100103
    … …  
    380383
    381384      CALL datareadnc(relief,phis,alb,surfith,z0S,
    382      &          zmeaS,zstdS,zsigS,zgamS,ztheS)
     385     &          zmeaS,zstdS,zsigS,zgamS,ztheS,
     386     &          hmonsS,summitS,zavgS)
    383387
    384388      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
    … …  
    392396      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zgamS,zgam)
    393397      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,ztheS,zthe)
     398      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,hmonsS,hmons)
     399!      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,summitS,summit)
     400!      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zavgS,zavg)
    394401
    395402      endif ! of if (choix_1.ne.1)
    … …  
    16421649      call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
    16431650     .              nqtot,dtphys,real(day_ini),0.0,
    1644      .              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
     1651     .              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe,
     1652     .              hmons)
    16451653      call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
    16461654     &              dtphys,hour_ini,

  • trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F

    r1861 r1974  
    77      SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
    88     &    pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
    9      &    nueffrad,QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d,
     9     &    QREFvis3d,QREFir3d,omegaREFir3d,
     10     &    totstormfract,clearatm,dsords,
    1011     &    clearsky,totcloudfrac)
    11                                                    
     12                                                        
    1213! to use  'getin'
    1314      USE ioipsl_getincom, only: getin
    1415      use tracer_mod, only: noms, igcm_h2o_ice, igcm_dust_mass,
    1516     &                      igcm_dust_submicron, rho_dust, rho_ice,
    16      &                      nqdust
     17     &                      nqdust, igcm_stormdust_mass
    1718      use geometry_mod, only: latitude ! grid point latitudes (rad)
    1819      use comgeomfi_h, only: sinlat ! sines of grid point latitudes
    … …  
    2627     &            iaerdust,tauvis,
    2728     &            iaer_dust_conrath,iaer_dust_doubleq,
    28      &            iaer_dust_submicron,iaer_h2o_ice
     29     &            iaer_dust_submicron,iaer_h2o_ice,
     30     &            iaer_stormdust_doubleq
     31      USE calcstormfract_mod
    2932       IMPLICIT NONE
    3033c=======================================================================
    … …  
    5659c   QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
    5760c   QREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths;
    58 c   omegaREFvis3d(ngrid,nlayer,naerkind) \ 3d single scat. albedo
    5961c   omegaREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths;
    6062c
    … …  
    6567c   aerosol      aerosol(ig,l,1) is the dust optical
    6668c                depth in layer l, grid point ig
    67 
     69c   taualldust   CW17 total opacity for all dust scatterer stormdust included
    6870c
    6971c=======================================================================
    70 #include "callkeys.h"
     72      include "callkeys.h"
    7173
    7274c-----------------------------------------------------------------------
    … …  
    7779c    Input/Output
    7880c    ------------
    79       INTEGER ngrid,nlayer,nq
    80 
    81       REAL ls,zday,expfactor   
    82       REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
    83       REAL pq(ngrid,nlayer,nq)
    84       REAL tauref(ngrid), tau(ngrid,naerkind)
    85       REAL aerosol(ngrid,nlayer,naerkind)
    86       REAL dsodust(ngrid,nlayer)
    87       REAL reffrad(ngrid,nlayer,naerkind)
    88       REAL nueffrad(ngrid,nlayer,naerkind)
    89       REAL QREFvis3d(ngrid,nlayer,naerkind)
    90       REAL QREFir3d(ngrid,nlayer,naerkind)
    91       REAL omegaREFvis3d(ngrid,nlayer,naerkind)
    92       REAL omegaREFir3d(ngrid,nlayer,naerkind)
    93       real,intent(in) :: totcloudfrac(ngrid) ! total cloud fraction
    94       logical,intent(in) :: clearsky ! true for part without clouds,
    95                                      ! false for part with clouds (total or sub-grid clouds)
    96       real :: CLFtot ! total cloud fraction
     81      INTEGER, INTENT(IN) :: ngrid,nlayer,nq
     82      REAL, INTENT(IN) :: ls,zday
     83      REAL, INTENT(IN) ::  pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
     84      REAL, INTENT(IN) ::  pq(ngrid,nlayer,nq)
     85      REAL, INTENT(OUT) :: tauref(ngrid)
     86      REAL, INTENT(OUT) :: tau(ngrid,naerkind)
     87      REAL, INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind)
     88      REAL, INTENT(OUT) :: dsodust(ngrid,nlayer)
     89      REAL, INTENT(OUT) ::  dsords(ngrid,nlayer) !dso of stormdust 
     90      REAL, INTENT(INOUT) :: reffrad(ngrid,nlayer,naerkind)
     91      REAL, INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind)
     92      REAL, INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
     93      REAL, INTENT(IN) ::  omegaREFir3d(ngrid,nlayer,naerkind)
     94      LOGICAL, INTENT(IN) :: clearatm
     95      REAL, INTENT(IN) :: totstormfract(ngrid)
     96      REAL, INTENT(OUT) ::  tauscaling(ngrid) ! Scaling factor for qdust and Ndust
     97      REAL,INTENT(IN) :: totcloudfrac(ngrid) ! total cloud fraction
     98      LOGICAL,INTENT(IN) :: clearsky ! true for part without clouds,false for part with clouds (total or sub-grid clouds)
    9799c
    98100c    Local variables :
    99101c    -----------------
     102      REAL CLFtot ! total cloud fraction
     103      real expfactor
    100104      INTEGER l,ig,iq,i,j
    101105      INTEGER iaer           ! Aerosol index
    … …  
    110114c       for dust or water ice particles in the radiative transfer
    111115c       (see callradite.F for more information).
    112       REAL taudusttmp(ngrid)! Temporary dust opacity
    113                                !   used before scaling
    114       REAL tauscaling(ngrid) ! Scaling factor for qdust and Ndust
     116      REAL taudusttmp(ngrid)! Temporary dust opacity used before scaling
     117      REAL taubackdusttmp(ngrid)! Temporary background dust opacity used before scaling
     118      REAL taualldust(ngrid)! dust opacity all dust
     119      REAL taudust(ngrid)! dust opacity dust doubleq
     120      REAL taustormdust(ngrid)! dust opacity stormdust doubleq
     121      REAL taustormdusttmp(ngrid)! dust opacity stormdust doubleq before tauscaling
    115122      REAL taudustvis(ngrid) ! Dust opacity after scaling
    116123      REAL taudusttes(ngrid) ! Dust opacity at IR ref. wav. as
    … …  
    157164
    158165      tau(1:ngrid,1:naerkind)=0
     166      dsords(:,:)=0. !CW17: initialize dsords
     167      dsodust(:,:)=0.
    159168
    160169! identify tracers
    … …  
    166175        DO iaer=1,naerkind
    167176          txt=name_iaer(iaer)
    168           IF (txt(1:4).eq."dust") THEN
     177        ! CW17: choice tauscaling for stormdust or not
     178          IF ((txt(1:4).eq."dust").OR.(txt(1:5).eq."storm")) THEN
    169179            naerdust=naerdust+1
    170180            iaerdust(naerdust)=iaer
    … …  
    213223        call getin("tauvis",tauvis)
    214224
    215         IF (freedust) THEN
     225        IF (freedust.or.rdstorm) THEN ! if rdstorm no need to held opacity constant at the first levels
    216226          cstdustlevel = 1
    217227        ELSE
    218           cstdustlevel = cstdustlevel0
     228          cstdustlevel = cstdustlevel0 !Opacity in the first levels is held constant to
     229                                       !avoid unrealistic values due to constant lifting
    219230        ENDIF
    220231
    … …  
    374385     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
    375386                ! DENSITY SCALED OPACITY IN INFRARED:
    376                 dsodust(ig,l) =
    377      &          (  0.75 * QREFir3d(ig,cstdustlevel,iaer) /
    378      &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
    379      &          pq(ig,cstdustlevel,igcm_dust_mass)
     387!                dsodust(ig,l) =
     388!     &          (  0.75 * QREFir3d(ig,cstdustlevel,iaer) /
     389!     &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
     390!     &          pq(ig,cstdustlevel,igcm_dust_mass)
    380391              ENDDO
    381392            ELSE
    … …  
    387398     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
    388399                ! DENSITY SCALED OPACITY IN INFRARED:
    389                 dsodust(ig,l) =
    390      &          (  0.75 * QREFir3d(ig,l,iaer) /
    391      &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
    392      &          pq(ig,l,igcm_dust_mass)
     400!                dsodust(ig,l) =
     401!     &          (  0.75 * QREFir3d(ig,l,iaer) /
     402!     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
     403!     &          pq(ig,l,igcm_dust_mass)
    393404              ENDDO
    394405            ENDIF
    … …  
    473484        ENDIF ! end (clearsky)
    474485
    475 c       3. Outputs -- Now done in physiq.F
    476 !        IF (ngrid.NE.1) THEN
    477 !          CALL WRITEDIAGFI(ngrid,'tauVIS','tauext VIS refwvl',
    478 !     &      ' ',2,taucloudvis)
    479 !          CALL WRITEDIAGFI(ngrid,'tauTES','tauabs IR refwvl',
    480 !     &      ' ',2,taucloudtes)
    481 !          IF (callstats) THEN
    482 !            CALL wstats(ngrid,'tauVIS','tauext VIS refwvl',
    483 !     &        ' ',2,taucloudvis)
    484 !            CALL wstats(ngrid,'tauTES','tauabs IR refwvl',
    485 !     &        ' ',2,taucloudtes)
    486 !          ENDIF
    487 !        ELSE
    488 !         CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
    489 !        ENDIF
     486c==================================================================
     487        CASE("stormdust_doubleq") aerkind ! CW17 : Two-moment scheme for
     488c       stormdust  (transport of mass and number mixing ratio)
     489c==================================================================
     490c       aerosol is calculated twice : once within the dust storm (clearatm=false)
     491c       and once in the part of the mesh without dust storm (clearatm=true)
     492        aerosol(1:ngrid,1:nlayer,iaer) = 0.
     493        IF (clearatm) THEN  ! considering part of the mesh without storm
     494          aerosol(1:ngrid,1:nlayer,iaer)=1.e-25
     495        ELSE  ! part of the mesh with concentred dust storm
     496          DO l=1,nlayer
     497             IF (l.LE.cstdustlevel) THEN
     498c          Opacity in the first levels is held constant to
     499c           avoid unrealistic values due to constant lifting:
     500               DO ig=1,ngrid
     501                  aerosol(ig,l,iaer) =
     502     &           (  0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
     503     &           ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
     504     &           pq(ig,cstdustlevel,igcm_stormdust_mass) *
     505     &           ( pplev(ig,l) - pplev(ig,l+1) ) / g
     506               ENDDO
     507             ELSE
     508              DO ig=1,ngrid
     509                 aerosol(ig,l,iaer) =
     510     &           (  0.75 * QREFvis3d(ig,l,iaer) /
     511     &           ( rho_dust * reffrad(ig,l,iaer) )  ) *
     512     &           pq(ig,l,igcm_stormdust_mass) *
     513     &           ( pplev(ig,l) - pplev(ig,l+1) ) / g
     514              ENDDO
     515             ENDIF
     516          ENDDO
     517        ENDIF
    490518c==================================================================
    491519        END SELECT aerkind
    … …  
    498526c   is originally scaled to an equivalent odpref Pa pressure surface.
    499527c -----------------------------------------------------------------
     528
    500529
    501530#ifdef DUSTSTORM
    … …  
    640669      IF (freedust) THEN
    641670          tauscaling(:) = 1.
    642 
     671c        opacity obtained with stormdust
     672        IF (rdstorm) THEN
     673           taustormdusttmp(1:ngrid)=0.
     674           DO l=1,nlayer
     675             DO ig=1,ngrid
     676                taustormdusttmp(ig) = taustormdusttmp(ig)+
     677     &            aerosol(ig,l,iaerdust(2))
     678             ENDDO
     679           ENDDO
     680           !opacity obtained with background dust only
     681           taubackdusttmp(1:ngrid)=0. 
     682           DO l=1,nlayer
     683             DO ig=1,ngrid
     684                taubackdusttmp(ig) = taubackdusttmp(ig)+
     685     &           aerosol(ig,l,iaerdust(1))
     686             ENDDO
     687           ENDDO
     688        ENDIF !rdsstorm
    643689      ELSE
    644690c       Temporary scaling factor
    … …  
    669715     &                aerosol(ig,l,iaerdust(iaer))* tauscaling(ig))
    670716          ENDDO
    671         ENDDO
    672       ENDDO
    673 
    674       DO l=1,nlayer
    675         DO ig=1,ngrid
    676           dsodust(ig,l) = max(1E-20,dsodust(ig,l)* tauscaling(ig))
    677717        ENDDO
    678718      ENDDO
    … …  
    700740        tauref(:) = tauref(:) * odpref / pplev(:,1)
    701741      ENDIF
    702      
    703 c output for debug
    704 c        IF (ngrid.NE.1) THEN
    705 c             CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
    706 c     &      '#',2,taudusttmp)
    707 c             CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
    708 c     &      '#',2,tauscaling)
    709 c        ELSE
    710 c             CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
    711 c     &      '#',0,taudusttmp)
    712 c             CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
    713 c     &      '#',0,tauscaling)
    714 c        ENDIF
     742
    715743c -----------------------------------------------------------------
    716744c Column integrated visible optical depth in each point
    … …  
    724752      ENDDO
    725753
     754c     for diagnostics: opacity for all dust scatterers stormdust included
     755      taualldust(1:ngrid)=0.
     756      DO iaer=1,naerdust
     757        DO l=1,nlayer
     758          DO ig=1,ngrid
     759            taualldust(ig) = taualldust(ig) +
     760     &                         aerosol(ig,l,iaerdust(iaer))
     761          ENDDO
     762        ENDDO
     763      ENDDO
     764     
     765      IF (rdstorm) THEN
     766 
     767c     for diagnostics: opacity for dust in background only 
     768       taudust(1:ngrid)=0.
     769        DO l=1,nlayer
     770         DO ig=1,ngrid
     771           taudust(ig) = taudust(ig) +
     772     &                       aerosol(ig,l,iaer_dust_doubleq)
     773         ENDDO
     774        ENDDO
     775
     776c     for diagnostics: opacity for dust in storm only 
     777       taustormdust(1:ngrid)=0.
     778        DO l=1,nlayer
     779         DO ig=1,ngrid
     780           taustormdust(ig) = taustormdust(ig) +
     781     &                       aerosol(ig,l,iaer_stormdust_doubleq)
     782         ENDDO
     783        ENDDO
     784 
     785      ENDIF
     786     
     787
    726788#ifdef DUSTSTORM
    727789      IF (firstcall) THEN
    … …  
    733795c Density scaled opacity and column opacity output
    734796c -----------------------------------------------------------------
    735 c      dsodust(1:ngrid,1:nlayer) = 0.
    736 c      DO iaer=1,naerdust
    737 c        DO l=1,nlayer
    738 c          DO ig=1,ngrid
    739 c            dsodust(ig,l) = dsodust(ig,l) +
    740 c     &                      aerosol(ig,l,iaerdust(iaer)) * g /
    741 c     &                      (pplev(ig,l) - pplev(ig,l+1))
    742 c          ENDDO
    743 c        ENDDO
    744 c        IF (ngrid.NE.1) THEN
    745 c          write(txt2,'(i1.1)') iaer
    746 c          call WRITEDIAGFI(ngrid,'taudust'//txt2,
    747 c     &                    'Dust col opacity',
    748 c     &                    ' ',2,tau(1,iaerdust(iaer)))
    749 c          IF (callstats) THEN
    750 c            CALL wstats(ngrid,'taudust'//txt2,
    751 c     &                 'Dust col opacity',
    752 c     &                 ' ',2,tau(1,iaerdust(iaer)))
    753 c          ENDIF
    754 c        ENDIF
    755 c      ENDDO
    756 
    757 c      IF (ngrid.NE.1) THEN
    758 c       CALL WRITEDIAGFI(ngrid,'dsodust','tau*g/dp',
    759 c    &                    'm2.kg-1',3,dsodust)
    760 c        IF (callstats) THEN
    761 c          CALL wstats(ngrid,'dsodust',
    762 c     &               'tau*g/dp',
    763 c     &               'm2.kg-1',3,dsodust)
    764 c        ENDIF
    765 c      ELSE
    766 c        CALL WRITEDIAGFI(ngrid,"dsodust","dsodust","m2.kg-1",1,
    767 c     &                   dsodust)
    768 c      ENDIF ! of IF (ngrid.NE.1)
    769 c -----------------------------------------------------------------
     797          IF (rdstorm) then
     798            DO l=1,nlayer
     799              IF (l.LE.cstdustlevel) THEN
     800                DO ig=1,ngrid
     801                  dsodust(ig,l)=dsodust(ig,l) +
     802     &                      aerosol(ig,l,iaer_dust_doubleq) * g /
     803     &                      (pplev(ig,l) - pplev(ig,l+1))
     804
     805                  dsords(ig,l) = dsords(ig,l) +
     806     &              aerosol(ig,l,iaer_stormdust_doubleq)* g/
     807     &              (pplev(ig,l) - pplev(ig,l+1))
     808                ENDDO
     809              ELSE
     810                DO ig=1,ngrid
     811                  dsodust(ig,l) =dsodust(ig,l) +
     812     &                      aerosol(ig,l,iaer_dust_doubleq) * g /
     813     &                      (pplev(ig,l) - pplev(ig,l+1))
     814                   dsords(ig,l) = dsords(ig,l) +
     815     &                        aerosol(ig,l,iaer_stormdust_doubleq)* g/
     816     &                        (pplev(ig,l) - pplev(ig,l+1))
     817                ENDDO
     818              ENDIF
     819            ENDDO
     820          ENDIF
     821
     822c -----------------------------------------------------------------
     823c -----------------------------------------------------------------
     824c aerosol/X for stormdust to prepare calculation of radiative transfer
     825c -----------------------------------------------------------------
     826      if (rdstorm) then
     827         DO l=1,nlayer
     828           DO ig=1,ngrid
     829               aerosol(ig,l,iaer_stormdust_doubleq) =
     830     &           aerosol(ig,l,iaer_stormdust_doubleq)/totstormfract(ig)
     831           ENDDO
     832         ENDDO
     833      endif
     834
    770835
    771836      END SUBROUTINE aeropacity

  • trunk/LMDZ.MARS/libf/phymars/callkeys.h

    r1818 r1974  
    1313     &   ,lifting,freedust,callddevil,scavenging,sedimentation          &
    1414     &   ,activice,water,tifeedback,microphys,supersat,caps,photochem   &
    15      &   ,calltherm,callrichsl,callslope,tituscap,callyamada4,co2clouds,&
    16      &   co2useh2o,meteo_flux,CLFvaryingCO2,spantCO2,CLFvarying,        &
    17      &   satindexco2
     15     &   ,calltherm,callrichsl,callslope,tituscap,callyamada4,co2clouds &
     16     &   ,co2useh2o,meteo_flux,CLFvaryingCO2,spantCO2,CLFvarying        &
     17     &   ,satindexco2,rdstorm
    1818     
    1919      COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche    &
    20      &   ,dustbin,nltemodel,nircorr,solvarmod,solvaryear
     20     &   ,dustbin,nltemodel,nircorr,solvarmod,solvaryear,dustinjection
    2121     
    2222      COMMON/callkeys_r/topdustref,semi,alphan,euveff,                  &
    23      &   tke_heat_flux,dustrefir,fixed_euv_value,CLFfixval
     23     &   tke_heat_flux,dustrefir,fixed_euv_value,CLFfixval,             &
     24     &   coeff_injection
    2425     
    2526      LOGICAL callrad,calldifv,calladj,callcond,callsoil,               &
    … …  
    4243      real euveff
    4344      real tke_heat_flux
     45      real coeff_injection ! dust injection scheme coefficient
    4446      real CLFfixval
    4547
    … …  
    5355      integer solvarmod   ! model for solar EUV variation
    5456      integer solvaryear  ! mars year for realisticly varying solar EUV
     57      integer dustinjection ! dust injection scheme number
    5558
    5659      logical rayleigh
    5760      logical tracer
    5861      integer dustbin
    59       logical freedust
     62      logical freedust  
    6063      logical active,doubleq,submicron,lifting,callddevil,scavenging
     64      logical rdstorm ! rocket dust storm parametrization
    6165      logical sedimentation
    6266      logical activice,tifeedback,supersat,caps

  • trunk/LMDZ.MARS/libf/phymars/callradite_mod.F

    r1969 r1974  
    77      SUBROUTINE callradite(icount,ngrid,nlayer,nq,zday,ls,pq,albedo,
    88     $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
    9      $     dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
    10      &     tauref,tau,aerosol,dsodust,tauscaling,taucloudtes,rdust,rice,
    11      &     nuice,co2ice,clearsky,totcloudfrac)
     9     $     dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,
     10     $     fluxtop_sw,tauref,tau,aerosol,dsodust,tauscaling,
     11     $     taucloudtes,rdust,rice,nuice,co2ice,rstormdust,
     12     $     totstormfract,clearatm,dsords,
     13     $     clearsky,totcloudfrac)
    1214
    1315      use aeropacity_mod, only: aeropacity
    … …  
    1618      use dimradmars_mod, only: naerkind, name_iaer,
    1719     &            iaer_dust_conrath,iaer_dust_doubleq,
    18      &            iaer_dust_submicron,iaer_h2o_ice
     20     &            iaer_dust_submicron,iaer_h2o_ice,
     21     &            iaer_stormdust_doubleq
    1922      use yomlw_h, only: gcp, nlaylte
    2023      use comcstfi_h, only: g,cpp
    … …  
    154157c                         the "naerkind" kind of aerosol optical
    155158c                         properties.
    156 
    157159c=======================================================================
    158160c
    … …  
    160162c    -------------
    161163c
    162 #include "callkeys.h"
     164      include "callkeys.h"
    163165
    164166c-----------------------------------------------------------------------
    … …  
    170172
    171173      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
    172       REAL,INTENT(IN) :: tauscaling(ngrid) ! Conversion factor for
     174      REAL,INTENT(INOUT) :: tauscaling(ngrid) ! Conversion factor for
    173175                               ! qdust and Ndust
    174176      REAL,INTENT(IN) :: albedo(ngrid,2),emis(ngrid)
    … …  
    194196      REAL,INTENT(OUT) :: nuice(ngrid,nlayer)  ! Estimated effective variance
    195197      REAL,INTENT(IN) :: co2ice(ngrid)           ! co2 ice surface layer (kg.m-2)
     198
     199c     rocket dust storm
     200      LOGICAL,INTENT(IN) :: clearatm ! true for background dust
     201      REAL,INTENT(IN) :: totstormfract(ngrid) ! dust storm mesh fraction
     202      REAL,INTENT(OUT) :: rstormdust(ngrid,nlayer)  ! Storm dust geometric mean radius (m)
     203      REAL dsords(ngrid,nlayer) ! density scaled opacity for rocket dust storm dust
     204
    196205c     sub-grid scale water ice clouds
    197       real,intent(out) :: totcloudfrac(ngrid)
    198       logical,intent(in) :: clearsky
     206      LOGICAL,INTENT(IN) :: clearsky
     207      REAL,INTENT(IN) :: totcloudfrac(ngrid)
    199208
    200209c
    … …  
    292301         iaer_dust_submicron=0
    293302         iaer_h2o_ice=0
     303         iaer_stormdust_doubleq=0
    294304
    295305         aer_count=0
    … …  
    326336           enddo
    327337         endif
     338         if (rdstorm.AND.active) then
     339           do iaer=1,naerkind
     340             if (name_iaer(iaer).eq."stormdust_doubleq") then
     341               iaer_stormdust_doubleq = iaer
     342               aer_count = aer_count + 1
     343             endif
     344           enddo
     345         end if
    328346
    329347c        Check that we identified all tracers:
    … …  
    372390c     Updating aerosol size distributions:
    373391      CALL updatereffrad(ngrid,nlayer,
    374      &                rdust,rice,nuice,
     392     &                rdust,rstormdust,rice,nuice,
    375393     &                reffrad,nueffrad,
    376394     &                pq,tauscaling,tau,pplay)
    … …  
    385403c     Computing aerosol optical depth in each layer:
    386404      CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
    387      &     pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
    388      &     nueffrad,QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d,
    389      &     clearsky,totcloudfrac)
     405     &    pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
     406     &    QREFvis3d,QREFir3d,omegaREFir3d,
     407     &    totstormfract,clearatm,dsords,
     408     &    clearsky,totcloudfrac)
    390409
    391410c     Starting loop on sub-domain

  • trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F

    r1962 r1974  
    44
    55      CONTAINS
    6      
    7       SUBROUTINE callsedim(ngrid,nlay, ptimestep,
    8      &                pplev,zlev, zlay, pt, pdt, rdust, rice,
     6
     7      SUBROUTINE callsedim(ngrid,nlay,ptimestep,
     8     &                pplev,zlev,zlay,pt,pdt,rdust,rstormdust,rice,
    99     &                rsedcloud,rhocloud,
    10      &                pq, pdqfi, pdqsed,pdqs_sed,nq,
     10     &                pq,pdqfi,pdqsed,pdqs_sed,nq,
    1111     &                tau,tauscaling)
    12 ! to use  'getin'
     12
    1313      USE ioipsl_getincom, only: getin
    1414      USE updaterad, only: updaterdust,updaterice_micro,updaterice_typ
    … …  
    1717     &                      igcm_ccn_mass, igcm_ccn_number,
    1818     &                      igcm_h2o_ice, nuice_sed, nuice_ref,
    19      &                      igcm_ccnco2_mass, igcm_ccnco2_number,
    20      &                      igcm_co2_ice
     19     &                      igcm_ccnco2_mass,igcm_ccnco2_number,
     20     &                      igcm_co2_ice, igcm_stormdust_mass,
     21     &                      igcm_stormdust_number
    2122      USE newsedim_mod, ONLY: newsedim
    2223      USE comcstfi_h, ONLY: g
    … …  
    6162c    Aerosol radius provided by the water ice microphysical scheme:
    6263      real,intent(out) :: rdust(ngrid,nlay) ! Dust geometric mean radius (m)
     64      real,intent(out) :: rstormdust(ngrid,nlay) ! Stormdust geometric mean radius (m)
    6365      real,intent(out) :: rice(ngrid,nlay)  ! H2O Ice geometric mean radius (m)
    6466c     Sedimentation radius of water ice
    … …  
    8991      real r0dust(ngrid,nlay) ! geometric mean radius used for
    9092                                    !   dust (m)
    91 !      real r0ccn(ngrid,nlay)  ! geometric mean radius used for
     93      real r0stormdust(ngrid,nlay) ! Geometric mean radius used for stormdust (m)
    9294!                                    !   CCNs (m)
    9395      real,save :: beta ! correction for the shape of the ice particles (cf. newsedim)
    94 
    9596c     for ice radius computation
    9697      REAL Mo,No
    … …  
    133134      INTEGER,SAVE :: iccn_number ! index of tracer containing CCN number
    134135                                  !   mix. ratio
     136      INTEGER,SAVE :: istormdust_mass  !  index of tracer containing
     137                                       !stormdust mass mix. ratio
     138      INTEGER,SAVE :: istormdust_number !  index of tracer containing
     139                                        !stormdust number mix. ratio                     
    135140      INTEGER,SAVE :: iccnco2_number ! index of tracer containing CCN number
    136141      INTEGER,SAVE :: iccnco2_mass ! index of tracer containing CCN number
    … …  
    140145      LOGICAL,SAVE :: firstcall=.true.
    141146
     147
     148
    142149c    ** un petit test de coherence
    143150c       --------------------------
    144 
    145151      ! AS: firstcall OK absolute
    146152      IF (firstcall) THEN
    … …  
    240246       ENDIF                    !of if (co2clouds)
    241247
    242 
    243         IF (water) THEN
     248       IF (water) THEN
    244249         write(*,*) "correction for the shape of the ice particles ?"
    245250         beta=0.75 ! default value
    … …  
    251256            write(*,*) "water_param nueff Radiative:", nuice_ref
    252257          ENDIF
    253         ENDIF
     258       ENDIF
     259
     260       IF (rdstorm) THEN ! identifying stormdust tracers for sedimentation
     261           istormdust_mass=0      ! dummy initialization
     262           istormdust_number=0    ! dummy initialization
     263
     264           do iq=1,nq
     265             if (noms(iq).eq."stormdust_mass") then
     266               istormdust_mass=iq
     267               write(*,*)"callsedim: istormdust_mass=",istormdust_mass
     268             endif
     269             if (noms(iq).eq."stormdust_number") then
     270               istormdust_number=iq
     271               write(*,*)"callsedim: istormdust_number=",
     272     &                                           istormdust_number
     273             endif
     274           enddo
     275
     276           ! check that we did find the tracers
     277           if ((istormdust_mass.eq.0).or.(istormdust_number.eq.0)) then
     278             write(*,*) 'callsedim: error! could not identify'
     279             write(*,*) ' tracers for stormdust mass and number mixing'
     280             write(*,*) ' ratio and rdstorm is activated!'
     281             stop
     282           endif
     283       ENDIF !of if (rdstorm)
    254284     
    255285        firstcall=.false.
    … …  
    268298      zt(1:ngrid,1:nlay)=pt(1:ngrid,1:nlay)
    269299     &                         +pdt(1:ngrid,1:nlay)*ptimestep
    270 
    271300
    272301c    Computing the different layer properties
    … …  
    295324        end do
    296325      endif
    297 
    298 
     326      ! rocket dust storm
     327      if (rdstorm) then
     328        do l=1,nlay
     329          do ig=1, ngrid
     330     
     331         call updaterdust(zqi(ig,l,igcm_stormdust_mass),
     332     &               zqi(ig,l,igcm_stormdust_number),r0stormdust(ig,l),
     333     &               tauscaling(ig))
     334         
     335          end do
     336        end do
     337      endif
    299338c =================================================================
    300339      do iq=1,nq
    … …  
    307346c -----------------------------------------------------------------
    308347          if ((doubleq.and.
    309      &        ((iq.eq.idust_mass).or.
    310      &         (iq.eq.idust_number)))) then
     348     &     ((iq.eq.idust_mass).or.(iq.eq.idust_number).or.
     349     &     (iq.eq.istormdust_mass).or.(iq.eq.istormdust_number)))) then
    311350     
    312351c           Computing size distribution:
    313352c           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    314353
    315 c            if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
     354            if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
    316355              do  l=1,nlay
    317356                do ig=1, ngrid
    … …  
    319358                end do
    320359              end do
    321               sigma0 = varian
     360            else if ((iq.eq.istormdust_mass).or.
     361     &                                (iq.eq.istormdust_number)) then
     362              do  l=1,nlay
     363                do ig=1, ngrid
     364                  r0(ig,l)=r0stormdust(ig,l)
     365                end do
     366              end do
     367            endif
     368            sigma0 = varian
    322369
    323370c        Computing mass mixing ratio for each particle size
    324371c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    325           IF ((iq.EQ.idust_mass).or.(iq.EQ.iccn_mass)) then
     372          IF ((iq.EQ.idust_mass).or.(iq.EQ.istormdust_mass)) then
    326373            radpower = 2
    327374          ELSE  ! number
    … …  
    376423
    377424          do ir=1,nr
    378          
    379425               call newsedim(ngrid,nlay,1,1,ptimestep,
    380426     &         pplev,masse,epaisseur,zt,rd(ir),(/rho_dust/),qr(1,1,ir),
    … …  
    391437                 zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir)
    392438               ENDDO
    393              ENDDO
     439             ENDDO           
    394440          enddo ! of do ir=1,nr
    395441c -----------------------------------------------------------------
    … …  
    457503       ENDDO
    458504      endif ! of if (doubleq)
    459      
     505
     506      if (rdstorm) then
     507       DO l = 1, nlay
     508        DO ig=1,ngrid
     509         call updaterdust(zqi(ig,l,igcm_stormdust_mass),
     510     &                zqi(ig,l,igcm_stormdust_number),rstormdust(ig,l),
     511     &                tauscaling(ig))   
     512        ENDDO
     513       ENDDO
     514      endif ! of if (rdstorm)
     515 
    460516c     Update the ice particle size "rice"
    461517c     -------------------------------------
    … …  
    490546     
    491547      END MODULE callsedim_mod
     548

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

    r1770 r1974  
    99  real,save,allocatable :: mu0(:)
    1010  real,save,allocatable :: fract(:)
     11  real,save,allocatable :: local_time(:) ! local solar time as fraction of day (0,1)
    1112
    1213contains
    … …  
    1920    allocate(mu0(ngrid))
    2021    allocate(fract(ngrid))
     22    allocate(local_time(ngrid))
    2123 
    2224  end subroutine ini_comsaison_h
    … …  
    2931    if (allocated(mu0)) deallocate(mu0)
    3032    if (allocated(fract)) deallocate(fract)
     33    if (allocated(local_time)) deallocate(local_time)
    3134
    3235  end subroutine end_comsaison_h

  • trunk/LMDZ.MARS/libf/phymars/conf_phys.F

    r1918 r1974  
    280280         write(*,*)" calllott = ",calllott
    281281
     282! rocket dust storm injection scheme
     283         write(*,*)"call rocket dust storm and slope lifting",
     284     &              " parametrization"
     285         rdstorm=.false. ! default value
     286         call getin("rdstorm",rdstorm)
     287         write(*,*)" rdstorm = ",rdstorm
    282288
    283289         write(*,*)"rad.transfer is computed every iradia",
    … …  
    359365         endif
    360366
     367! dust injection scheme
     368        dustinjection=0 ! default: no injection scheme
     369        call getin("dustinjection",dustinjection)
     370        write(*,*)" dustinjection = ",dustinjection
     371! dust injection scheme coefficient       
     372        coeff_injection=1. ! default value
     373        call getin("coeff_injection",coeff_injection)
     374        write(*,*)" coeff_in,jection = ",coeff_injection
     375
    361376! free evolving dust
    362377! freedust=true just says that there is no lifting and no dust opacity scaling.
    … …  
    372387         ! this test is valid in GCM case
    373388         ! ... not in mesoscale case, for we want to activate mesoscale lifting
    374          if (freedust.and.lifting) then
    375            print*,'if freedust is used, then lifting should not be used'
    376            print*,'lifting forced to false !!'
    377            lifting=.false.
     389         if (freedust.and.dustinjection.eq.0)then
     390           if(lifting) then
     391             print*,'if freedust is used and dustinjection = 0,
     392     &      then lifting should not be used'
     393             stop
     394           endif
    378395         endif
    379396#endif
    380 
     397         if (dustinjection.eq.1)then
     398           if(.not.lifting) then
     399             print*,"if dustinjection=1, then lifting should be true"
     400             stop
     401           endif
     402           if(.not.freedust) then
     403             print*,"if dustinjection=1, then freedust should be true"
     404             stop
     405           endif
     406         endif
     407! rocket dust storm
     408! Test of incompatibility:
     409! if rdstorm is used, then doubleq should be true
     410         if (rdstorm.and..not.doubleq) then
     411           print*,'if rdstorm is used, then doubleq should be used !'
     412           stop
     413         endif
     414         if (rdstorm.and..not.active) then
     415           print*,'if rdstorm is used, then active should be used !'
     416           stop
     417         endif
     418         if (rdstorm.and..not.lifting) then
     419           print*,'if rdstorm is used, then lifting should be used !'
     420           stop
     421         endif
     422         if (rdstorm.and..not.freedust) then
     423           print*,'if rdstorm is used, then freedust should be used !'
     424           stop
     425         endif
     426         if (rdstorm.and.(dustinjection.eq.0)) then
     427           print*,'if rdstorm is used, then dustinjection should
     428     &             be used !'
     429           stop
     430         endif
    381431! Dust IR opacity
    382432         write(*,*)" Wavelength for infrared opacity of dust ?"
    … …  
    642692           WRITE(*,*) 'equal to 2.'
    643693           CALL ABORT
    644          ELSE IF ( (.NOT.activice).AND.(naerkind.GT.1) ) THEN
    645            WRITE(*,*) 'naerkind is greater than unity, but'
    646            WRITE(*,*) 'activice has not been set to .true.'
    647            WRITE(*,*) 'in callphys.def; this is not logical!'
    648            CALL ABORT
    649694         ENDIF
    650695
    … …  
    666711          ! and picky compilers who know name_iaer(2) is out of bounds
    667712          j=2
     713         IF (rdstorm.AND..NOT.activice) name_iaer(2) =
     714     &                       "stormdust_doubleq" !! storm dust two-moment scheme
     715         IF (rdstorm.AND.water.AND.activice) name_iaer(3) =
     716     &                                             "stormdust_doubleq"
    668717         IF (water.AND.activice) name_iaer(j) = "h2o_ice"      !! radiatively-active clouds
    669718         IF (submicron.AND.active) name_iaer(j) = "dust_submicron" !! JBM experimental stuff
    670719         endif ! of if (nq.gt.1)
     720
    671721c        ----------------------------------------------------------
    672722

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

    r1771 r1974  
    3232                              ! submicron population of dust
    3333                              ! particles
     34  integer iaer_stormdust_doubleq ! Storm dust profile is given by the
     35                              ! mass mixing ratio of the two moment scheme
     36                              ! method (doubleq)
    3437  integer iaer_h2o_ice ! Water ice particles
    3538

  • trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F

    r1944 r1974  
    1111     &                     albedice, iceradius, dtemisice, z0,
    1212     &                     zmea, zstd, zsig, zgam, zthe, phisfi,
    13      &                     watercaptag
     13     &                     watercaptag, hmons
    1414      use slope_mod, only: theta_sl, psi_sl
    1515      use phyredem, only: physdem0,physdem1
    … …  
    557557      zgam(1)=0.E+0
    558558      zthe(1)=0.E+0
     559c  for rocket dust storm scheme     
     560      hmons(1)=0.E+0
    559561
    560562
    … …  
    716718      call physdem0("startfi.nc",longitude,latitude,nsoilmx,ngrid,llm,
    717719     &              nq,dtphys,float(day0),time,cell_area,
    718      &              albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
     720     &              albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,hmons)
    719721      call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,
    720722     &              dtphys,time,

  • trunk/LMDZ.MARS/libf/phymars/initracer.F

    r1720 r1974  
    2424
    2525
    26 #include "callkeys.h"
     26      include "callkeys.h"
    2727
    2828      integer,intent(in) :: ngrid ! number of atmospheric columns
    … …  
    3232      integer iq,ig,count
    3333      real r0_lift , reff_lift, nueff_lift
     34      real r0_storm,reff_storm
    3435c     Ratio of small over large dust particles (used when both
    3536c       doubleq and the submicron mode are active); In Montmessin
    3637c       et al. (2002), a value of 25 has been deduced;
    3738      real, parameter :: popratio = 25.
    38       character(len=20) :: txt ! to store some text
     39      character(len=30) :: txt ! to store some text
    3940
    4041c-----------------------------------------------------------------------
    … …  
    6970      igcm_h2o_vap=0
    7071      igcm_h2o_ice=0
     72      igcm_stormdust_mass=0
     73      igcm_stormdust_number=0
    7174      igcm_co2=0
    7275      igcm_co=0
    … …  
    147150        enddo
    148151      endif ! of if (submicron)
     152       if (rdstorm) then
     153        do iq=1,nq
     154          if (noms(iq).eq."stormdust_mass") then
     155            igcm_stormdust_mass=iq
     156            count=count+1
     157          endif
     158          if (noms(iq).eq."stormdust_number") then
     159            igcm_stormdust_number=iq
     160            count=count+1
     161          endif
     162        enddo
     163      endif ! of if (rdstorm)     
    149164      ! 2. find chemistry and water tracers
    150165      do iq=1,nq
    … …  
    422437#else
    423438        alpha_lift(igcm_dust_mass)=1.e-6 !1.e-6 !Lifted mass coeff
     439        IF (dustinjection.ge.1) THEN
     440                reff_lift = 3.0e-6 ! Effective radius of lifted dust (m)
     441                alpha_lift(igcm_dust_mass)=(4/3.)*reff_lift*rho_dust
     442     &                                          /2.4
     443        ENDIF
    424444#endif
    425445
    … …  
    437457        write(*,*) "initracer: doubleq_param alpha_lift:",
    438458     &    alpha_lift(igcm_dust_mass)
     459!c ----------------------------------------------------------------------
     460!c rocket dust storm scheme
     461!c lifting tracer stormdust using same distribution than
     462!c normal dust
     463        if (rdstorm) then
     464          reff_storm=3.e-6 ! reff_lift !3.e-6
     465          r0_storm=reff_storm/ref_r0
     466          rho_q(igcm_stormdust_mass)=rho_dust
     467          rho_q(igcm_stormdust_number)=rho_dust
     468
     469          alpha_devil(igcm_stormdust_mass)=9.e-9   ! dust devil lift mass coeff
     470          alpha_lift(igcm_stormdust_mass)=4./3./2.4*reff_storm*rho_dust
     471
     472          write(*,*) 'alpha_lift(rds):',alpha_lift(igcm_stormdust_mass)
     473 
     474          alpha_devil(igcm_stormdust_number)=r3n_q*
     475     &                      alpha_devil(igcm_stormdust_mass)/r0_storm**3
     476          alpha_lift(igcm_stormdust_number)=r3n_q*
     477     &                       alpha_lift(igcm_stormdust_mass)/r0_storm**3
     478 
     479          radius(igcm_stormdust_mass) = reff_storm
     480          radius(igcm_stormdust_number) = reff_storm
     481        end if !(rdstorm)
     482!c ----------------------------------------------------------------------
     483     
    439484      else
    440485
    … …  
    618663          endif
    619664       endif
    620        
     665 
     666       if (rdstorm) then
     667       ! verify that we indeed have stormdust_mass and stormdust_number tracers
     668         if (igcm_stormdust_mass.eq.0) then
     669           write(*,*) "initracer: error !!"
     670           write(*,*) "  cannot use rdstorm option without ",
     671     &                "a stormdust_mass tracer !"
     672           stop
     673         endif
     674         if (igcm_stormdust_number.eq.0) then
     675           write(*,*) "initracer: error !!"
     676           write(*,*) "  cannot use rdstorm option without ",
     677     &                "a stormdust_number tracer !"
     678           stop
     679         endif
     680       endif
     681     
    621682       if (callnlte) then ! NLTE requirements
    622683         if (nltemodel.ge.1) then

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

    r1944 r1974  
    1212  use tracer_mod, only: noms ! tracer names
    1313  use surfdat_h, only: phisfi, albedodat, z0, z0_default,&
    14                        zmea, zstd, zsig, zgam, zthe
     14                       zmea, zstd, zsig, zgam, zthe, hmons
    1515  use iostart, only: nid_start, open_startphy, close_startphy, &
    1616                     get_field, get_var, inquire_field, &
    … …  
    1919
    2020  implicit none
     21 
     22  include "callkeys.h"
    2123!======================================================================
    2224! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
    … …  
    179181endif
    180182
     183! hmons
     184call get_field("hmons",hmons,found)
     185if (.not.found) then
     186  write(*,*) "WARNING: phyetat0: Failed loading <hmons>"
     187  if (rdstorm) then
     188    call abort
     189  else
     190    write(*,*) "will continue anyway..."
     191    write(*,*) "because you may not need it."
     192    hmons(:)=0.
     193  end if
     194else
     195  do ig=1,ngrid
     196    if (hmons(ig) .eq. -999999.)  hmons(ig)=0.
     197  enddo
     198  write(*,*) "phyetat0: <hmons> range:", &
     199             minval(hmons), maxval(hmons)
     200endif
    181201     
    182202! Time axis
    … …  
    277297if (.not.found) then
    278298  write(*,*) "phyetat0: <tauscaling> not in file"
    279   tauscaling(:) = -1
     299  tauscaling(:) = 1
    280300else
    281301  write(*,*) "phyetat0: dust conversion factor <tauscaling> range:", &
    … …  
    296316call get_field("wstar",wstar,found,indextime)
    297317if (.not.found) then
    298   write(*,*) "phyetat0: <totcloudfrac> not in file! Set to zero"
     318  write(*,*) "phyetat0: <wstar> not in file! Set to zero"
    299319  wstar(:)=0
    300320else

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

    r1944 r1974  
    1 !
    2 ! $Id: $
    3 !
    41module phyredem
    52
    … …  
    107subroutine physdem0(filename,lonfi,latfi,nsoil,ngrid,nlay,nq, &
    118                         phystep,day_ini,time,airefi, &
    12                          alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe)
     9                         alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe, &
     10                         phmons)
    1311! create physics restart file and write time-independent variables
    1412  use comsoil_h, only: inertiedat, volcapa, mlayer
    … …  
    1614  use surfdat_h, only: zmea, zstd, zsig, zgam, zthe, &
    1715                       z0_default, albedice, emisice, emissiv, &
    18                        iceradius, dtemisice, phisfi, z0
     16                       iceradius, dtemisice, phisfi, z0,   &
     17                       hmons
    1918  use dimradmars_mod, only: tauvis
    2019  use iostart, only : open_restartphy, close_restartphy, &
    … …  
    4645  real,intent(in) :: pzgam(ngrid)
    4746  real,intent(in) :: pzthe(ngrid)
    48  
     47  real,intent(in) :: phmons(ngrid)
     48
    4949  real :: tab_cntrl(length) ! nb "length=100" defined in iostart module
    5050 
    … …  
    131131  call put_field("ZGAM","Relief: gamma parameter",zgam)
    132132  call put_field("ZTHE","Relief: theta parameter",zthe)
    133  
     133  ! Write hmons if it was read in startfi
     134  if (maxval(hmons).ne.minval(hmons)) then
     135    call put_field("hmons","Relief: hmons parameter",hmons)
     136  end if
     137     
    134138  ! Soil thermal inertia
    135139  call put_field("inertiedat","Soil thermal inertia",ith)

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

    r1922 r1974  
    4747      use time_phylmdz_mod, only: init_time
    4848      use co2cloud_mod, only: ini_co2cloud,end_co2cloud
     49      use rocketduststorm_mod, only: ini_rocketduststorm_mod, &
     50                                     end_rocketduststorm_mod
    4951
    5052      IMPLICIT NONE
    … …  
    112114      call ini_co2cloud(ngrid,nlayer)
    113115     
     116      ! allocate arrays in "rocketduststorm_mod":
     117      call end_rocketduststorm_mod
     118      call ini_rocketduststorm_mod(ngrid)
     119     
    114120      END SUBROUTINE phys_state_var_init
    115121

  • trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

    r1973 r1974  
    1717      use co2cloud_mod, only: co2cloud, mem_Mccn_co2, mem_Mh2o_co2,
    1818     &                        mem_Nccn_co2
    19       use aeropacity_mod
    20       use callradite_mod
    21       use callsedim_mod
     19      use callradite_mod, only: callradite
     20      use callsedim_mod, only: callsedim
     21      use rocketduststorm_mod, only: rocketduststorm, dustliftday
     22      use calcstormfract_mod, only: calcstormfract
    2223      use tracer_mod, only: noms, mmol, igcm_co2, igcm_n2, igcm_co2_ice,
    2324     &                      igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice,
    … …  
    2728     &                      igcm_dust_mass, igcm_dust_number, igcm_h2o2,
    2829     &                      nuice_ref, rho_ice, rho_dust, ref_r0,
    29      &                      igcm_he
     30     &                      igcm_he, igcm_stormdust_mass,
     31     &                      igcm_stormdust_number
    3032      use comsoil_h, only: inertiedat, ! soil thermal inertia
    3133     &                     tsoil, nsoilmx ! number of subsurface layers
    32       use geometry_mod, only: longitude, latitude, cell_area
     34      use geometry_mod, only: longitude, latitude, cell_area,
     35     &                        longitude_deg
    3336      use comgeomfi_h, only: sinlon, coslon, sinlat, coslat
    3437      use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam,
    … …  
    3639     &                     frost_albedo_threshold,
    3740     &                     tsurf, co2ice, emis,
    38      &                     capcal, fluxgrd, qsurf
    39       use comsaison_h, only: dist_sol, declin, mu0, fract
     41     &                     capcal, fluxgrd, qsurf,
     42     &                     hmons
     43      use comsaison_h, only: dist_sol, declin, mu0, fract, local_time
    4044      use slope_mod, only: theta_sl, psi_sl
    4145      use conc_mod, only: rnew, cpnew, mmean
    … …  
    4852      use planete_h, only: aphelie, periheli, year_day, peri_day,
    4953     &                     obliquit
    50       USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad
     54      USE comcstfi_h, only: r, cpp, mugaz, g, rcp, pi, rad 
    5155      USE calldrag_noro_mod, ONLY: calldrag_noro
    5256      USE vdifc_mod, ONLY: vdifc
    … …  
    5458     &                      fill_data_thermos, allocate_param_thermos
    5559      use iono_h, only: allocate_param_iono
     60      use compute_dtau_mod, only: compute_dtau
    5661#ifdef MESOSCALE
    5762      use comsoil_h, only: mlayer,layer
    … …  
    6469      use eofdump_mod, only: eofdump
    6570      USE vertical_layers_mod, ONLY: ap,bp,aps,bps
     71     
    6672#endif
    6773
    … …  
    207213      REAL,INTENT(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s)
    208214
     215     
    209216
    210217c Local saved variables:
    … …  
    212219      INTEGER,SAVE :: day_ini ! Initial date of the run (sol since Ls=0)
    213220      INTEGER,SAVE :: icount     ! counter of calls to physiq during the run.
     221
    214222#ifdef DUSTSTORM
    215223      REAL pq_tmp(ngrid, nlayer, 2) ! To compute tendencies due the dust bomb
    216224#endif
     225           
    217226c     Variables used by the water ice microphysical scheme:
    218227      REAL rice(ngrid,nlayer)    ! Water ice geometric mean radius (m)
    … …  
    257266      REAL fluxtop_sw(ngrid,2)     !Outgoing SW (solar) flux to space (W.m-2)
    258267      REAL tauref(ngrid)           ! Reference column optical depth at odpref
     268c     rocket dust storm
     269      REAL totstormfract(ngrid)     ! fraction of the mesh where the dust storm is contained
     270      logical clearatm              ! clearatm used to calculate twice the radiative
     271                                    ! transfer when rdstorm is active :
     272                                    !            - in a mesh with stormdust and background dust (false)
     273                                    !            - in a mesh with background dust only (true)
     274     
    259275      real,parameter :: odpref=610. ! DOD reference pressure (Pa)
    260276      REAL tau(ngrid,naerkind)     ! Column dust optical depth at each point
    261277                                   ! AS: TBD: this one should be in a module !
    262       REAL dsodust(ngrid,nlayer)   ! density-scaled opacity (in infrared)
    263278      REAL zls                       !  solar longitude (rad)
    264279      REAL zday                      ! date (time since Ls=0, in martian days)
    … …  
    271286      REAL zdtlw(ngrid,nlayer)     ! (K/s)
    272287      REAL zdtsw(ngrid,nlayer)     ! (K/s)
    273 !      REAL cldtlw(ngrid,nlayer)     ! (K/s) LW heating rate for clear area
    274 !      REAL cldtsw(ngrid,nlayer)     ! (K/s) SW heating rate for clear area
     288      REAL pdqrds(ngrid,nlayer,nq) ! tendency for dust after rocketduststorm
     289
    275290      REAL zdtnirco2(ngrid,nlayer) ! (K/s)
    276291      REAL zdtnlte(ngrid,nlayer)   ! (K/s)
    … …  
    315330      PARAMETER (length=100)
    316331
     332c     Variables for the total dust for diagnostics
     333      REAL qdusttotal(ngrid,nlayer) !it equals to dust + stormdust
     334
     335      INTEGER iaer
     336
    317337c local variables only used for diagnostic (output in file "diagfi" or "stats")
    318338c -----------------------------------------------------------------------------
    … …  
    320340      REAL zu(ngrid,nlayer),zv(ngrid,nlayer)
    321341      REAL zq(ngrid,nlayer,nq)
     342
    322343      REAL fluxtop_sw_tot(ngrid), fluxsurf_sw_tot(ngrid)
    323344      character*2 str2
    … …  
    325346      real zdtdif(ngrid,nlayer), zdtadj(ngrid,nlayer)
    326347      real rdust(ngrid,nlayer) ! dust geometric mean radius (m)
     348      real rstormdust(ngrid,nlayer) ! stormdust geometric mean radius (m)
    327349      integer igmin, lmin
    328350      logical tdiag
    … …  
    340362      real mass(nq)             ! global mass of tracers (g)
    341363      REAL mtot(ngrid)          ! Total mass of water vapor (kg/m2)
     364      REAL mstormdtot(ngrid)    ! Total mass of stormdust tracer (kg/m2)
     365      REAL mdusttot(ngrid)    ! Total mass of dust tracer (kg/m2)
    342366      REAL icetot(ngrid)        ! Total mass of water ice (kg/m2)
    343367      REAL mtotco2(ngrid)       ! Total mass of co2 vapor (kg/m2)
    … …  
    360384      REAL nccn(ngrid,nlayer)  ! true n ccn (kg/kg)
    361385      REAL qccn(ngrid,nlayer)  ! true q ccn (kg/kg)
     386c     definition tendancies of stormdust tracers
     387      REAL rdsdqdustsurf(ngrid) ! surface q stormdust flux (kg/m2/s)
     388      REAL rdsdndustsurf(ngrid) ! surface n stormdust flux (number/m2/s)
     389      REAL rdsndust(ngrid,nlayer) ! true n stormdust (kg/kg)
     390      REAL rdsqdust(ngrid,nlayer) ! true q stormdust (kg/kg)
     391      REAL wspeed(ngrid,nlayer+1) ! vertical speed tracer stormdust
     392      REAL dsodust(ngrid,nlayer)
     393      REAL dsords(ngrid,nlayer)
     394
    362395      REAL nccnco2(ngrid,nlayer)   ! true n ccnco2 (kg/kg)
    363396      REAL qccnco2(ngrid,nlayer)  ! true q ccnco2 (kg/kg)
    … …  
    396429      REAL tstar(ngrid)  ! friction velocity and friction potential temp
    397430      REAL L_mo(ngrid),vhf(ngrid),vvv(ngrid)
    398 !      REAL zu2(ngrid)
     431      real qdustrds0(ngrid,nlayer),qdustrds1(ngrid,nlayer)
     432      real qstormrds0(ngrid,nlayer),qstormrds1(ngrid,nlayer)   
     433      real qdusttotal0(ngrid),qdusttotal1(ngrid)
    399434
    400435c sub-grid scale water ice clouds (A. Pottier 2013)
    … …  
    543578     &                 nsoilmx,ngrid,nlayer,nq,
    544579     &                 ptimestep,pday,time_phys,cell_area,
    545      &                 albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
     580     &                 albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
     581     &                 hmons)
    546582         endif
    547583#endif
    … …  
    565601      zdtsurf(:)=0
    566602      dqsurf(:,:)=0
     603     
    567604#ifdef DUSTSTORM
    568605      pq_tmp(:,:,:)=0
    … …  
    572609
    573610      zday=pday+ptime ! compute time, in sols (and fraction thereof)
    574    
     611      ! Compute local time at each grid point
     612      DO ig=1,ngrid
     613       local_time(ig)=modulo(1.+(zday-INT(zday))
     614     &                +(longitude_deg(ig)/15)/24,1.)
     615      ENDDO
    575616c     Compute Solar Longitude (Ls) :
    576617c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    702743     &                             zdtnlte(1:ngrid,1:nlayer)/86400.
    703744              endif
    704            else
     745           ELSE
    705746              zdtnlte(:,:)=0.
    706            endif
     747           ENDIF !end callnlte
    707748
    708749c          Find number of layers for LTE radiation calculations
    709750           IF(MOD(iphysiq*(icount-1),day_step).EQ.0)
    710751     &          CALL nlthermeq(ngrid,nlayer,zplev,zplay)
     752
     753c          rocketstorm : compute dust storm mesh fraction
     754           IF (rdstorm) THEN
     755              CALL calcstormfract(ngrid,nlayer,nq,pq,
     756     &                 totstormfract)
     757           ENDIF
    711758
    712759c          Note: Dustopacity.F has been transferred to callradite.F
    … …  
    724771c          Transfer through CO2 (except NIR CO2 absorption)
    725772c            and aerosols (dust and water ice)
    726 
     773           ! callradite for background dust
     774           clearatm=.true.
    727775c          Radiative transfer
    728776c          ------------------
    … …  
    730778           clearsky=.false. ! part with clouds for both cases CLFvarying true/false
    731779           CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
    732      $     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
    733      $     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
    734      $     tauref,tau,aerosol,dsodust,tauscaling,taucloudtes,rdust,rice,
    735      $     nuice,co2ice,clearsky,totcloudfrac)
     780     &     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
     781     &     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,
     782     &     fluxtop_sw,tauref,tau,aerosol,dsodust,tauscaling,
     783     &     taucloudtes,rdust,rice,nuice,co2ice,rstormdust,
     784     &     totstormfract,clearatm,dsords,
     785     &     clearsky,totcloudfrac)
     786
    736787           ! case of sub-grid water ice clouds: callradite for the clear case
    737788            IF (CLFvarying) THEN
    … …  
    746797     &              fluxsurf_swclf,fluxtop_lwclf,fluxtop_swclf,tauref,
    747798     &              tau,aerosol,dsodust,tauscaling,taucloudtesclf,rdust,
    748      &              rice,nuice,co2ice, clearsky, totcloudfrac)
     799     &              rice,nuice,co2ice,rstormdust,totstormfract,
     800     &              clearatm,dsords,clearsky,totcloudfrac)
    749801               clearsky = .false.  ! just in case.
    750802               ! Sum the fluxes and heating rates from cloudy/clear
    … …  
    774826
    775827            ENDIF ! (CLFvarying)
    776 
     828           
     829           ! Dustinjection
     830           if (dustinjection.gt.0) then
     831             CALL compute_dtau(ngrid,nlayer,
     832     &                          zday,pplev,tauref,
     833     &                          ptimestep,dustliftday,local_time)
     834           endif
     835c============================================================================
     836           
    777837#ifdef DUSTSTORM
    778838!! specific case: compute the added quantity of dust for perturbation
    779839       if (firstcall) then
    780840         pdq(1:ngrid,1:nlayer,igcm_dust_mass)=
    781      &      pdq(1:ngrid,1:nlayer,igcm_dust_mass)
     841     &      pdq(1:ngrid,1:nlayer,igcm_dust_mass)     
    782842     &      - pq_tmp(1:ngrid,1:nlayer,1)
    783843     &      + pq(1:ngrid,1:nlayer,igcm_dust_mass)
    … …  
    897957      ENDIF ! of IF (callrad)
    898958
     959c     3. Rocket dust storm
     960c    -------------------------------------------
     961      IF (rdstorm) THEN
     962         clearatm=.false.
     963         pdqrds(:,:,:)=0.
     964         qdusttotal0(:)=0.
     965         qdusttotal1(:)=0.
     966         do ig=1,ngrid
     967           do l=1,nlayer
     968            zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l) ! updated potential
     969                                             ! temperature tendency
     970            ! for diagnostics
     971            qdustrds0(ig,l)=pq(ig,l,igcm_dust_mass)+
     972     &           pdq(ig,l,igcm_dust_mass)*ptimestep
     973            qstormrds0(ig,l)=pq(ig,l,igcm_stormdust_mass)+
     974     &           pdq(ig,l,igcm_stormdust_mass)*ptimestep
     975            qdusttotal0(ig)=qdusttotal0(ig)+(qdustrds0(ig,l)+
     976     &          qstormrds0(ig,l))*(zplev(ig,l)-
     977     &           zplev(ig,l+1))/g
     978           enddo
     979         enddo
     980         call writediagfi(ngrid,'qdustrds0','qdust before rds',
     981     &           'kg/kg ',3,qdustrds0)
     982         call writediagfi(ngrid,'qstormrds0','qstorm before rds',
     983     &           'kg/kg ',3,qstormrds0)
     984
     985         CALL rocketduststorm(ngrid,nlayer,nq,ptime,ptimestep,
     986     &                      pq,pdq,pt,pdt,zplev,zplay,zzlev,
     987     &                      zzlay,zdtsw,zdtlw,
     988c               for radiative transfer
     989     &                      clearatm,icount,zday,zls,
     990     &                      tsurf,igout,totstormfract,
     991c               input sub-grid scale cloud
     992     &                      clearsky,totcloudfrac,
     993c               output
     994     &                      pdqrds,wspeed,dsodust,dsords)
     995                   
     996c      update the tendencies of both dust after vertical transport
     997         DO l=1,nlayer
     998          DO ig=1,ngrid
     999           pdq(ig,l,igcm_stormdust_mass)=
     1000     &              pdq(ig,l,igcm_stormdust_mass)+
     1001     &                      pdqrds(ig,l,igcm_stormdust_mass)
     1002           pdq(ig,l,igcm_stormdust_number)=
     1003     &              pdq(ig,l,igcm_stormdust_number)+
     1004     &                      pdqrds(ig,l,igcm_stormdust_number)
     1005
     1006           pdq(ig,l,igcm_dust_mass)=
     1007     &       pdq(ig,l,igcm_dust_mass)+ pdqrds(ig,l,igcm_dust_mass)
     1008           pdq(ig,l,igcm_dust_number)=
     1009     &           pdq(ig,l,igcm_dust_number)+
     1010     &                  pdqrds(ig,l,igcm_dust_number)
     1011
     1012          ENDDO
     1013         ENDDO
     1014         do l=1,nlayer
     1015           do ig=1,ngrid
     1016             qdustrds1(ig,l)=pq(ig,l,igcm_dust_mass)+
     1017     &           pdq(ig,l,igcm_dust_mass)*ptimestep
     1018             qstormrds1(ig,l)=pq(ig,l,igcm_stormdust_mass)+
     1019     &           pdq(ig,l,igcm_stormdust_mass)*ptimestep
     1020             qdusttotal1(ig)=qdusttotal1(ig)+(qdustrds1(ig,l)+
     1021     &          qstormrds1(ig,l))*(zplev(ig,l)-
     1022     &           zplev(ig,l+1))/g
     1023           enddo
     1024         enddo
     1025
     1026c        for diagnostics
     1027         call writediagfi(ngrid,'qdustrds1','qdust after rds',
     1028     &           'kg/kg ',3,qdustrds1)
     1029         call writediagfi(ngrid,'qstormrds1','qstorm after rds',
     1030     &           'kg/kg ',3,qstormrds1)
     1031         
     1032         call writediagfi(ngrid,'qdusttotal0','q sum before rds',
     1033     &           'kg/kg ',2,qdusttotal0)
     1034         call writediagfi(ngrid,'qdusttotal1','q sum after rds',
     1035     &           'kg/kg ',2,qdusttotal1)
     1036
     1037      ENDIF ! end of if(rdstorm)
     1038
    8991039c-----------------------------------------------------------------------
    900 c    3. Gravity wave and subgrid scale topography drag :
     1040c    4. Gravity wave and subgrid scale topography drag :
    9011041c    -------------------------------------------------
    9021042
    … …  
    9171057
    9181058c-----------------------------------------------------------------------
    919 c    4. Vertical diffusion (turbulent mixing):
     1059c    5. Vertical diffusion (turbulent mixing):
    9201060c    -----------------------------------------
    9211061
    … …  
    9681108     $        zdum1,zdum2,zdh,pdq,zflubid,
    9691109     $        zdudif,zdvdif,zdhdif,zdtsdif,q2,
    970      &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,sensibFlux)
    971 
     1110     &        zdqdif,zdqsdif,wstar,zcdv,zcdh,hfmax_th,sensibFlux,
     1111     &        dustliftday,local_time)
    9721112
    9731113          DO ig=1,ngrid
    … …  
    10181158             zdqdif(1:ngrid,2:nlayer,1:nq) = 0.
    10191159             DO iq=1, nq
    1020               IF ((iq .ne. igcm_dust_mass)
    1021      &        .and. (iq .ne. igcm_dust_number)) THEN
    1022                 zdqsdif(:,iq)=0.
    1023               ENDIF
     1160               IF ((iq .ne. igcm_dust_mass)
     1161     &          .and. (iq .ne. igcm_dust_number)) THEN
     1162                 zdqsdif(:,iq)=0.
     1163               ENDIF
    10241164             ENDDO
    10251165           ELSE
    … …  
    10411181
    10421182c-----------------------------------------------------------------------
    1043 c   5. Thermals :
     1183c   6. Thermals :
    10441184c   -----------------------------
    10451185
    … …  
    10851225        lmax_th_out(:)=0.
    10861226      end if
    1087 
    10881227c-----------------------------------------------------------------------
    1089 c   5. Dry convective adjustment:
     1228c   7. Dry convective adjustment:
    10901229c   -----------------------------
    10911230
    … …  
    11071246     $                zduadj,zdvadj,zdhadj,
    11081247     $                zdqadj)
    1109 
    11101248
    11111249         DO l=1,nlayer
    … …  
    11331271
    11341272c-----------------------------------------------------------------------
    1135 c   6. Specific parameterizations for tracers
     1273c   8. Specific parameterizations for tracers
    11361274c:   -----------------------------------------
    11371275
    11381276      if (tracer) then
    11391277
    1140 c   6a. Water and ice
     1278c   8a. Water and ice
    11411279c     ---------------
    11421280
    … …  
    12191357         END IF  ! of IF (water)
    12201358
    1221 c   6a bis. CO2 clouds (CL & JA)
     1359c   8a bis. CO2 clouds (CL & JA)
    12221360c        ---------------------------------------
    12231361c        CO2 ice cloud condensation in the atmosphere
    … …  
    13421480
    13431481
    1344 c   6b. Aerosol particles
     1482c   8b. Aerosol particles
    13451483c     -------------------
    1346 
    13471484c        ----------
    13481485c        Dust devil :
    … …  
    13711508c        -------------
    13721509c        Sedimentation :   acts also on water ice
    1373 c        ------------- 
     1510c        -------------
    13741511         IF (sedimentation) THEN
    13751512           zdqsed(1:ngrid,1:nlayer,1:nq)=0
    … …  
    13781515c Sedimentation for co2 clouds tracers are inside co2cloud microtimestep
    13791516c Zdqssed isn't
    1380            call callsedim(ngrid,nlayer, ptimestep,
    1381      &                zplev,zzlev, zzlay, pt, pdt, rdust, rice,
    1382      &                rsedcloud,rhocloud,
    1383      &                pq, pdq, zdqsed, zdqssed,nq,
     1517           call callsedim(ngrid,nlayer,ptimestep,
     1518     &                zplev,zzlev,zzlay,pt,pdt,rdust,rstormdust,
     1519     &                rice,rsedcloud,rhocloud,
     1520     &                pq,pdq,zdqsed,zdqssed,nq,
    13841521     &                tau,tauscaling)
    13851522c Flux at the surface of co2 ice computed in co2cloud microtimestep
    … …  
    13871524              zdqssed(1:ngrid,igcm_co2_ice)=zdqssed_co2(1:ngrid)
    13881525           ENDIF
     1526
     1527           IF (rdstorm) THEN
     1528c             Storm dust cannot sediment to the surface
     1529              DO ig=1,ngrid 
     1530                 zdqsed(ig,1,igcm_stormdust_mass)=
     1531     &             zdqsed(ig,1,igcm_stormdust_mass)+
     1532     &             zdqssed(ig,igcm_stormdust_mass) /
     1533     &             ((pplev(ig,1)-pplev(ig,2))/g)
     1534                 zdqsed(ig,1,igcm_stormdust_number)=
     1535     &             zdqsed(ig,1,igcm_stormdust_number)+
     1536     &             zdqssed(ig,igcm_stormdust_number) /
     1537     &               ((pplev(ig,1)-pplev(ig,2))/g)
     1538                 zdqssed(ig,igcm_stormdust_mass)=0.
     1539                 zdqssed(ig,igcm_stormdust_number)=0.
     1540              ENDDO
     1541           ENDIF !rdstorm
    13891542
    13901543           DO iq=1, nq
    … …  
    14001553             ENDDO
    14011554           ENDDO
     1555
    14021556         END IF   ! of IF (sedimentation)
    1403        
     1557
    14041558c Add lifted dust to tendancies after sedimentation in the LES (AC)
    14051559      IF (turb_resolved) THEN
    … …  
    14171571           ENDDO
    14181572      ENDIF
    1419 
    14201573c
    1421 c   6c. Chemical species
     1574c   8c. Chemical species
    14221575c     ------------------
    14231576
    … …  
    14931646#endif
    14941647
    1495 c   6d. Updates
     1648c   8d. Updates
    14961649c     ---------
    14971650
    … …  
    15111664#ifndef MESOSCALE
    15121665c-----------------------------------------------------------------------
    1513 c   7. THERMOSPHERE CALCULATION
     1666c   9. THERMOSPHERE CALCULATION
    15141667c-----------------------------------------------------------------------
    15151668
    … …  
    15351688      endif ! of if (callthermos)
    15361689#endif
    1537 
    15381690c-----------------------------------------------------------------------
    1539 c   8. Carbon dioxide condensation-sublimation:
     1691c   10. Carbon dioxide condensation-sublimation:
    15401692c     (should be the last atmospherical physical process to be computed)
    15411693c   -------------------------------------------
    … …  
    16071759#endif
    16081760     
    1609       ENDIF  ! of IF (callcond)
    1610      
    1611 
     1761      ENDIF  ! of IF (callcond)     
    16121762c-----------------------------------------------------------------------
    1613 c   9. Surface  and sub-surface soil temperature
     1763c   11. Surface  and sub-surface soil temperature
    16141764c-----------------------------------------------------------------------
    16151765c
    16161766c
    1617 c   9.1 Increment Surface temperature:
     1767c   11.1 Increment Surface temperature:
    16181768c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    16191769
    … …  
    16561806
    16571807c
    1658 c   9.2 Compute soil temperatures and subsurface heat flux:
     1808c   11.2 Compute soil temperatures and subsurface heat flux:
    16591809c   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    16601810      IF (callsoil) THEN
    … …  
    16691819        ENDIF
    16701820      ENDIF
    1671      
     1821
     1822c     To avoid negative values
     1823      IF (rdstorm) THEN
     1824           where (pq(:,:,igcm_stormdust_mass) +
     1825     &      ptimestep*pdq(:,:,igcm_stormdust_mass) < 0.)
     1826             pdq(:,:,igcm_stormdust_mass) =
     1827     &       - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30
     1828             pdq(:,:,igcm_stormdust_number) =
     1829     &       - pq(:,:,igcm_stormdust_number)/ptimestep + 1.e-30
     1830           end where
     1831           where (pq(:,:,igcm_stormdust_number) +
     1832     &      ptimestep*pdq(:,:,igcm_stormdust_number) < 0.)
     1833             pdq(:,:,igcm_stormdust_mass) =
     1834     &       - pq(:,:,igcm_stormdust_mass)/ptimestep + 1.e-30
     1835             pdq(:,:,igcm_stormdust_number) =
     1836     &       - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
     1837           end where
     1838
     1839            where (pq(:,:,igcm_dust_mass) +
     1840     &      ptimestep*pdq(:,:,igcm_dust_mass) < 0.)
     1841             pdq(:,:,igcm_dust_mass) =
     1842     &       - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
     1843             pdq(:,:,igcm_dust_number) =
     1844     &       - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
     1845           end where
     1846           where (pq(:,:,igcm_dust_number) +
     1847     &      ptimestep*pdq(:,:,igcm_dust_number) < 0.)
     1848             pdq(:,:,igcm_dust_mass) =
     1849     &       - pq(:,:,igcm_dust_mass)/ptimestep + 1.e-30
     1850             pdq(:,:,igcm_dust_number) =
     1851     &       - pq(:,:,igcm_dust_number)/ptimestep + 1.e-30
     1852           end where
     1853      ENDIF !(rdstorm)   
    16721854
    16731855c-----------------------------------------------------------------------
    1674 c  10. Write output files
     1856c  12. Write output files
    16751857c  ----------------------
    16761858
    … …  
    18852067              endif
    18862068           endif ! of (doubleq)
     2069
     2070           if (rdstorm) then   ! diagnostics of stormdust tendancies for 1D and 3D
     2071              mstormdtot(:)=0
     2072              mdusttot(:)=0
     2073              qdusttotal(:,:)=0
     2074              do ig=1,ngrid
     2075                rdsdqdustsurf(ig) =
     2076     &                zdqssed(ig,igcm_stormdust_mass)*tauscaling(ig)
     2077                rdsdndustsurf(ig) =
     2078     &                zdqssed(ig,igcm_stormdust_number)*tauscaling(ig)
     2079                rdsndust(ig,:) =
     2080     &                pq(ig,:,igcm_stormdust_number)*tauscaling(ig)
     2081                rdsqdust(ig,:) =
     2082     &                pq(ig,:,igcm_stormdust_mass)*tauscaling(ig)
     2083                do l=1,nlayer
     2084                    mstormdtot(ig) = mstormdtot(ig) +
     2085     &                      zq(ig,l,igcm_stormdust_mass) *
     2086     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
     2087                    mdusttot(ig) = mdusttot(ig) +
     2088     &                      zq(ig,l,igcm_dust_mass) *
     2089     &                      (zplev(ig,l) - zplev(ig,l+1)) / g
     2090                    qdusttotal(ig,l) = qdust(ig,l)+rdsqdust(ig,l) !calculate total dust
     2091                enddo
     2092              enddo
     2093           endif !(rdstorm)
     2094
    18872095           if (co2clouds) then
    18882096              do ig=1,ngrid
    … …  
    19062114              enddo
    19072115           endif ! of if (co2clouds)
    1908            
    1909                    
     2116                             
    19102117           if (water) then
    19112118             mtot(:)=0
    … …  
    21382345               call wstats(ngrid,'dustN','Dust number',
    21392346     &                        'part/kg',3,ndust)
     2347               if (rdstorm) then
     2348                 call wstats(ngrid,'reffstormdust','reffdust',
     2349     &                          'm',3,rstormdust*ref_r0)
     2350                 call wstats(ngrid,'rdsdustq','Dust mass mr',
     2351     &                          'kg/kg',3,rdsqdust)
     2352                 call wstats(ngrid,'rdsdustN','Dust number',
     2353     &                        'part/kg',3,rdsndust)
     2354               end if
    21402355             else
    21412356               do iq=1,dustbin
    … …  
    23452560c        call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
    23462561c    &                   'w.m-2',3,zdtlw)
    2347   
     2562 
    23482563            if (.not.activice) then
    23492564               CALL WRITEDIAGFI(ngrid,'tauTESap',
    … …  
    24342649     &                       qsurf(1:ngrid,igcm_h2o_ice))
    24352650#endif
    2436 
    24372651            CALL WRITEDIAGFI(ngrid,'mtot',
    24382652     &                       'total mass of water vapor',
    … …  
    24432657            vmr=zq(1:ngrid,1:nlayer,igcm_h2o_ice)
    24442658     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_ice)
    2445             call WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr',
     2659            CALL WRITEDIAGFI(ngrid,'vmr_h2oice','h2o ice vmr',
    24462660     &                       'mol/mol',3,vmr)
    24472661            vmr=zq(1:ngrid,1:nlayer,igcm_h2o_vap)
    24482662     &      *mmean(1:ngrid,1:nlayer)/mmol(igcm_h2o_vap)
    2449             call WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr',
     2663            CALL WRITEDIAGFI(ngrid,'vmr_h2ovap','h2o vap vmr',
    24502664     &                       'mol/mol',3,vmr)
    24512665            CALL WRITEDIAGFI(ngrid,'reffice',
    … …  
    25552769     &                        'part/kg',3,ndust)
    25562770             call WRITEDIAGFI(ngrid,'dsodust',
    2557      &                        'dust density scaled opacity',
    2558      &                        'm2.kg-1',3,dsodust)
    2559 c             call WRITEDIAGFI(ngrid,"tauscaling",
    2560 c     &                    "dust conversion factor"," ",2,tauscaling)
     2771     &                       'density scaled optcial depth',
     2772     &                       'm2.kg-1',3,dsodust)
     2773             call WRITEDIAGFI(ngrid,'dso',
     2774     &                       'density scaled optcial depth',
     2775     &                       'm2.kg-1',3,dsodust+dsords)
    25612776           else
    25622777             do iq=1,dustbin
    … …  
    25692784           endif ! (doubleq)
    25702785
     2786           if (rdstorm) then  ! writediagfi tendencies stormdust tracers
     2787             call WRITEDIAGFI(ngrid,'reffdust','reffdust',
     2788     &                        'm',3,rdust*ref_r0)
     2789             call WRITEDIAGFI(ngrid,'reffstormdust','reffstormdust',
     2790     &                        'm',3,rstormdust*ref_r0)
     2791             call WRITEDIAGFI(ngrid,'mstormdtot',
     2792     &                        'total mass of stormdust only',
     2793     &                        'kg.m-2',2,mstormdtot)
     2794             call WRITEDIAGFI(ngrid,'mdusttot',
     2795     &                        'total mass of dust only',
     2796     &                        'kg.m-2',2,mdusttot)
     2797             call WRITEDIAGFI(ngrid,'rdsdqsdust',
     2798     &                        'deposited surface stormdust mass',
     2799     &                        'kg.m-2.s-1',2,rdsdqdustsurf)
     2800             call WRITEDIAGFI(ngrid,'rdsdustq','storm Dust mass mr',
     2801     &                        'kg/kg',3,rdsqdust)
     2802             call WRITEDIAGFI(ngrid,'rdsdustqmodel','storm Dust massmr',
     2803     &                        'kg/kg',3,pq(:,:,igcm_stormdust_mass))
     2804             call WRITEDIAGFI(ngrid,'rdsdustN','storm Dust number',
     2805     &                        'part/kg',3,rdsndust)
     2806             call WRITEDIAGFI(ngrid,"stormfract",
     2807     &                "fraction of the mesh, with stormdust","none",
     2808     &                                               2,totstormfract)
     2809             call WRITEDIAGFI(ngrid,'qsurf',
     2810     &                 'stormdust injection',
     2811     &                 'kg.m-2',2,qsurf(:,igcm_stormdust_mass))
     2812             call WRITEDIAGFI(ngrid,'pdqsurf',
     2813     &                  'tendancy stormdust mass at surface',
     2814     &                  'kg.m-2',2,dqsurf(:,igcm_stormdust_mass))
     2815             call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust',
     2816     &                        'm/s',3,wspeed(:,1:nlayer))
     2817             call WRITEDIAGFI(ngrid,'zdqsed_dustq'
     2818     &          ,'sedimentation q','kg.m-2.s-1',3,
     2819     &          zdqsed(:,:,igcm_dust_mass))
     2820             call WRITEDIAGFI(ngrid,'zdqssed_dustq'
     2821     &          ,'sedimentation q','kg.m-2.s-1',2,
     2822     &          zdqssed(:,igcm_dust_mass))
     2823             call WRITEDIAGFI(ngrid,'zdqsed_rdsq'
     2824     &          ,'sedimentation q','kg.m-2.s-1',3,
     2825     &          zdqsed(:,:,igcm_stormdust_mass))
     2826             call WRITEDIAGFI(ngrid,'rdust','rdust',
     2827     &                        'm',3,rdust)
     2828             call WRITEDIAGFI(ngrid,'rstormdust','rstormdust',
     2829     &                        'm',3,rstormdust)
     2830             call WRITEDIAGFI(ngrid,'totaldustq','total dust mass',
     2831     &                        'kg/kg',3,qdusttotal)
     2832             call WRITEDIAGFI(ngrid,'dsords',
     2833     &                       'density scaled opacity of stormdust',
     2834     &                       'm2.kg-1',3,dsords)
     2835           endif ! (rdstorm)
     2836
    25712837           if (scavenging) then
    25722838             call WRITEDIAGFI(ngrid,'ccnq','CCN mass mr',
    … …  
    25742840             call WRITEDIAGFI(ngrid,'ccnN','CCN number',
    25752841     &                        'part/kg',3,nccn)
     2842             call WRITEDIAGFI(ngrid,'surfccnq','Surf nuclei mass mr',
     2843     &                        'kg.m-2',2,qsurf(1,igcm_ccn_mass))
     2844             call WRITEDIAGFI(ngrid,'surfccnN','Surf nuclei number',
     2845     &                        'kg.m-2',2,qsurf(1,igcm_ccn_number))
    25762846           endif ! (scavenging)
    25772847
    … …  
    28043074             call WRITEDIAGFI(ngrid,'rdust','rdust',
    28053075     &                        'm',1,rdust)
    2806            endif
     3076           endif ! doubleq 1D
     3077           if (rdstorm) then
     3078             call writediagfi(1,'aerosol_dust','opacity of env. dust',''
     3079     &                        ,1,aerosol(:,:,igcm_dust_mass))
     3080             call writediagfi(1,'aerosol_stormdust',
     3081     &                         'opacity of storm dust',''
     3082     &                        ,1,aerosol(:,:,igcm_stormdust_mass))
     3083             call WRITEDIAGFI(ngrid,'dqsdifdustq','diffusion',
     3084     &                       'kg.m-2.s-1',0,zdqsdif(1,igcm_dust_mass))
     3085             call WRITEDIAGFI(ngrid,'dqsdifrdsq','diffusion',
     3086     &                  'kg.m-2.s-1',0,zdqsdif(1,igcm_stormdust_mass))
     3087             call WRITEDIAGFI(ngrid,'mstormdtot',
     3088     &                        'total mass of stormdust only',
     3089     &                        'kg.m-2',0,mstormdtot)
     3090             call WRITEDIAGFI(ngrid,'mdusttot',
     3091     &                        'total mass of dust only',
     3092     &                        'kg.m-2',0,mdusttot)
     3093             call WRITEDIAGFI(ngrid,'tauref',
     3094     &                    'Dust ref opt depth','NU',0,tauref)
     3095             call WRITEDIAGFI(ngrid,'rdsdqsdust',
     3096     &                        'deposited surface stormdust mass',
     3097     &                        'kg.m-2.s-1',0,rdsdqdustsurf)
     3098             call WRITEDIAGFI(ngrid,'rdsdustq','storm Dust mass mr',
     3099     &                        'kg/kg',1,rdsqdust)
     3100             call WRITEDIAGFI(ngrid,"stormfract",
     3101     &                         "fraction of the mesh,with stormdust",
     3102     &                         "none",0,totstormfract)
     3103             call WRITEDIAGFI(ngrid,'rdsqsurf',
     3104     &                 'stormdust at surface',
     3105     &                 'kg.m-2',0,qsurf(:,igcm_stormdust_mass))
     3106             call WRITEDIAGFI(ngrid,'qsurf',
     3107     &                  'dust mass at surface',
     3108     &                  'kg.m-2',0,qsurf(:,igcm_dust_mass))
     3109             call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust',
     3110     &                        'm/s',1,wspeed)
     3111             call WRITEDIAGFI(ngrid,'totaldustq','total dust mass',
     3112     &                        'kg/kg',1,qdusttotal)
     3113             call WRITEDIAGFI(ngrid,'dsords',
     3114     &                       'density scaled opacity of stormdust',
     3115     &                       'm2.kg-1',1,dsords)
     3116             call WRITEDIAGFI(ngrid,'zdqsed_dustq'
     3117     &          ,'sedimentation q','kg.m-2.s-1',1,
     3118     &          zdqsed(1,:,igcm_dust_mass))
     3119             call WRITEDIAGFI(ngrid,'zdqsed_rdsq'
     3120     &          ,'sedimentation q','kg.m-2.s-1',1,
     3121     &          zdqsed(1,:,igcm_stormdust_mass))
     3122           endif !(rdstorm 1D)
     3123
    28073124           if (water.AND.tifeedback) then
    28083125             call WRITEDIAGFI(ngrid,"soiltemp",
    … …  
    29753292
    29763293      icount=icount+1
     3294       
    29773295
    29783296      END SUBROUTINE physiq

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

    r1918 r1974  
    66                    iaer_dust_conrath,iaer_dust_doubleq,&
    77                    iaer_dust_submicron,iaer_h2o_ice,&
     8                    iaer_stormdust_doubleq,&
    89                    file_id,radiustab, gvis, omegavis, &
    910                    QVISsQREF, gIR, omegaIR, &
    … …  
    153154!       Reference wavelength in the infrared:
    154155        longrefir(iaer)=12.1E-6  ! 825cm-1 TES/MGS
     156!==================================================================
     157        CASE("stormdust_doubleq") aerkind   ! Two-moment scheme for stormdust - radiative properties
     158!==================================================================
     159!       Visible domain:
     160        file_id(iaer,1) = 'optprop_dustvis_TM_n50.dat' !T-Matrix
     161!       Infrared domain:
     162        file_id(iaer,2) = 'optprop_dustir_n50.dat'     !Mie
     163!       Reference wavelength in the visible:
     164        longrefvis(iaer)=0.67E-6
     165!       If not equal to 0.67e-6 -> change readtesassim accordingly;
     166!       Reference wavelength in the infrared:
     167        longrefir(iaer)=dustrefir
    155168!==================================================================
    156169        END SELECT aerkind

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

    r1770 r1974  
    2020  real,save :: iceradius(2) , dtemisice(2)
    2121  real,save,allocatable :: zmea(:),zstd(:),zsig(:),zgam(:),zthe(:)
     22  real,save,allocatable :: hmons(:)
    2223  real,save,allocatable :: z0(:) ! surface roughness length (m)
    2324  real,save :: z0_default ! default (constant over planet) surface roughness (m)
    … …  
    5556    allocate(capcal(ngrid))
    5657    allocate(fluxgrd(ngrid))
    57  
     58    allocate(hmons(ngrid))
     59   
    5860  end subroutine ini_surfdat_h
    5961
    … …  
    7981    if (allocated(capcal))      deallocate(capcal)
    8082    if (allocated(fluxgrd))     deallocate(fluxgrd)
     83    if (allocated(hmons))       deallocate(hmons)
    8184
    8285  end subroutine end_surfdat_h

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

    r1941 r1974  
    66      integer,save :: nqmx ! initialized in conf_phys
    77   
    8       character*20,allocatable,save ::  noms(:)  ! name of the tracer
     8      character*30,allocatable,save ::  noms(:)  ! name of the tracer
    99      real,allocatable,save :: mmol(:)           ! mole mass of tracer (g/mol-1)
    1010      real,allocatable,save :: radius(:)   ! dust and ice particle radius (m)
    … …  
    2323      real,save :: nuiceco2_sed   ! Sedimentation effective variance of the co2 ice dist.
    2424      real,save :: nuiceco2_ref   ! Effective variance of the co2 ice dist.
    25 
     25     
    2626      real,save :: ccn_factor  ! ratio of nuclei for water ice particles
    2727
    … …  
    4242      integer,save :: igcm_ccn_number ! CCN number mixing ratio
    4343      integer,save :: igcm_dust_submicron ! submicron dust mixing ratio
    44    
     44      integer,save :: igcm_stormdust_mass   !  storm dust mass mixing ratio
     45      integer,save :: igcm_stormdust_number !  storm dust number mixing ratio
     46   
    4547      integer,save :: igcm_ccnco2_mass   ! CCN (dust and/or water ice) for CO2 mass mixing ratio
    4648      integer,save :: igcm_ccnco2_number ! CCN (dust and/or water ice) for CO2 number mixing ratio

  • trunk/LMDZ.MARS/libf/phymars/updatereffrad_mod.F

    r1969 r1974  
    66     
    77      SUBROUTINE updatereffrad(ngrid,nlayer,
    8      &                rdust,rice,nuice,
     8     &                rdust,rstormdust,rice,nuice,
    99     &                reffrad,nueffrad,
    1010     &                pq,tauscaling,tau,pplay)
    … …  
    1414     &                       igcm_h2o_ice, igcm_ccn_mass, radius,
    1515     &                       igcm_ccn_number, nuice_ref, varian,
    16      &                       ref_r0, igcm_dust_submicron
     16     &                       ref_r0, igcm_dust_submicron,
     17     &                       igcm_stormdust_mass,igcm_stormdust_number
    1718       USE dimradmars_mod, only: nueffdust,naerkind,
    1819     &            name_iaer,
    1920     &            iaer_dust_conrath,iaer_dust_doubleq,
    20      &            iaer_dust_submicron,iaer_h2o_ice
     21     &            iaer_dust_submicron,iaer_h2o_ice,
     22     &            iaer_stormdust_doubleq
    2123
    2224       IMPLICIT NONE
    … …  
    4547
    4648c-----------------------------------------------------------------------
    47 c     Inputs:
     49c     Inputs/outputs:
    4850c     ------
    4951
    50       INTEGER ngrid,nlayer
     52      INTEGER, INTENT(in) :: ngrid,nlayer
    5153c     Ice geometric mean radius (m)
    52       REAL :: rice(ngrid,nlayer)
     54      REAL, INTENT(out) :: rice(ngrid,nlayer)
    5355c     Estimated effective variance of the size distribution (n.u.)
    54       REAL :: nuice(ngrid,nlayer)
     56      REAL, INTENT(out) :: nuice(ngrid,nlayer)
    5557c     Tracer mass mixing ratio (kg/kg)
    56       REAL pq(ngrid,nlayer,nqmx)
    57       REAL rdust(ngrid,nlayer) ! Dust geometric mean radius (m)
    58      
    59       REAL pplay(ngrid,nlayer) ! altitude at the middle of the layers
    60       REAL tau(ngrid,naerkind)
    61 
    62 
    63 c     Outputs:
    64 c     -------
    65 
     58      REAL, INTENT(in) :: pq(ngrid,nlayer,nqmx)
     59      REAL, INTENT(out) :: rdust(ngrid,nlayer) ! Dust geometric mean radius (m)
     60      REAL, INTENT(out) :: rstormdust(ngrid,nlayer) ! Dust geometric mean radius (m)
     61      REAL, INTENT(in) :: pplay(ngrid,nlayer) ! altitude at the middle of the layers
     62      REAL, INTENT(in) :: tau(ngrid,naerkind)
    6663c     Aerosol effective radius used for radiative transfer (meter)
    67       REAL :: reffrad(ngrid,nlayer,naerkind)
     64      REAL, INTENT(out) :: reffrad(ngrid,nlayer,naerkind)
    6865c     Aerosol effective variance used for radiative transfer (n.u.)
    69       REAL :: nueffrad(ngrid,nlayer,naerkind)
    70 
     66      REAL, INTENT(out) :: nueffrad(ngrid,nlayer,naerkind)
     67      REAL, INTENT(in) :: tauscaling(ngrid)         ! Convertion factor for qccn and Nccn
     68     
    7169c     Local variables:
    7270c     ---------------
    … …  
    8583      REAL Mo,No                       ! Mass and number of ccn
    8684      REAL rhocloud(ngrid,nlayer)  ! Cloud density (kg.m-3)
    87       REAL tauscaling(ngrid)         ! Convertion factor for qccn and Nccn
    8885
    8986      LOGICAL,SAVE :: firstcall=.true.
    … …  
    114111          ENDDO
    115112        ENDIF
     113
     114        ! updating radius of stormdust particles
     115        IF (rdstorm.AND.active) THEN
     116          DO l=1,nlayer
     117            DO ig=1, ngrid
     118              call updaterdust(pq(ig,l,igcm_stormdust_mass),
     119     &                 pq(ig,l,igcm_stormdust_number),rstormdust(ig,l))
     120              nueffdust(ig,l) = exp(varian**2.)-1.
     121             ENDDO
     122           ENDDO
     123        ENDIF
    116124       
    117125c       1.2 Water-ice particles
    … …  
    126134
    127135          IF (firstcall) THEN
    128             IF (minval(tauscaling).lt.0) tauscaling(:) = 1.e-3 ! default value when non-read in startfi is -1
    129             IF (freedust)                tauscaling(:) = 1.    ! if freedust, enforce no rescaling at all
     136            !IF (minval(tauscaling).lt.0) tauscaling(:) = 1.e-3 ! default value when non-read in startfi is -1
     137            !IF (freedust)                tauscaling(:) = 1.    ! if freedust, enforce no rescaling at all
    130138            firstcall = .false.
    131139          ENDIF
    … …  
    206214          ENDDO
    207215c==================================================================
     216        CASE("stormdust_doubleq") aerkind! Two-moment scheme for
     217c       stormdust; same distribution than normal dust
     218c==================================================================
     219          DO l=1,nlayer
     220            DO ig=1,ngrid
     221              reffrad(ig,l,iaer) = rstormdust(ig,l) * ref_r0
     222              nueffrad(ig,l,iaer) = nueffdust(ig,l)
     223            ENDDO
     224          ENDDO
     225c==================================================================
    208226        END SELECT aerkind
    209227      ENDDO ! iaer (loop on aerosol kind)

  • trunk/LMDZ.MARS/libf/phymars/vdifc_mod.F

    r1969 r1974  
    1212     $                pdudif,pdvdif,pdhdif,pdtsrf,pq2,
    1313     $                pdqdif,pdqsdif,wstar,zcdv_true,zcdh_true,
    14      $                hfmax,sensibFlux)
     14     $                hfmax,sensibFlux,dustliftday,local_time)
    1515
    1616      use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
    1717     &                      igcm_dust_submicron, igcm_h2o_vap,
    18      &                      igcm_h2o_ice, alpha_lift
     18     &                      igcm_h2o_ice, alpha_lift,
     19     &                      igcm_stormdust_mass, igcm_stormdust_number
    1920      use surfdat_h, only: watercaptag, frost_albedo_threshold, dryness
    2021      USE comcstfi_h, ONLY: cpp, r, rcp, g
    2122      use turb_mod, only: turb_resolved, ustar, tstar
     23      use compute_dtau_mod, only: ti_injection,tf_injection
    2224      IMPLICIT NONE
    2325
    … …  
    7981      real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq)
    8082      real,intent(out) :: pdqdif(ngrid,nlay,nq)
    81       real,intent(out) :: pdqsdif(ngrid,nq)
     83      real,intent(out) :: pdqsdif(ngrid,nq)
     84      REAL,INTENT(in) :: dustliftday(ngrid)
     85      REAL,INTENT(in) :: local_time(ngrid)
    8286     
    8387c   local:
    … …  
    691695             end do
    692696           else if (doubleq) then
    693              do ig=1,ngrid
     697             if (dustinjection.eq.0) then !injection scheme 0 (old)
     698                                          !or 2 (injection in CL)
     699              do ig=1,ngrid
    694700               if(co2ice(ig).lt.1) then ! pas de soulevement si glace CO2
    695701                 pdqsdif(ig,igcm_dust_mass) =
    … …  
    698704     &             -alpha_lift(igcm_dust_number)
    699705               end if
    700              end do
     706              end do
     707             elseif(dustinjection.eq.1)then ! dust injection scheme = 1 injection from surface
     708              do ig=1,ngrid
     709               if(co2ice(ig).lt.1) then ! pas de soulevement si glace CO2
     710                IF((ti_injection.LE.local_time(ig)).and.
     711     &                  (local_time(ig).LE.tf_injection)) THEN
     712                  if (rdstorm) then !Rocket dust storm scheme
     713                        pdqsdif(ig,igcm_stormdust_mass) =
     714     &                          -alpha_lift(igcm_stormdust_mass)
     715     &                          *dustliftday(ig)
     716                        pdqsdif(ig,igcm_stormdust_number) =
     717     &                          -alpha_lift(igcm_stormdust_number)
     718     &                          *dustliftday(ig)
     719                        pdqsdif(ig,igcm_dust_mass)= 0.
     720                        pdqsdif(ig,igcm_dust_number)= 0.
     721                  else
     722                        pdqsdif(ig,igcm_dust_mass)=
     723     &                          -dustliftday(ig)*
     724     &                          alpha_lift(igcm_dust_mass)               
     725                        pdqsdif(ig,igcm_dust_number)=
     726     &                          -dustliftday(ig)*
     727     &                          alpha_lift(igcm_dust_number)
     728                  endif
     729                  if (submicron) then
     730                        pdqsdif(ig,igcm_dust_submicron) = 0.
     731                  endif ! if (submicron)
     732                ELSE ! outside dust injection time frame
     733                        pdqsdif(ig,igcm_dust_mass)= 0.
     734                        pdqsdif(ig,igcm_dust_number)= 0.
     735                        pdqsdif(ig,igcm_stormdust_mass)= 0.
     736                        pdqsdif(ig,igcm_stormdust_number)= 0.
     737                ENDIF
     738               
     739                end if ! of if(co2ice(ig).lt.1)
     740              end do
     741             endif ! end if dustinjection
    701742           else if (submicron) then
    702743             do ig=1,ngrid
    … …  
    723764c      Inversion pour l'implicite sur q
    724765c       --------------------------------
    725         do iq=1,nq
     766        do iq=1,nq  !for all tracers including stormdust
    726767          CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
    727768
Note: See TracChangeset for help on using the changeset viewer.