Ignore:
Timestamp:
Jun 21, 2022, 11:16:11 AM (3 years ago)
Author:
aslmd
Message:

step 1 of moving initracer.F to initracer.F90: get all changes from initracer.F90 in recent commits to initracer.F

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/libf/phystd/initracer.F

    r2550 r2729  
    66      USE recombin_corrk_mod, ONLY: ini_recombin
    77      USE mod_phys_lmdz_para, only: is_master, bcast
     8      use generic_cloud_common_h
    89      IMPLICIT NONE
    9 c=======================================================================
    10 c   subject:
    11 c   --------
    12 c   Initialization related to tracer
    13 c   (transported dust, water, chemical species, ice...)
    14 c
    15 c   Name of the tracer
    16 c
    17 c   Test of dimension :
    18 c   Initialize COMMON tracer in tracer.h, using tracer names provided
    19 c   by the argument nametrac
    20 c
    21 c   author: F.Forget
    22 c   ------
    23 c            Ehouarn Millour (oct. 2008) identify tracers by their names
    24 c            Y Jaziri & J. Vatant d'Ollone (2020) : Modern traceur.def
    25 c=======================================================================
     10!=======================================================================
     11!   subject:
     12!   --------
     13!   Initialization related to tracer
     14!   (transported dust, water, chemical species, ice...)
     15!
     16!   Name of the tracer
     17!
     18!   Test of dimension :
     19!   Initialize COMMON tracer in tracer.h, using tracer names provided
     20!   by the argument nametrac
     21!
     22!   author: F.Forget
     23!   ------
     24!            Ehouarn Millour (oct. 2008) identify tracers by their names
     25!            Y Jaziri & J. Vatant d'Ollone (2020) : Modern traceur.def
     26!=======================================================================
    2627
    2728      integer,intent(in) :: ngrid,nq
    … …  
    2930
    3031      character(len=500) :: tracline   ! to read traceur.def lines
    31       character(len=30)  :: txt        ! to store some text
     32      ! character(len=30)  :: txt        ! to store some text
    3233      integer iq,ig,count,ierr
    3334
    34 
    35 
    36 c-----------------------------------------------------------------------
    37 c  radius(nq)      ! aerosol particle radius (m)
    38 c  rho_q(nq)       ! tracer densities (kg.m-3)
    39 c  qext(nq)        ! Single Scat. Extinction coeff at 0.67 um
    40 c  alpha_lift(nq)  ! saltation vertical flux/horiz flux ratio (m-1)
    41 c  alpha_devil(nq) ! lifting coeeficient by dust devil
    42 c  rho_dust          ! Mars dust density
    43 c  rho_ice           ! Water ice density
    44 c  doubleq           ! if method with mass (iq=1) and number(iq=2) mixing ratio
    45 c  varian            ! Characteristic variance of log-normal distribution
    46 c-----------------------------------------------------------------------
     35!-----------------------------------------------------------------------
     36!  radius(nq)      ! aerosol particle radius (m)
     37!  rho_q(nq)       ! tracer densities (kg.m-3)
     38!  qext(nq)        ! Single Scat. Extinction coeff at 0.67 um
     39!  alpha_lift(nq)  ! saltation vertical flux/horiz flux ratio (m-1)
     40!  alpha_devil(nq) ! lifting coeeficient by dust devil
     41!  rho_dust          ! Mars dust density
     42!  rho_ice           ! Water ice density
     43!  doubleq           ! if method with mass (iq=1) and number(iq=2) mixing ratio
     44!  varian            ! Characteristic variance of log-normal distribution
     45!-----------------------------------------------------------------------
    4746
    4847      if (is_master) then ! only the master proc/thread needs do this
    … …  
    5049        moderntracdef=.false. ! For modern traceur.def (default false, old type)
    5150
    52         open(407, form = 'formatted', status = 'old',
    53      $          file = 'traceur.def', iostat=ierr)
     51        open(407, form = 'formatted', status = 'old', &
     52             file = 'traceur.def', iostat=ierr)
    5453        if (ierr /=0) then
    55           call abort_physic('initracer',
    56      $    'Problem in opening traceur.def',1)
     54          ! call abort_physic('initracer',&
     55          ! 'Problem in opening traceur.def',1)
     56          print*,'Problem in opening traceur.def'
     57          stop
    5758        end if
    5859!! - Modif. by JVO and YJ for modern planetary traceur.def ---------------
    … …  
    6970               ! Temporary not implemented solution
    7071               if (nqtot/=nq) then
    71                  call abort_physic('initracer','Different number of '//
    72      & 'tracers in dynamics and physics not managed yet',1)
     72      !            call abort_physic('initracer','Different number of '// &
     73      ! 'tracers in dynamics and physics not managed yet',1)
     74                print*,'!= nbr oftracers in dynamics and physics not managed yet'
     75                stop
    7376               endif
    7477               EXIT
    7578             ENDIF
    7679           ELSE ! If pb, or if reached EOF without having found nqtot
    77              call abort_physic('initracer','Unable to read numbers '//
    78      & 'of tracers in traceur.def',1)
     80      !        call abort_physic('initracer','Unable to read numbers '// &
     81      !  'of tracers in traceur.def',1)
     82            print*,"unable to read numbers of tracer in tracer.def"
     83            stop
    7984           ENDIF
    8085         ENDDO
    … …  
    112117         ALLOCATE(is_recomb_qotf(nqtot))
    113118       ENDIF
     119       IF (.NOT. allocated(is_condensable)) allocate(is_condensable(nq)) !LT
     120       IF (.NOT. allocated(constants_mass)) allocate(constants_mass(nq))
     121       IF (.NOT. allocated(constants_delta_vapH)) allocate(constants_delta_vapH(nq))
     122       IF (.NOT. allocated(constants_Tref)) allocate(constants_Tref(nq))
     123       IF (.NOT. allocated(constants_Pref)) allocate(constants_Pref(nq))
     124       IF (.NOT. allocated(constants_epsi_generic)) allocate(constants_epsi_generic(nq))
     125       IF (.NOT. allocated(constants_RLVTT_generic)) allocate(constants_RLVTT_generic(nq))
     126       IF (.NOT. allocated(constants_metallicity_coeff)) allocate(constants_metallicity_coeff(nq))
     127
    114128       !! initialization
    115129       alpha_lift(:)     = 0.
    … …  
    128142       radius(:)=0.
    129143       qext(:)=0.
     144
     145       ! For condensable tracers, by Lucas Teinturier and Noé Clément (2022)
     146
     147       is_condensable(:)= 0
     148
     149       constants_mass(:)=0
     150       constants_delta_vapH(:)=0
     151       constants_Tref(:)=0
     152       constants_Pref(:)=0
     153       constants_epsi_generic(:)=0
     154       constants_RLVTT_generic(:)=0
     155       constants_metallicity_coeff(:)=0
     156
     157       rho_q(:) = 0. !need to be init here if we want to read it from modern traceur with get_tracdat
    130158
    131159
    … …  
    186214      igcm_hcaer=0
    187215
    188 
    189 
    190216      ! 1. find dust tracers
    191217      count=0
    … …  
    387413          count=count+1
    388414        endif
     415
     416      enddo ! of do iq=1,nq
     417
     418      ! 3. find condensable traceurs different from h2o and co2
     419      do iq=1,nq
     420        if ((index(noms(iq),"vap") .ne. 0) .and. (index(noms(iq),"h2o") .eq. 0) .and. (index(noms(iq),"co2") .eq. 0)) then
     421          count=count+1
     422        endif
     423        if ((index(noms(iq),"ice") .ne. 0) .and. (index(noms(iq),"h2o") .eq. 0) .and. (index(noms(iq),"co2") .eq. 0)) then
     424          count=count+1
     425        endif
     426
    389427      enddo ! of do iq=1,nq
    390428     
    … …  
    413451      if (is_master) close(407)
    414452
     453      ! Get specific data of condensable tracers
     454      do iq=1,nq
     455        if((is_condensable(iq)==1)) then
     456                write(*,*) "There is a specie which is condensable, for generic condensation : ", noms(iq)
     457                write(*,*) 'looking specie parameters for : ',noms(iq)(1:len(trim(noms(iq)))-4)
     458                call specie_parameters_table(noms(iq)(1:len(trim(noms(iq)))-4))
     459                constants_mass(iq)=m
     460                constants_delta_vapH(iq)=delta_vapH
     461                constants_Tref(iq)=Tref
     462                constants_Pref(iq)=Pref
     463                constants_epsi_generic(iq)=epsi_generic
     464                constants_RLVTT_generic(iq)=RLVTT_generic
     465                constants_metallicity_coeff(iq)=metallicity_coeff
     466        else
     467                write(*,*) "This tracer is not condensable, for generic condensation :  : ", noms(iq)
     468                write(*,*) "We keep condensable constants at zero"
     469        endif !(is_condensable(iq)==1) .and. (index(noms(iq),"vap") .ne. 0))
     470      enddo ! iq=1,nq
     471
    415472      ! Calculate number of species in the chemistry
    416473      nesp = sum(is_chim)
    417474      write(*,*) 'Number of species in the chemistry nesp = ',nesp
    418475
     476      ! Calculate number of generic tracers
     477      ngt = sum(is_condensable)
     478      write(*,*) 'Number of generic tracer is  ngt = ',ngt
    419479!     Processing modern traceur options
    420480      if(moderntracdef) then
    … …  
    422482      endif
    423483
    424 c------------------------------------------------------------
    425 c     Initialisation tracers ....
    426 c------------------------------------------------------------
    427       rho_q(1:nq)=0
     484!------------------------------------------------------------
     485!     Initialisation tracers ....
     486!------------------------------------------------------------
     487      ! rho_q(1:nq)=0
    428488
    429489      rho_dust=2500.  ! Mars dust density (kg.m-3)
    … …  
    431491      rho_co2=1620.   ! CO2 ice density (kg.m-3)
    432492
    433 c     Initialization for water vapor
    434 c     ------------------------------
     493!     Initialization for water vapor
     494!     ------------------------------
    435495      if(water) then
    436          radius(igcm_h2o_vap)=0.
    437          Qext(igcm_h2o_vap)=0.
    438          alpha_lift(igcm_h2o_vap) =0.
    439          alpha_devil(igcm_h2o_vap)=0.
    440         qextrhor(igcm_h2o_vap)= 0.
     496        radius(igcm_h2o_vap)=0.
     497        Qext(igcm_h2o_vap)=0.
     498        alpha_lift(igcm_h2o_vap) =0.
     499        alpha_devil(igcm_h2o_vap)=0.
     500        qextrhor(igcm_h2o_vap)= 0.
    441501
    442502         !! this is defined in surfdat_h.F90
    … …  
    455515!           alpha_lift(igcm_h2o_ice) =0.
    456516!           alpha_devil(igcm_h2o_ice)=0.
    457            qextrhor(igcm_h2o_ice)= (3./4.)*Qext(igcm_h2o_ice)
    458      $       / (rho_ice*radius(igcm_h2o_ice))
     517           qextrhor(igcm_h2o_ice)= (3./4.)*Qext(igcm_h2o_ice) &
     518            / (rho_ice*radius(igcm_h2o_ice))
    459519
    460520
    … …  
    469529         if (igcm_h2o_vap.eq.0) then
    470530           write(*,*) "initracer: error !!"
    471            write(*,*) "  cannot use water option without ",
    472      &                "an h2o_vap tracer !"
     531           write(*,*) "  cannot use water option without ",&
     532                     "an h2o_vap tracer !"
    473533           stop
    474534         endif
    475535         if (igcm_h2o_ice.eq.0) then
    476536           write(*,*) "initracer: error !!"
    477            write(*,*) "  cannot use water option without ",
    478      &                "an h2o_ice tracer !"
     537           write(*,*) "  cannot use water option without ",&
     538                     "an h2o_ice tracer !"
    479539           stop
    480540         endif
    … …  
    482542
    483543
    484 c     Output for records:
    485 c     ~~~~~~~~~~~~~~~~~~
     544!     Output for records:
     545!     ~~~~~~~~~~~~~~~~~~
    486546      write(*,*)
    487547      Write(*,*) '******** initracer : dust transport parameters :'
    … …  
    511571          ! option mmol
    512572          if (index(tracline,'mmol='   ) /= 0) then
    513               read(tracline(index(tracline,'mmol=')+len('mmol='):),*)
    514      $             mmol(iq)
    515               write(*,*) ' Parameter value (traceur.def) : mmol=',
    516      $             mmol(iq)
    517           else
    518               write(*,*) ' Parameter value (default)     : mmol=',
    519      $             mmol(iq)
     573              read(tracline(index(tracline,'mmol=')+len('mmol='):),*)&
     574                  mmol(iq)
     575              write(*,*) ' Parameter value (traceur.def) : mmol=', &
     576                  mmol(iq)
     577          else
     578              write(*,*) ' Parameter value (default)     : mmol=', &
     579                 mmol(iq)
    520580          end if
    521581          ! option aki
    522582          if (index(tracline,'aki='   ) /= 0) then
    523               read(tracline(index(tracline,'aki=')+len('aki='):),*)
    524      $             aki(iq)
    525               write(*,*) ' Parameter value (traceur.def) : aki=',
    526      $             aki(iq)
    527           else
    528               write(*,*) ' Parameter value (default)     : aki=',
    529      $             aki(iq)
     583              read(tracline(index(tracline,'aki=')+len('aki='):),*) &
     584                  aki(iq)
     585              write(*,*) ' Parameter value (traceur.def) : aki=', &
     586                  aki(iq)
     587          else
     588              write(*,*) ' Parameter value (default)     : aki=', &
     589                  aki(iq)
    530590          end if
    531591          ! option cpi
    532592          if (index(tracline,'cpi='   ) /= 0) then
    533               read(tracline(index(tracline,'cpi=')+len('cpi='):),*)
    534      $             cpi(iq)
    535               write(*,*) ' Parameter value (traceur.def) : cpi=',
    536      $             cpi(iq)
    537           else
    538               write(*,*) ' Parameter value (default)     : cpi=',
    539      $             cpi(iq)
     593              read(tracline(index(tracline,'cpi=')+len('cpi='):),*) &
     594                  cpi(iq)
     595              write(*,*) ' Parameter value (traceur.def) : cpi=', &
     596                  cpi(iq)
     597          else
     598              write(*,*) ' Parameter value (default)     : cpi=', &
     599                  cpi(iq)
    540600          end if
    541601          ! option is_chim
    542602          if (index(tracline,'is_chim='   ) /= 0) then
    543           read(tracline(index(tracline,'is_chim=')+len('is_chim='):),*)
    544      $             is_chim(iq)
    545               write(*,*) ' Parameter value (traceur.def) : is_chim=',
    546      $             is_chim(iq)
    547           else
    548               write(*,*) ' Parameter value (default)     : is_chim=',
    549      $             is_chim(iq)
     603          read(tracline(index(tracline,'is_chim=')+len('is_chim='):),*) &
     604                  is_chim(iq)
     605              write(*,*) ' Parameter value (traceur.def) : is_chim=', &
     606                  is_chim(iq)
     607          else
     608              write(*,*) ' Parameter value (default)     : is_chim=', &
     609                  is_chim(iq)
    550610          end if
    551611          ! option is_rad
    552612          if (index(tracline,'is_rad=') /= 0) then
    553               read(tracline(index(tracline,'is_rad=')
    554      $         +len('is_rad='):),*) is_rad(iq)
    555               write(*,*) ' Parameter value (traceur.def) : is_rad=',
    556      $         is_rad(iq)
    557           else
    558               write(*,*) ' Parameter value (default)     : is_rad=',
    559      $         is_rad(iq)
     613              read(tracline(index(tracline,'is_rad=') &
     614              +len('is_rad='):),*) is_rad(iq)
     615              write(*,*) ' Parameter value (traceur.def) : is_rad=', &
     616              is_rad(iq)
     617          else
     618              write(*,*) ' Parameter value (default)     : is_rad=',  &
     619              is_rad(iq)
    560620          end if
    561621          ! option is_recomb
    562622          if (index(tracline,'is_recomb=') /= 0) then
    563               read(tracline(index(tracline,'is_recomb=')
    564      $         +len('is_recomb='):),*) is_recomb(iq)
    565               write(*,*) ' Parameter value (traceur.def) : is_recomb=',
    566      $         is_recomb(iq)
    567           else
    568               write(*,*) ' Parameter value (default)     : is_recomb=',
    569      $         is_recomb(iq)
     623              read(tracline(index(tracline,'is_recomb=') &
     624              +len('is_recomb='):),*) is_recomb(iq)
     625              write(*,*) ' Parameter value (traceur.def) : is_recomb=', &
     626              is_recomb(iq)
     627          else
     628              write(*,*) ' Parameter value (default)     : is_recomb=', &
     629              is_recomb(iq)
    570630          end if
    571631          ! option is_recomb_qset
    572632          if (index(tracline,'is_recomb_qset=') /= 0) then
    573               read(tracline(index(tracline,'is_recomb_qset=')
    574      $         +len('is_recomb_qset='):),*) is_recomb_qset(iq)
    575               write(*,*) ' Parameter value (traceur.def) :'//
    576      $         ' is_recomb_qset=',
    577      $         is_recomb_qset(iq)
    578           else
    579               write(*,*) ' Parameter value (default)     :'//
    580      $         ' is_recomb_qset=',
    581      $         is_recomb_qset(iq)
     633              read(tracline(index(tracline,'is_recomb_qset=') &
     634              +len('is_recomb_qset='):),*) is_recomb_qset(iq)
     635              write(*,*) ' Parameter value (traceur.def) :'// &
     636              ' is_recomb_qset=', &
     637              is_recomb_qset(iq)
     638          else
     639              write(*,*) ' Parameter value (default)     :'// &
     640              ' is_recomb_qset=', &
     641              is_recomb_qset(iq)
    582642          endif
    583643          ! option is_recomb_qotf
    584644          if (index(tracline,'is_recomb_qotf=') /= 0) then
    585               read(tracline(index(tracline,'is_recomb_qotf=')
    586      $         +len('is_recomb_qotf='):),*) is_recomb_qotf(iq)
    587               write(*,*) ' Parameter value (traceur.def) :'//
    588      $         ' is_recomb_qotf=',
    589      $         is_recomb_qotf(iq)
    590           else
    591               write(*,*) ' Parameter value (default)     :'//
    592      $         ' is_recomb_qotf=',
    593      $         is_recomb_qotf(iq)
     645              read(tracline(index(tracline,'is_recomb_qotf=') &
     646              +len('is_recomb_qotf='):),*) is_recomb_qotf(iq)
     647              write(*,*) ' Parameter value (traceur.def) :'// &
     648              ' is_recomb_qotf=', &
     649              is_recomb_qotf(iq)
     650          else
     651              write(*,*) ' Parameter value (default)     :'// &
     652              ' is_recomb_qotf=',  &
     653              is_recomb_qotf(iq)
    594654          end if
     655          !option is_condensable (LT)
     656          if (index(tracline,'is_condensable=') /=0) then
     657            read(tracline(index(tracline,'is_condensable=') &
     658              +len('is_condensable='):),*) is_condensable(iq)
     659            write(*,*) ' Parameter value (traceur.def) :'// &
     660              ' is_condensable=', is_condensable(iq)       
     661          else
     662              write(*,*) ' Parameter value (default)     :'// &
     663              ' is_condensable=', is_condensable(iq)       
     664          endif
     665          !option radius
     666          if (index(tracline,'radius=') .ne. 0) then
     667            read(tracline(index(tracline,'radius=') &
     668              +len('radius='):),*) radius(iq)
     669            write(*,*)'Parameter value (traceur.def) :'// &
     670              "radius=",radius(iq), "m"
     671          else
     672            write(*,*) ' Parameter value (default)     :'// &
     673            ' radius=', radius(iq)," m"     
     674          endif
     675          !option rho
     676          if (index(tracline,'rho=') .ne. 0) then
     677            read(tracline(index(tracline,'rho=') &
     678              +len('rho='):),*) rho_q(iq)
     679            write(*,*)'Parameter value (traceur.def) :'// &
     680              "rho=",rho_q(iq)
     681          else
     682            write(*,*) ' Parameter value (default)     :'// &
     683              ' rho=', rho_q(iq)       
     684          endif
    595685      end subroutine get_tracdat
    596686
Note: See TracChangeset for help on using the changeset viewer.