Changeset 1795


Ignore:
Timestamp:
Oct 12, 2017, 12:26:18 PM (7 years ago)
Author:
jvatant
Message:

Making Titan's hazy again, part II
+ Added calmufi and inimufi routines that interface YAMMS model
+ Major changes of the tracer gestion in tracer_h (new query by name)
+ Update the deftank
JVO

Location:
trunk/LMDZ.TITAN
Files:
6 added
3 deleted
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/README

    r1789 r1795  
    13351335
    13361336== 28/09/2017 == JVO
    1337 Added the surface methane tank and put it in start files
     1337Added the surface methane tank and put it in start files
     1338
     1339== 06/10/2017 and 12/10/2017 : r1792-94 == JVO
     1340Added the microphysical model YAMMS from J.Burgalat in libf/muphytitan
     1341Added the routines calmufi and inimufi that interfaces this model
     1342For now on you can do microphysics but no clouds yet ( you can but you shouldn't !)
     1343Big modifs of the tracer gestion/init in the physiq with a new query by names (see tracer_h )

  • trunk/LMDZ.TITAN/deftank/callphys.def

    r1790 r1795  
    3939corrk      = .true.
    4040# folder in which correlated-k data is stored ?
    41   corrkdir   = 23x23_nocia_z55_T220
     41  corrkdir   = 23x23
    4242# call visible gaseous absorption in radiative transfer ?
    4343callgasvis = .true.
    … …  
    7777# call chemistry ? (must have tracer=true)
    7878callchim = .false.
     79# Activate zonal mean ? (must be true if callchim)
     80moyzon_ch = .false.
    7981# chemistry is computed every "ichim" phys. timestep
    80 ichim = 20
    81 # Activate zonal mean ? (must match callchim)
    82 moyzon_ch = .false.
     82ichim = 200
     83
     84## Microphysics
     85## ~~~~~~~~~~~~
     86# call microphysics ? (must have tracer=true)
     87callmufi = .false.
     88# Activate zonal mean ? (must be true if callmufi by now)
     89moyzon_mu = .false.
     90## If yes, clouds ?
     91callclouds = .false.
     92### If yes, number of ices ? (must be compatible with traceur.def AND microphysical model )
     93nices = 3
     94# ~~~~~ Parameters for microphysics ~~~~~
     95# Path to microphys. config file ?
     96config_mufi = ../../datagcm/microphysics/config.ne15.cfg
     97# Fractal dimension ?
     98df_mufi = 2.0
     99# Monomer radius (m) ?
     100rm_mufi = 6.66e-08
     101# Aerosol density (kg.m-3) ?
     102rho_aer_mufi = 1.e3
     103# Pressure level of aer. production (Pa) ?
     104p_prod = 1.0
     105# Aer. production rate (kg.m-2.s-1) ?
     106tx_prod = 3.5e-13
     107# Equivalent radius production (m) ?
     108rc_prod = 2.0e-8
     109# Radius of air (nitrogen) molecule (m) ?
     110air_rad = 1.75e-10
     111
    83112
    84113## Star type

  • trunk/LMDZ.TITAN/libf/phytitan/calchim.F90

    r1787 r1795  
    164164       enddo
    165165       rmil(NLEV) = rinter(NLEV)+(rinter(NLEV)-rinter(NLEV-1))/2.
     166! rmil et rinter ne servent que pour la diffusion (> plafond-100km) donc ok en moyenne planetaire
    166167
    167168! lecture de tcp.ver, une seule fois

  • trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90

    r1788 r1795  
    44      logical,save :: callrad,corrk,calldifv,UseTurbDiff
    55!$OMP THREADPRIVATE(callrad,corrk,calldifv,UseTurbDiff)
    6       logical,save :: calladj,callsoil,callchim
    7 !$OMP THREADPRIVATE(calladj,callsoil,callchim)
     6      logical,save :: calladj,callsoil
     7!$OMP THREADPRIVATE(calladj,callsoil)
    88      logical,save :: season,diurnal,tlocked,rings_shadow,lwrite
    99!$OMP THREADPRIVATE(season,diurnal,tlocked,rings_shadow,lwrite)
    … …  
    1414      logical,save :: strictboundcorrk                                     
    1515!$OMP THREADPRIVATE(strictboundcorrk)
     16
     17      logical,save :: callchim, callmufi, callclouds
     18!$OMP THREADPRIVATE(callchim,callmufi,callclouds)
    1619
    1720      logical,save :: enertest
    … …  
    3639     
    3740      integer,save :: ichim
     41!$OMP THREADPRIVATE(ichim)
     42
    3843      integer,save :: iddist
    3944      integer,save :: iradia
    4045      integer,save :: startype
    41 !$OMP THREADPRIVATE(ichim,iddist,iradia,startype)
     46!$OMP THREADPRIVATE(iddist,iradia,startype)
     47     
     48      real,save :: df_mufi, rm_mufi, rho_aer_mufi
     49      real,save :: p_prod, tx_prod, rc_prod
     50      real,save :: air_rad
     51!$OMP THREADPRIVATE(df_mufi, rm_mufi, rho_aer_mufi,p_prod,tx_prod,rc_prod,air_rad)
    4252
    4353      real,save :: Fat1AU

  • trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90

    r1788 r1795  
    55
    66      ! Main directory: 'datadir':
    7       ! Default for Berserker @ UChicago:
    8 !      character(len=300) :: datadir='/home/rwordsworth/datagcm'
    9       ! Default for Gnome Idataplex:
    10 !      character(len=300) :: datadir='/san/home/rdword/gcm/datagcm'
    117      ! Default for LMD machines:
    12       character(len=300),save :: datadir='/u/lmdz/WWW/planets/LMDZ.GENERIC/datagcm'
     8      character(len=300),save :: datadir='datagcm'
    139!$OMP THREADPRIVATE(datadir)
    1410     
    … …  
    1814      character(len=12),parameter :: surfdir="surface_data"
    1915     
     16      ! Default directory for microphysics
     17      character(LEN=300),save :: config_mufi ='datagcm/microphysics/config.cfg'
     18
    2019      end module datafile_mod
    2120!-----------------------------------------------------------------------

  • trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90

    r1788 r1795  
    1010
    1111  use radinc_h, only: ini_radinc_h
    12   use datafile_mod, only: datadir
     12  use datafile_mod, only: datadir,config_mufi
    1313  use comdiurn_h, only: sinlat, coslat, sinlon, coslon
    1414  use comgeomfi_h, only: totarea, totarea_planet
    … …  
    319319     call getin_p("ichim",ichim)
    320320     write(*,*)" ichim = ",ichim
     321
     322! Microphysics
     323
     324     write(*,*) "Run with or without microphysics?"
     325     callmufi=.false. ! default value
     326     call getin_p("callmufi",callmufi)
     327     write(*,*)" callmufi = ",callmufi
     328
     329     ! sanity check
     330     if (callmufi.and.(.not.tracer)) then
     331       print*,"You are running microphysics without tracer"
     332       print*,"Please start again with tracer =.true."
     333       stop
     334     endif
     335
     336     write(*,*) "Compute clouds?"
     337     callclouds=.false. ! default value
     338     call getin_p("callclouds",callclouds)
     339     write(*,*)" callclouds = ",callclouds
     340
     341     ! sanity check
     342     if (callclouds.and.(.not.callmufi)) then
     343       print*,"You are trying to make clouds without microphysics !"
     344       print*,"Please start again with callmufi =.true."
     345       stop
     346     endif
     347
     348     write(*,*) "Fractal dimension ?"
     349     df_mufi=2.0 ! default value
     350     call getin_p("df_mufi",df_mufi)
     351
     352     write(*,*) "Monomer radius (m) ?"
     353     rm_mufi=6.66e-08 ! default value
     354     call getin_p("rm_mufi",rm_mufi)
     355
     356     write(*,*) "Aerosol density (kg.m-3)?"
     357     rho_aer_mufi=1.e3 ! default value
     358     call getin_p("rho_aer_mufi",rho_aer_mufi)
     359
     360     write(*,*) "Pressure level of aer. production (Pa) ?"
     361     p_prod=1.0 ! default value
     362     call getin_p("p_prod",p_prod)
     363     
     364     write(*,*) "Aerosol production rate (kg.m-2.s-1) ?"
     365     tx_prod=3.5e-13 ! default value
     366     call getin_p("tx_prod",tx_prod)
     367
     368     write(*,*) "Equivalent radius production (m) ?"
     369     rc_prod=2.0e-8 ! default value
     370     call getin_p("rc_prod",rc_prod)
     371
     372     write(*,*) "Radius of air (nitrogen) molecule (m) ?"
     373     air_rad=1.75e-10 ! default value
     374     call getin_p("air_rad",air_rad)
     375
     376     write(*,*) "Path to microphys. config file ?"
     377     config_mufi='datagcm/microphysics/config.cfg' ! default value
     378     call getin_p("config_mufi",config_mufi)
    321379
    322380! Soil model

  • trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90

    r1789 r1795  
    2222      use geometry_mod, only: latitude, longitude, cell_area
    2323      USE comgeomfi_h, only: totarea, totarea_planet
    24       USE tracer_h, only: noms, mmol
     24      USE tracer_h
    2525      use time_phylmdz_mod, only: ecritphy, iphysiq, nday
    2626      use phyetat0_mod, only: phyetat0
    … …  
    3232      use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
    3333      use time_phylmdz_mod, only: daysec
    34       use logic_mod, only: moyzon_ch
     34      use logic_mod, only: moyzon_ch, moyzon_mu
    3535      use moyzon_mod, only: tmoy, playmoy, zphibar, zphisbar, zplevbar, &
    3636                            zplaybar, zzlevbar, zzlaybar, ztfibar
    … …  
    7373!      V. Tracers
    7474!         V.1. Chemistry
     75!         V.2. Microphysics
    7576!         V.3. Updates (pressure variations, surface budget).
    7677!         V.4. Surface Tracer Update.
    … …  
    159160!   -------
    160161
     162
    161163      integer,intent(in) :: ngrid             ! Number of atmospheric columns.
    162164      integer,intent(in) :: nlayer            ! Number of atmospheric layers.
    … …  
    277279     
    278280      real zdqchi(ngrid,nlayer,nq)    ! Chemical tendency ( chemistry routine ).
    279       real zdqmph(ngrid,nlayer,nq)    ! Microphysical tendency ( condensation routine only for now, no microphysical routines ).
     281      real zdqcond(ngrid,nlayer,nq)   ! Condensation tendency ( chemistry routine ).
     282     
     283      real zdqmufi(ngrid,nlayer,nq)   ! Microphysical tendency.
    280284                       
    281285      ! For Winds : (m/s/s)
    … …  
    349353! ----------------------------------------------------------------------------
    350354     
    351       integer,parameter  :: nmicro=0 ! Temporary ! To be put in start/def
     355      character*10,dimension(:),allocatable,save :: nomqy
     356
     357!$OMP THREADPRIVATE(nomqy)
    352358
    353359      real ctimestep ! Chemistry timestep (s)
    … …  
    356362      real temp_eq(nlayer), press_eq(nlayer)
    357363      real zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
    358 !      real zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer)
    359364      real ztemp(ngrid,nlayer)
    360365     
    361       real ychim(ngrid,nlayer,nq-nmicro)   
    362 
    363       real rat_mmol(nq) ! Molar fraction ratio
    364 
    365       ! 2D vmr tendencies ( chemistry and condensation )
    366       real,dimension(:,:,:),allocatable,save :: dycchi
    367       ! Must be saved since chemistry is not called every step
    368      
    369       real dycmph(ngrid,nlayer,nq-nmicro)       
    370       ! ----------------------------------------------------------
    371          
    372       real,dimension(:,:),allocatable,save       :: qysat
    373      
    374       character*10,dimension(:),allocatable,save :: nomqy
    375 !$OMP THREADPRIVATE(dycchi,qysat,nomqy)
     366      real , allocatable, dimension(:,:,:),save :: ychim
     367
     368      ! 2D vmr tendencies ( chemistry and condensation ) for chem. tracers
     369      real,dimension(:,:,:),allocatable,save :: dycchi ! Saved since chemistry is not called every step
     370      real dyccond(ngrid,nlayer,nq)       
     371         
     372      real,dimension(:,:),allocatable,save   :: qysat
     373     
     374!$OMP THREADPRIVATE(ychim,dycchi,qysat)
    376375
    377376      real,dimension(:),allocatable,save :: tankCH4 ! Depth of surface methane tank
    … …  
    379378
    380379!-----------------------------------------------------------------------------
     380    ! Interface to calmufi
     381    !   --> needed in order to pass assumed-shape arrays. Otherwise we must put calmufi in a module
     382    !       (to have an explicit generated by the compiler).
     383    !   Or one can put calmufi in MMP_GCM module (in muphytitan).
     384    INTERFACE
     385      SUBROUTINE calmufi(plev, zlev, play, zlay, temp, pq, zdq)
     386        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: plev  !! Pressure levels (Pa).
     387        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlev  !! Altitude levels (m).
     388        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: play  !! Pressure layers (Pa).
     389        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: zlay  !! Altitude at the center of each layer (m).
     390        REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp  !! Temperature at the center of each layer (K). 
     391        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: pq   !! Tracers (\(kg.kg^{-1}}\)).
     392        REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq  !! Tendency for tracers (\(kg.kg^{-1}}\)).
     393      END SUBROUTINE calmufi
     394    END INTERFACE
    381395     
    382396!==================================================================================================
    … …  
    390404! --------------------------------
    391405      if (firstcall) then
     406        ! initracer:
     407        !   initialize nmicro,
     408        call initracer2(nq,nametrac)
     409
     410        ! ----------------------------------------------------------------------------
    392411
    393412        ! Allocate saved arrays.
    … …  
    423442        ALLOCATE(dycchi(ngrid,nlayer,nq-nmicro))
    424443        ALLOCATE(qysat(nlayer,nq-nmicro))
    425         ALLOCATE(nomqy(nq-nmicro+1))
     444        ALLOCATE(nomqy(nq-nmicro+1)) ! +1 because of hv
    426445       
    427446        ! This is defined in comsaison_h
    … …  
    439458         zdtlw(:,:) = 0.0
    440459
    441      
    442 !        Initialize tracer names, indexes and properties.
    443 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    444          IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq)) ! (because noms is an argument of physdem1 whether or not tracer is on)
    445          if (tracer) then
    446             call initracer(ngrid,nq,nametrac)
    447          endif
    448 
    449          rat_mmol(:) = mmol(:) / mugaz
     460!        Initialize names, timestep and saturation profiles for chemistry
     461!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     462
     463         if ( callchim .and. (nq.gt.nmicro) ) then
     464
     465            ctimestep = ptimestep*REAL(ichim)
     466
     467            do iq=nmicro+1,nq
     468               nomqy(iq-nmicro) = nametrac(iq)
     469            enddo
     470
     471            nomqy(nq-nmicro+1) = "HV"
     472           
     473            ! qysat is taken at the equator ( small variations of t,p)
     474            temp_eq(:)  = tmoy(:)
     475            press_eq(:) = playmoy(:)/100. ! en mbar
     476           
     477            call inicondens(nq-nmicro,press_eq,temp_eq,nomqy,qysat)
     478         
     479            zdqchi(:,:,:)  = 0.0
     480            zdqcond(:,:,:) = 0.0
     481
     482         endif
     483
     484!        Initialize microphysics.
     485!        ~~~~~~~~~~~~~~~~~~~~~~~~
     486
     487         IF ( callmufi ) THEN
     488
     489           call inimufi(nq,ptimestep)
     490
     491         ENDIF
     492
    450493
    451494!        Read 'startfi.nc' file.
    … …  
    485528         call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
    486529
    487 
    488530         if(tlocked)then
    489531            print*,'Planet is tidally locked at resonance n=',nres
    … …  
    533575         endif
    534576
     577         ! Sanity check for microphysics
     578         if ( ((.not.moyzon_mu).and.(callmufi)) ) then
     579            print *, "moyzon_mu=",moyzon_mu," and callmufi=",callmufi
     580            print *, "Please activate zonal mean to run microphysics (for now) !"
     581            stop
     582         endif
     583
    535584         ! Sanity check for chemistry
    536          if ( ((moyzon_ch).and.(.not.callchim)) .or. ((.not.moyzon_ch).and.(callchim)) ) then
     585         if ( (.not.moyzon_ch) .and. (callchim) ) then
    537586            print *, "moyzon_ch=",moyzon_ch," and callchim=",callchim
    538             print *, "This is not compatible..."
     587            print *, "Please activate zonal mean to run chemistry !"
    539588            stop
    540589         endif
    541590
    542          ! Initialize names, timestep and saturation profiles for chemistry
    543 
    544          if ( callchim .and. (nq.gt.nmicro) ) then
    545 
    546             ctimestep = ptimestep*REAL(ichim)
    547 
    548             do iq=nmicro+1,nq
    549                nomqy(iq-nmicro) = nametrac(iq)
    550             enddo
    551 
    552             nomqy(nq-nmicro+1) = "HV"
    553            
    554             ! qysat is taken at the equator ( small variations of t,p)
    555             temp_eq(:)  = tmoy(:)
    556             press_eq(:) = playmoy(:)/100. ! en mbar
    557            
    558             call inicondens(nq-nmicro,press_eq,temp_eq,nomqy,qysat)
    559          
    560          endif
    561591         
    562592         ! XIOS outputs
    … …  
    647677      ! -------------------------------Taken from old Titan --------------------------
    648678      ! zonal averages needed
    649       if (moyzon_ch) then
     679      if (moyzon_ch .or. moyzon_mu) then
    650680         
    651681         zzlaybar(1,:)=(zphibar(1,:)+zphisbar(1))/g
    … …  
    681711
    682712      endif  ! moyzon
     713
    683714      ! -------------------------------------------------------------------------------------
    684      
    685715      ! Compute potential temperature
    686716      ! Note : Potential temperature calculation may not be the same in physiq and dynamic...
    … …  
    9811011  !   V.1. Chemistry
    9821012  ! -------------------------
     1013
    9831014         if (callchim) then
    9841015
    985             ! Utilisation de la moyenne zonale dans calchim
     1016            ! Using zonal mean for calchim
    9861017            zplev(:,:) = zplevbar(:,:)
    9871018            zplay(:,:) = zplaybar(:,:)
    … …  
    9901021            ztemp(:,:) = ztfibar(:,:)
    9911022
     1023            allocate(ychim(ngrid,nlayer,nq-nmicro))
     1024
    9921025            if (nq.gt.nmicro) then
    9931026               do iq = nmicro+1,nq
    … …  
    9971030
    9981031            ! Condensation tendency after the transport
    999             do iq=1,nq-nmicro
     1032            do iq=nmicro+1,nq
    10001033               do l=1,nlayer
    10011034                  do ig=1,ngrid
    1002                      if ( ychim(ig,l,iq).gt.qysat(l,iq) ) then
    1003                         dycmph(ig,l,nmicro+iq)= ( -ychim(ig,l,iq)+qysat(l,iq) ) / ptimestep
     1035                     if ( ychim(ig,l,iq-nmicro).gt.qysat(l,iq-nmicro) ) then
     1036                        dyccond(ig,l,iq)= ( -ychim(ig,l,iq-nmicro)+qysat(l,iq-nmicro) ) / ptimestep
    10041037                     endif
    10051038                  enddo
    … …  
    10071040            enddo
    10081041
    1009             if( mod(icount-1,ichim).eq.0.) then
     1042            if ( callclouds ) dyccond(:,:,ices_indx) = 0.0 ! Condensation will be calculated in the cloud microphysics
     1043
     1044            if( mod(icount-1,ichim).eq.0. ) then
    10101045               
    1011                print *, "On passe dans la chimie..."
    1012                
     1046               print *, "We enter in the chemistry ..."
    10131047               call calchim(ngrid,nq-nmicro,ychim,nomqy,declin,zls,ctimestep, &
    10141048                           ztemp,zplay,zplev,zzlay,zzlev,dycchi,nlayer+70)
    … …  
    10181052            endif
    10191053           
     1054            ! TEMPORARY COMMENT
     1055            ! These conversions as well as 2D/1D should be put in phytrac
     1056            ! ( GCM-friendly tracers and tendencies -> format for photochem and microphys )
     1057
    10201058            if (nq.gt.nmicro) then
    1021                ! We convert tendencies back to mass mixing ratio
     1059               ! We convert tendencies to mass mixing ratio
    10221060               do iq=nmicro+1,nq
    1023                   zdqchi(:,:,iq) = dycchi(:,:,iq-nmicro) / rat_mmol(iq)
    1024                   zdqmph(:,:,iq) = dycmph(:,:,iq-nmicro) / rat_mmol(iq)
     1061                  zdqchi(:,:,iq)  = dycchi(:,:,iq-nmicro) / rat_mmol(iq)
     1062                  zdqcond(:,:,iq) = dyccond(:,:,iq) / rat_mmol(iq)
    10251063               enddo
    10261064
    10271065               pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + &
    1028                     zdqchi(1:ngrid,1:nlayer,1:nq) + zdqmph(1:ngrid,1:nlayer,1:nq)
     1066                    zdqchi(1:ngrid,1:nlayer,1:nq) + zdqcond(1:ngrid,1:nlayer,1:nq)
    10291067               
    10301068            endif
    10311069           
    1032          endif
     1070         endif ! end of 'callchim'
     1071
     1072  ! -------------------
     1073  !   V.2 Microphysics
     1074  ! -------------------
     1075 
     1076         if (callmufi) then
     1077
     1078            ! Using zonal mean for microphysics
     1079            zplev(:,:) = zplevbar(:,:)
     1080            zplay(:,:) = zplaybar(:,:)
     1081            zzlev(:,:) = zzlevbar(:,:)
     1082            zzlay(:,:) = zzlaybar(:,:)
     1083            ztemp(:,:) = ztfibar(:,:)
     1084
     1085            ! Inside this routine we will split 2D->1D, intensive->extensive and separate different types of tracers
     1086            ! Should be put in phytrac
     1087
     1088            call calmufi(zplev,zzlev,zplay,zzlay,ztemp,pq,zdqmufi)
     1089
     1090            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmufi(1:ngrid,1:nlayer,1:nq)   
     1091
     1092         endif ! end of 'callmufi'
    10331093     
    10341094  ! ---------------
    … …  
    13461406     
    13471407         !call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
    1348          !call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
    13491408         call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
    13501409         call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
    13511410         call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
    1352          call writediagfi(ngrid,"shad","rings"," ", 2, fract)
    13531411         
    13541412!           call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
    … …  
    14031461      if (tracer) then
    14041462     
    1405          do iq=1,nq
    1406             call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
    1407             call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
    1408                              'kg m^-2',2,qsurf_hist(1,iq) )
    1409                          
    1410 !          call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
    1411 !                          'kg m^-2',2,qsurf(1,iq) )                         
    1412 
    1413          enddo ! end of 'nq' loop
    1414          
     1463         if (callmufi) then ! For now we assume an given order for tracers !
     1464            call writediagfi(ngrid,"mu_m0as","Spherical mode 0th order moment",'kg/kg',3,zq(:,:,1))
     1465            call writediagfi(ngrid,"mu_m3as","Spherical mode 3rd order moment",'kg/kg',3,zq(:,:,2))
     1466            call writediagfi(ngrid,"mu_m0af","Fractal mode 0th order moment",'kg/kg',3,zq(:,:,3))
     1467            call writediagfi(ngrid,"mu_m3af","Fractal mode 3rd order moment",'kg/kg',3,zq(:,:,4))
     1468         endif ! end of 'callmufi'
     1469
    14151470       endif ! end of 'tracer'
    1416 
    1417 
    1418       ! Output spectrum.
    1419       if(specOLR.and.corrk)then     
    1420          call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W/m^2/cm^-1",3,OLR_nu)
    1421          call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W/m^2/cm^-1",3,OSR_nu)
    1422       endif
    14231471
    14241472! XIOS outputs

  • trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90

    r1788 r1795  
    11
    2        module tracer_h
    3 
    4        implicit none
    5 
    6 ! nqtot : total number of tracers
    7        INTEGER, SAVE :: nqtot
    8 !$OMP THREADPRIVATE(nqtot)
    9 
    10        character*20, DIMENSION(:), ALLOCATABLE :: noms   ! name of the tracer
    11        real, DIMENSION(:), ALLOCATABLE :: mmol     ! mole mass of tracer (g/mol-1)
    12        real, DIMENSION(:), ALLOCATABLE :: rho_q    ! tracer densities (kg.m-3)
    13 
    14 !$OMP THREADPRIVATE(noms,mmol,rho_q)
    15 
    16 ! tracer indexes: these are initialized in initracer and should be 0 if the
    17 !                 corresponding tracer does not exist
    18 
    19       ! chemistry:
    20      
    21       integer :: igcm_h
    22       integer :: igcm_h2
    23       integer :: igcm_ch
    24       integer :: igcm_ch2s
    25       integer :: igcm_ch2
    26       integer :: igcm_ch3
    27       integer :: igcm_ch4
    28       integer :: igcm_c2
    29       integer :: igcm_c2h
    30       integer :: igcm_c2h2
    31       integer :: igcm_c2h3
    32       integer :: igcm_c2h4
    33       integer :: igcm_c2h5
    34       integer :: igcm_c2h6
    35       integer :: igcm_c3h3
    36       integer :: igcm_c3h5
    37       integer :: igcm_c3h6
    38       integer :: igcm_c3h7
    39       integer :: igcm_c4h
    40       integer :: igcm_c4h3
    41       integer :: igcm_c4h4
    42       integer :: igcm_c4h2s
    43       integer :: igcm_ch2cch2
    44       integer :: igcm_ch3cch
    45       integer :: igcm_c3h8
    46       integer :: igcm_c4h2
    47       integer :: igcm_c4h6
    48       integer :: igcm_c4h10
    49       integer :: igcm_ac6h6
    50       integer :: igcm_c3h2
    51       integer :: igcm_c4h5
    52       integer :: igcm_ac6h5
    53       integer :: igcm_n2
    54       integer :: igcm_n4s
    55       integer :: igcm_cn
    56       integer :: igcm_hcn
    57       integer :: igcm_h2cn
    58       integer :: igcm_chcn
    59       integer :: igcm_ch2cn
    60       integer :: igcm_ch3cn
    61       integer :: igcm_c3n
    62       integer :: igcm_hc3n
    63       integer :: igcm_nccn
    64       integer :: igcm_c4n2
    65      
    66      
    67 !$OMP THREADPRIVATE(igcm_h,igcm_h2,igcm_ch,igcm_ch2s,igcm_ch2,igcm_ch3,igcm_ch4,   &
    68       !$OMP igcm_c2,igcm_c2h,igcm_c2h2,igcm_c2h3,igcm_c2h4,igcm_c2h5,igcm_c2h6,    &
    69       !$OMP igcm_c3h3,igcm_c3h5,igcm_c3h6,igcm_c3h7,igcm_c4h,igcm_c4h3,igcm_c4h4,  &
    70       !$OMP igcm_c4h2s,igcm_ch2cch2,igcm_ch3cch,igcm_c3h8,igcm_c4h2,igcm_c4h6,     &
    71       !$OMP igcm_c4h10,igcm_ac6h6,igcm_c3h2,igcm_c4h5,igcm_ac6h5,igcm_n2,igcm_n4s, &
    72       !$OMP igcm_cn,igcm_hcn,igcm_h2cn,igcm_chcn,igcm_ch2cn,igcm_ch3cn,igcm_c3n,   &
    73       !$OMP igcm_hc3n,igcm_nccn,igcm_c4n2)
    74 
    75        end module tracer_h
    76 
     2MODULE tracer_h
     3  !! Stores data related to physics tracers.
     4  !!
     5  !! The module stores public global variables related to the number of tracers
     6  !! available in the physics and their kind:
     7  !!
     8  !! Currently, tracers can be used either for chemistry process (nchimi) or
     9  !! microphysics (nmicro).
     10  !!
     11  !! The subroutine "initracer2" initializes and performs sanity check of
     12  !! the tracer definitions given in traceur.def and the required tracers in physics
     13  !! (based on the run parameters).
     14  !!
     15  !! The module provides additional methods:
     16  !!
     17  !!   - indexOfTracer : search for the index of a tracer in the global table (tracers_h:noms) by name.
     18  !!   - nameOfTracer  : get the name of tracer from a given index (of the global table).
     19  !!   - dumpTracers   : print the names of all tracers indexes given in argument.
     20  !!
     21  IMPLICIT NONE
     22
     23  INTEGER, SAVE :: nqtot  = 0 !! Total number of tracers
     24  INTEGER, SAVE :: nmicro = 0 !! Number of microphysics tracers.
     25  INTEGER, SAVE :: nice   = 0 !! Number of microphysics ice tracers (subset of nmicro).
     26  INTEGER, SAVE :: nchimi = 0 !! Number of chemical (gaz species) tracers.
     27  !$OMP THREADPRIVATE(nqtot,nmicro,nice,nchimi)
     28
     29  INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: chimi_indx  !! Indexes of all chemical species tracers
     30  INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: micro_indx  !! Indexes of all microphysical tracers
     31  INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: ices_indx   !! Indexes of all ice microphysical tracers
     32
     33  CHARACTER(len=20), DIMENSION(:), ALLOCATABLE, SAVE :: noms      !! name of the tracer
     34  REAL, DIMENSION(:), ALLOCATABLE, SAVE              :: mmol      !! mole mass of tracer(g/mol-1)
     35  REAL, DIMENSION(:), ALLOCATABLE, SAVE              :: rat_mmol  !! molar mass ratio
     36  REAL, DIMENSION(:), ALLOCATABLE, SAVE              :: rho_q     !! tracer densities (kg.m-3)
     37  !$OMP THREADPRIVATE(noms,mmol,rat_mmol,rho_q)
     38
     39
     40
     41  ! tracer indexes: these are initialized in initracer and should be 0 if the
     42  !                 corresponding tracer does not exist
     43
     44
     45  CONTAINS
     46
     47  SUBROUTINE initracer2(nq,nametrac,talk_to_me)
     48    !! Initialize tracer names and attributes.
     49    !!
     50    !! The method initializes the list of tracer names used in the physics from their
     51    !! dynamics counterpart.
     52    !!
     53    !! In addition, it initializes arrays of indexes for the different sub-processs of the physics:
     54    !!
     55    !!   - tracers_h:micro_indxs, the array of tracers indexes used for the microphysics.
     56    !!   - tracers_h:chimi_indxs, the array of tracers indexes used for the chemistry.
     57    !!
     58    !! The method also initializes the molar mass array (tracers_h:mmol) for the chemistry and the
     59    !! molar mass ratio (tracers_h:rat_mmol).
     60    !!
     61    !! @note
     62    !! Strict checking of chemical species name is performed here if the chemistry is activated
     63    !! (see callchim variable). All the values of 'cnames' must be found in the tracers names
     64    !! related to chemistry.
     65    !! @note
     66    !! Tests are more permissive for the microphysics and is only based on the mimimum number of
     67    !! tracers expected. Strict name checking is performed in inimufi.
     68    USE callkeys_mod
     69    USE comcstfi_mod, only: mugaz
     70    IMPLICIT NONE
     71
     72    INTEGER, INTENT(in)                          :: nq         !! Total number of tracers (fixed at compile time)
     73    character(len=20), DIMENSION(nq), INTENT(in) :: nametrac   !! name of the tracer from dynamics (from 'traceurs.def')
     74    LOGICAL, INTENT(in), OPTIONAL                :: talk_to_me !! Enable verbose mode.
     75
     76    LOGICAL                                      :: verb,found
     77    CHARACTER(len=20)                            :: str
     78    !! Hard-coded chemical species for Titan chemistry
     79    CHARACTER(len=10), DIMENSION(44), PARAMETER  :: cnames = &
     80      (/"H         ", "H2        ", "CH        ", "CH2s      ", "CH2       ", "CH3       ", &
     81        "CH4       ", "C2        ", "C2H       ", "C2H2      ", "C2H3      ", "C2H4      ", &
     82        "C2H5      ", "C2H6      ", "C3H3      ", "C3H5      ", "C3H6      ", "C3H7      ", &
     83        "C4H       ", "C4H3      ", "C4H4      ", "C4H2s     ", "CH2CCH2   ", "CH3CCH    ", &
     84        "C3H8      ", "C4H2      ", "C4H6      ", "C4H10     ", "AC6H6     ", "C3H2      ", &
     85        "C4H5      ", "AC6H5     ", "N2        ", "N4S       ", "CN        ", "HCN       ", &
     86        "H2CN      ", "CHCN      ", "CH2CN     ", "CH3CN     ", "C3N       ", "HC3N      ", &
     87        "NCCN      ", "C4N2      "/)
     88    !! Hard-coded chemical species molar mass (g.mol-1), shares the same indexing than cnames.
     89    REAL, DIMENSION(44), PARAMETER               :: cmmol = (/ &
     90        1.01   , 2.0158, 13.02, 14.03, 14.03, 15.03, 16.04  , 24.02, 25.03, 26.04  , 27.05  , &
     91        28.05  , 29.06 , 30.07, 39.06, 41.07, 42.08, 43.09  , 49.05, 51.07, 52.08  , 50.06  , &
     92        40.07  , 40.07 , 44.11, 50.06, 54.09, 58.13, 78.1136, 38.05, 53.07, 77.1136, 28.0134, &
     93        14.01  , 26.02 , 27.04, 28.05, 39.05, 40.04, 41.05  , 50.04, 51.05, 52.04  , 76.1   /)
     94
     95    INTEGER :: i,j,n
     96
     97    verb = .true. ; IF (PRESENT(talk_to_me)) verb = talk_to_me
     98
     99    ! nqtot could be used everywhere in the physic :)
     100    nqtot=nq
     101
     102    IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq))
     103    noms(:)=nametrac(:)
     104
     105    ALLOCATE(rho_q(nq)) ! Defined for all tracers, currently initialized to 0.0
     106    rho_q(:) = 0.0
     107
     108    ALLOCATE(mmol(nq),rat_mmol(nq))  ! Defined for all tracers, (actually) initialized only for chemical tracers
     109    mmol(:)  = 0.0
     110    rat_mmol(:) = 1.0
     111
     112    ! Compute number of microphysics tracers:
     113    ! By convention they all have the prefix "mu_" (case sensitive !)
     114    nmicro = 0
     115    IF (callmufi) THEN
     116      DO i=1,nq
     117        str = noms(i)
     118        IF (str(1:3) == "mu_") nmicro = nmicro+1
     119      ENDDO
     120      ! Checking the expected number of tracers:
     121      !   no cloud:  4 ; w cloud :  4 + 2 + (1+)
     122      ! Note that we do not make assumptions on the number of chemical species for clouds, this
     123      ! will be checked in inimufi.
     124      IF (callclouds) THEN
     125        IF (nmicro < 7) THEN
     126          WRITE(*,'((a),I3,(a))') "initracer2:error: Inconsistent number of microphysical tracers &
     127            &(expected at least 7 tracers,",nmicro," given)"
     128          CALL abort_gcm("initracer2", "inconsistent number of tracers", 42)
     129          STOP
     130        ENDIF
     131      ELSE IF (nmicro < 4) THEN
     132          WRITE(*,'((a),I3,(a))') "initracer2:error: Inconsistent number of microphysical tracers &
     133            &(expected at least 4 tracers,",nmicro," given)"
     134          CALL abort_gcm("initracer2", "inconsistent number of tracers", 42)
     135      ELSE IF  (nmicro > 4) THEN
     136        WRITE(*,'(a)') "initracer2:info: I was expecting only four tracers, you gave me &
     137          &more. I'll just pretend nothing happen !"
     138      ENDIF
     139      ! microphysics indexes share the same values than original tracname.
     140      ALLOCATE(micro_indx(nmicro))
     141      j = 1
     142      DO i=1,nq
     143        str = noms(i)
     144        IF (str(1:3) == "mu_") THEN
     145          micro_indx(j) = i
     146          j=j+1
     147        ENDIF
     148      ENDDO
     149    ELSE
     150      ALLOCATE(micro_indx(nmicro))
     151    ENDIF
     152
     153    ! Compute number of chemical species:
     154    !   simply assume that all other tracers ARE chemical species
     155    nchimi = nqtot - nmicro
     156
     157    ! Titan chemistry requires exactly 44 tracers:
     158    ! Test should be in callchim condition
     159
     160    IF (callchim) THEN
     161      IF (nchimi < SIZE(cnames)) THEN
     162        WRITE(*,*) "initracer2:error: Inconsistent number of chemical species given (",SIZE(cnames)," expected)"
     163        CALL abort_gcm("initracer2", "inconsistent number of tracers", 42)
     164      ENDIF
     165      ALLOCATE(chimi_indx(nchimi))
     166      n = 0 ! counter on chimi_indx
     167      DO j=1,SIZE(cnames)
     168        found = .false.
     169        DO i=1,nq
     170          IF (TRIM(cnames(j)) == TRIM(noms(i))) THEN
     171            n = n + 1
     172            chimi_indx(n) = i
     173            mmol(i) = cmmol(j)
     174            rat_mmol(i) = cmmol(j)/mugaz
     175            found = .true.
     176            EXIT
     177          ENDIF
     178        ENDDO
     179        IF (.NOT.found) THEN
     180          WRITE(*,*) "initracer2:error: "//TRIM(cnames(j))//" is missing from tracers definition file."
     181        ENDIF
     182      ENDDO
     183      IF (n < SIZE(cnames)) THEN
     184        WRITE(*,*) "initracer2:error: Inconsistent number of chemical species given (",SIZE(cnames)," expected)"
     185        CALL abort_gcm("initracer2", "inconsistent number of tracers", 42)
     186      ENDIF
     187    ELSE
     188      ALLOCATE(chimi_indx(0))
     189    ENDIF
     190    IF (verb) THEN
     191      IF (callmufi.OR.callchim) WRITE(*,*) "===== INITRACER2 SPEAKING ====="
     192      IF (callmufi) THEN
     193         WRITE(*,*) "Found ",nmicro, "microphysical tracers"
     194        call dumpTracers(micro_indx)
     195        WRITE(*,*) "-------------------------------"
     196      ENDIF
     197      IF (callchim) THEN
     198        WRITE(*,*) "Found ",nchimi, "chemical tracers"
     199        call dumpTracers(chimi_indx)
     200        WRITE(*,*) "-------------------------------"
     201      ENDIF
     202    ENDIF
     203
     204  END SUBROUTINE initracer2
     205
     206  FUNCTION indexOfTracer(name, sensitivity) RESULT(idx)
     207    !! Get the index of a tracer by name.
     208    !!
     209    !! The function searches in the global tracer table (tracer_h:noms)
     210    !! for the given name and returns the first index matching "name".
     211    !!
     212    !! If no name in the table matches the given one, -1 is returned !
     213    !!
     214    !! @warning
     215    !! initracers must be called before any use of this function.
     216    IMPLICIT NONE
     217    CHARACTER(len=*), INTENT(in)  :: name         !! Name of the tracer to search.
     218    LOGICAL, OPTIONAL, INTENT(in) :: sensitivity  !! Case sensitivity (true by default).
     219    INTEGER :: idx                                !! Index of the first tracer matching name or -1 if not found.
     220    LOGICAL                  :: zsens
     221    INTEGER                  :: j
     222    CHARACTER(len=LEN(name)) :: zname
     223    zsens = .true. ; IF(PRESENT(sensitivity)) zsens = sensitivity
     224    idx = -1
     225    IF (.NOT.ALLOCATED(noms)) RETURN
     226    IF (zsens) THEN
     227      DO j=1,SIZE(noms)
     228        IF (TRIM(noms(j)) == TRIM(name)) THEN
     229          idx = j ; RETURN
     230        ENDIF
     231      ENDDO
     232    ELSE
     233      zname = to_lower(name)
     234      DO j=1,SIZE(noms)
     235        IF (TRIM(to_lower(noms(j))) == TRIM(zname)) THEN
     236          idx = j ; RETURN
     237        ENDIF
     238      ENDDO
     239    ENDIF
     240
     241    CONTAINS
     242
     243    FUNCTION to_lower(istr) RESULT(ostr)
     244      !! Lower case conversion function.
     245      IMPLICIT NONE
     246      CHARACTER(len=*), INTENT(in) :: istr
     247      CHARACTER(len=LEN(istr)) :: ostr
     248      INTEGER :: i, ic
     249      ostr = istr
     250      DO i = 1, LEN_TRIM(istr)
     251        ic = ICHAR(istr(i:i))
     252        IF (ic >= 65 .AND. ic < 90) ostr(i:i) = char(ic + 32)
     253      ENDDO
     254    END FUNCTION to_lower
     255  END FUNCTION indexOfTracer
     256
     257  FUNCTION nameOfTracer(indx) RESULT(name)
     258    !! Get the name of a tracer by index.
     259    !!
     260    !! The function searches in the global tracer table (tracer_h:noms)
     261    !! and returns the name of the tracer at given index.
     262    !!
     263    !! If the index is out of range an empty string is returned.
     264    !!
     265    !! @warning
     266    !! initracers must be called before any use of this function.
     267    IMPLICIT NONE
     268    INTEGER, INTENT(in) :: indx   !! Index of the tracer name to retrieve.
     269    CHARACTER(len=20)   :: name   !! Name of the tracer at given index.
     270    name = ''
     271    IF (.NOT.ALLOCATED(noms)) RETURN
     272    IF (indx <= 0 .OR. indx > SIZE(noms)) RETURN
     273    name = noms(indx)
     274  END FUNCTION nameOfTracer
     275
     276  SUBROUTINE dumpTracers(indexes)
     277    !! Print the names of the given list of tracers indexes.
     278    INTEGER, DIMENSION(:), INTENT(in) :: indexes
     279    INTEGER :: i,idx,nt
     280    CHARACTER(len=:), ALLOCATABLE :: suffix
     281   
     282    IF (.NOT.ALLOCATED(noms)) THEN
     283      WRITE(*,'(a)') "[tracers_h:dump_tracers] warning: 'noms' is not allocated, tracers_h:initracer2 has not be called yet"
     284      RETURN
     285    ENDIF
     286    nt = size(noms)
     287    WRITE(*,"(a)") "local -> global : name"
     288    DO i=1,size(indexes)
     289      idx = indexes(i)
     290      IF (idx < 1 .OR. idx >= nt) THEN
     291        ! WRITE(*,'((a),I3,(a),I3,(a))') "index out of range (",idx,"/",nt,")"
     292        CYCLE
     293      ENDIF
     294      IF (ANY(chimi_indx == idx)) THEN
     295        suffix = ' (chimi)'
     296      ELSE IF (ANY(micro_indx == idx)) THEN
     297        suffix = ' (micro'
     298        IF (ALLOCATED(ices_indx)) THEN
     299          IF (ANY(ices_indx == idx)) suffix=suffix//", ice"
     300        ENDIF
     301        suffix=suffix//")"
     302      ELSE
     303        suffix=" ()"
     304      ENDIF
     305      WRITE(*,'(I5,(a),I6,(a))') i," -> ",idx ," : "//TRIM(noms(i))//suffix
     306    ENDDO
     307  END SUBROUTINE dumpTracers
     308
     309
     310END MODULE tracer_h
     311
Note: See TracChangeset for help on using the changeset viewer.