Ignore:
Timestamp:
Jul 18, 2013, 8:33:04 AM (12 years ago)
Author:
emillour
Message:

Mars GCM:

  • Bug fix: when running with photochemistry, ccns did not sediment! Fixed initracer.F. Also added that callsedim/newsedim use updated temperatures.
  • Converted run0 and run_mcd scripts to bash.

EM

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r740 r1005  
    11      SUBROUTINE callsedim(ngrid,nlay, ptimestep,
    2      &                pplev,zlev, zlay, pt, rdust, rice,
     2     &                pplev,zlev, zlay, pt, pdt, rdust, rice,
    33     &                rsedcloud,rhocloud,
    44     &                pq, pdqfi, pdqsed,pdqs_sed,nq,
    … …  
    66! to use  'getin'
    77      USE ioipsl_getincom
    8       USE updaterad
     8      USE updaterad,only: updaterdust,updaterice_micro,updaterice_typ
    99      IMPLICIT NONE
    1010
    … …  
    3535c   ----------
    3636
    37       INTEGER ngrid              ! number of horizontal grid points
    38       INTEGER nlay               ! number of atmospheric layers
    39       REAL ptimestep             ! physics time step (s)
    40       REAL pplev(ngrid,nlay+1)   ! pressure at inter-layers (Pa)
    41       REAL pt(ngrid,nlay)        ! temperature at mid-layer (K)
    42       REAL zlev(ngrid,nlay+1)    ! altitude at layer boundaries
     37      integer,intent(in) :: ngrid  ! number of horizontal grid points
     38      integer,intent(in) :: nlay   ! number of atmospheric layers
     39      real,intent(in) :: ptimestep ! physics time step (s)
     40      real,intent(in) :: pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa)
     41      real,intent(in) :: zlev(ngrid,nlay+1) ! altitude at layer boundaries
     42      real,intent(in) :: zlay(ngrid,nlay)   ! altitude at the middle of the layers
     43      real,intent(in) :: pt(ngrid,nlay) ! temperature at mid-layer (K)
     44      real,intent(in) :: pdt(ngrid,nlay) ! tendency on temperature, from
     45                                         ! previous processes (K/s)
    4346c    Aerosol radius provided by the water ice microphysical scheme:
    44       REAL rdust(ngrid,nlay)     ! Dust geometric mean radius (m)
    45       REAL rice(ngrid,nlay)      ! Ice geometric mean radius (m)
    46 
     47      real,intent(out) :: rdust(ngrid,nlay) ! Dust geometric mean radius (m)
     48      real,intent(out) :: rice(ngrid,nlay)  ! H2O Ice geometric mean radius (m)
     49c     Sedimentation radius of water ice
     50      real,intent(in) :: rsedcloud(ngridmx,nlayermx)
     51c     Cloud density (kg.m-3)
     52      real,intent(inout) :: rhocloud(ngridmx,nlayermx)
    4753c    Traceurs :
    48       integer nq             ! number of tracers
    49       real pq(ngrid,nlay,nq)  ! tracers (kg/kg)
    50       real pdqfi(ngrid,nlay,nq)  ! tendency before sedimentation (kg/kg.s-1)
    51       real pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
    52       real pdqs_sed(ngrid,nq)    ! flux at surface (kg.m-2.s-1)
     54      real,intent(in) :: pq(ngrid,nlay,nq)  ! tracers (kg/kg)
     55      real,intent(in) :: pdqfi(ngrid,nlay,nq)  ! tendency before sedimentation (kg/kg.s-1)
     56      real,intent(out) :: pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
     57      real,intent(out) :: pdqs_sed(ngrid,nq)    ! flux at surface (kg.m-2.s-1)
     58      integer,intent(in) :: nq  ! number of tracers
     59      real,intent(in) :: tau(ngrid,nlay) ! dust opacity
     60      real,intent(in) :: tauscaling(ngrid)
    5361     
    5462c   local:
    … …  
    5765      INTEGER l,ig, iq
    5866
    59       real zqi(ngridmx,nlayermx,nqmx) ! to locally store tracers
    60       real masse (ngridmx,nlayermx) ! Layer mass (kg.m-2)
    61       real epaisseur (ngridmx,nlayermx) ! Layer thickness (m)
    62       real wq(ngridmx,nlayermx+1) ! displaced tracer mass (kg.m-2)
    63       real r0(ngridmx,nlayermx) ! geometric mean radius used for
     67      real zqi(ngrid,nlay,nq) ! to locally store tracers
     68      real zt(ngrid,nlay) ! to locally store temperature
     69      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
     70      real epaisseur (ngrid,nlay) ! Layer thickness (m)
     71      real wq(ngrid,nlay+1) ! displaced tracer mass (kg.m-2)
     72      real r0(ngrid,nlay) ! geometric mean radius used for
    6473                                !   sedimentation (m)
    65       real r0dust(ngridmx,nlayermx) ! geometric mean radius used for
     74      real r0dust(ngrid,nlay) ! geometric mean radius used for
    6675                                    !   dust (m)
    67 !      real r0ccn(ngridmx,nlayermx)  ! geometric mean radius used for
     76!      real r0ccn(ngrid,nlay)  ! geometric mean radius used for
    6877!                                    !   CCNs (m)
    69 c     Sedimentation radius of water ice
    70       real rsedcloud(ngridmx,nlayermx)
    71       real beta      ! correction for the shape of the ice particles (cf. newsedim)
    72       save beta
    73 c     Cloud density (kg.m-3)
    74       real rhocloud(ngridmx,nlayermx)
    75      
     78      real,save :: beta ! correction for the shape of the ice particles (cf. newsedim)
     79
    7680c     for ice radius computation
    7781      REAL Mo,No
    78       REAL tau(ngrid,nlay), tauscaling(ngrid)
    79       REAL zlay(ngrid,nlay)   ! altitude at the middle of the layers
    8082      REAl ccntyp
    8183
    … …  
    8688c       1) Parameters used to represent the changes in fall
    8789c          velocity as a function of particle size;
    88       integer nr,ir
    89       parameter (nr=12) !(nr=7) ! number of bins
    90       real rd(nr),qr(ngridmx,nlayermx,nr)
    91       real rdi(nr+1)    ! extreme and intermediate radii
    92       real Sq(ngridmx,nlayermx)
    93       real rdmin,rdmax,rdimin,rdimax
    94       data rdmin/1.e-8/ !/1.e-7/
    95       data rdmax/30.e-6/
    96       data rdimin/1.e-10/
    97       data rdimax/1e-4/
    98       save rd, rdi
     90      integer ir
     91      integer,parameter :: nr=12 !(nr=7) ! number of bins
     92      real,save :: rd(nr)
     93      real qr(ngrid,nlay,nr)
     94      real,save :: rdi(nr+1)    ! extreme and intermediate radii
     95      real Sq(ngrid,nlay)
     96      real,parameter :: rdmin=1.e-8
     97      real,parameter :: rdmax=30.e-6
     98      real,parameter :: rdimin=1.e-8 ! 1.e-7
     99      real,parameter :: rdimax=1.e-4
    99100
    100101c       2) Second size distribution for the log-normal integration
    101102c          (the mass mixing ratio is computed for each radius)
    102103
    103       integer ninter, iint
    104       parameter (ninter=4)  ! nombre de point entre chaque rayon rdi
    105       real rr(ninter,nr)
    106       save rr
     104      integer iint
     105      integer,parameter :: ninter=4 ! number of points between each rdi radii
     106      real,save :: rr(ninter,nr)
    107107      integer radpower
    108108      real sigma0
    … …  
    110110c       3) Other local variables used in doubleq
    111111
    112       INTEGER idust_mass   ! index of tracer containing dust mass
    113                            !   mix. ratio
    114       INTEGER idust_number ! index of tracer containing dust number
    115                            !   mix. ratio
    116       INTEGER iccn_mass    ! index of tracer containing CCN mass
    117                            !   mix. ratio
    118       INTEGER iccn_number  ! index of tracer containing CCN number
    119                            !   mix. ratio
    120       SAVE idust_mass,idust_number
    121       SAVE iccn_mass,iccn_number
    122      
    123 
    124 c     Firstcall:
    125 
    126       LOGICAL firstcall
    127       SAVE firstcall
    128       DATA firstcall/.true./
     112      INTEGER,SAVE :: idust_mass  ! index of tracer containing dust mass
     113                                  !   mix. ratio
     114      INTEGER,SAVE :: idust_number ! index of tracer containing dust number
     115                                   !   mix. ratio
     116      INTEGER,SAVE :: iccn_mass  ! index of tracer containing CCN mass
     117                                 !   mix. ratio
     118      INTEGER,SAVE :: iccn_number ! index of tracer containing CCN number
     119                                  !   mix. ratio
     120
     121      LOGICAL,SAVE :: firstcall=.true.
    129122
    130123c    ** un petit test de coherence
    … …  
    169162          if (noms(iq).eq."dust_mass") then
    170163            idust_mass=iq
     164            write(*,*)"callsedim: idust_mass=",idust_mass
    171165          endif
    172166          if (noms(iq).eq."dust_number") then
    173167            idust_number=iq
     168            write(*,*)"callsedim: idust_number=",idust_number
    174169          endif
    175170        enddo
    … …  
    190185            if (noms(iq).eq."ccn_mass") then
    191186              iccn_mass=iq
     187              write(*,*)"callsedim: iccn_mass=",iccn_mass
    192188            endif
    193189            if (noms(iq).eq."ccn_number") then
    194190              iccn_number=iq
     191              write(*,*)"callsedim: iccn_number=",iccn_number
    195192            endif
    196193          enddo
    … …  
    223220c    -----------------
    224221
    225       zqi(1:ngrid,1:nlay,1:nqmx) = 0.
    226 c     Updating the mass mixing ratio with the tendencies coming
     222!      zqi(1:ngrid,1:nlay,1:nqmx) = 0.
     223c     Update the mass mixing ratio and temperature with the tendencies coming
    227224c       from other parameterizations:
    228225c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    229 
    230       do iq=1,nq
    231         do l=1,nlay
    232           do ig=1,ngrid
    233             zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
    234           enddo
    235         enddo
    236       enddo
     226      zqi(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq)
     227     &                         +pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep
     228      zt(1:ngrid,1:nlay)=pt(1:ngrid,1:nlay)
     229     &                         +pdt(1:ngrid,1:nlay)*ptimestep
     230
    237231
    238232c    Computing the different layer properties
    … …  
    342336         
    343337               call newsedim(ngrid,nlay,1,1,ptimestep,
    344      &         pplev,masse,epaisseur,pt,rd(ir),rho_dust,qr(1,1,ir),
     338     &         pplev,masse,epaisseur,zt,rd(ir),rho_dust,qr(1,1,ir),
    345339     &         wq,0.5)
    346340
    … …  
    365359              ! water ice sedimentation
    366360              call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
    367      &        ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rhocloud,
     361     &        ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rhocloud,
    368362     &        zqi(1,1,iq),wq,beta)
    369363            else
    370364              ! water ice sedimentation
    371365              call newsedim(ngrid,nlay,ngrid*nlay,1,
    372      &        ptimestep,pplev,masse,epaisseur,pt,rsedcloud,rho_q(iq),
     366     &        ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rho_q(iq),
    373367     &        zqi(1,1,iq),wq,beta)
    374368            endif ! of if (microphys)
    … …  
    383377          else
    384378            call newsedim(ngrid,nlay,1,1,ptimestep,
    385      &      pplev,masse,epaisseur,pt,radius(iq),rho_q(iq),
     379     &      pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),
    386380     &      zqi(1,1,iq),wq,1.0)
    387381c           Tendencies
    … …  
    405399c =================================================================
    406400      enddo ! of do iq=1,nq
    407  
     401
    408402c     Update the dust particle size "rdust"
    409403c     -------------------------------------
    … …  
    451445      endif ! of if (water)
    452446
    453       RETURN
    454447      END
    455448

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

    r740 r1005  
    1414c   by the dynamics in "advtrac.h"
    1515c
    16 c   Old conventions: (not used any more)
    17 c
    18 c   If water=T : q(iq=nqmx) is the water mass mixing ratio
    19 c     and q(iq=nqmx-1) is the ice mass mixing ratio
    20 
    21 c   If there is transported dust, it uses iq=1 to iq=dustbin
    22 c   If there is no transported dust : dustbin=0
    23 c   If doubleq=T : q(iq=1) is the dust mass mixing ratio
    24 c                  q(iq=2) is the dust number mixing ratio
    25 
    2616c
    2717c   author: F.Forget
    … …  
    380370
    381371c------------------------------------------------------------
    382 c     Initialisation tracers ....
     372c     Initialize tracers .... (in tracer.h)
    383373c------------------------------------------------------------
    384       call zerophys(nqmx,rho_q)
    385 
     374      ! start by setting everything to (default) zero
     375      rho_q(1:nqmx)=0     ! tracer density (kg.m-3)
     376      radius(1:nqmx)=0.   ! tracer particle radius (m)
     377      alpha_lift(1:nqmx) =0.  ! tracer saltation vertical flux/horiz flux ratio (m-1)
     378      alpha_devil(1:nqmx)=0.  ! tracer lifting coefficient by dust devils
     379
     380
     381      ! some reference values
    386382      rho_dust=2500.  ! Mars dust density (kg.m-3)
    387383      rho_ice=920.    ! Water ice density (kg.m-3)
    … …  
    497493      end if  ! (submicron)
    498494
    499 c     Initialization for photochemistry:
    500 c     ---------------------------------
    501       if (photochem) then
    502       ! initialize chemistry+water (water will be correctly initialized below)
    503       ! by initializing everything which is not dust ...
    504         do iq=1,nqmx
    505           txt=noms(iq)
    506           if (txt(1:4).ne."dust") then
    507             radius(iq)=0.
    508             alpha_lift(iq) =0.
    509             alpha_devil(iq)=0.
    510           endif
    511         enddo ! do iq=1,nqmx
    512       endif
    513 
    514495c     Initialization for water vapor
    515496c     ------------------------------

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

    r999 r1005  
    11721172
    11731173           call callsedim(ngrid,nlayer, ptimestep,
    1174      &                zplev,zzlev, zzlay, pt, rdust, rice,
     1174     &                zplev,zzlev, zzlay, pt, pdt, rdust, rice,
    11751175     &                rsedcloud,rhocloud,
    11761176     &                pq, pdq, zdqsed, zdqssed,nq,

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

    r410 r1005  
    9292      integer,save :: n_nom_def
    9393      integer :: n
    94       integer,parameter :: n_nom_def_max=99
    95       character(len=20),save :: nom_def(n_nom_def_max)
     94      integer,parameter :: n_nom_def_max=199
     95      character(len=120),save :: nom_def(n_nom_def_max)
    9696      logical,save :: firstcall=.true.
    9797     
Note: See TracChangeset for help on using the changeset viewer.