Changeset 3848 for trunk/LMDZ.MARS/libf/phymars

Timestamp:

Jul 15, 2025, 12:13:09 PM (4 weeks ago)

Author:

jbclement

Message:

Mars PCM 1D:
In the 1D model, if no input profile found, the 'dust_mass' and 'dust_number' tracers are initialized according to "Conrath dust" in line with what it's done for CO2 and HDO, instead of being set to 0 which caused bugs.
JBC

Location:

trunk/LMDZ.MARS/libf/phymars/dyn1d

Files:

• init_testphys1d_mod.F90 (modified) (4 diffs)
• read_profile_mod.F90 (modified) (9 diffs)

trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90

@@ -319 +319 @@
 write(*,*) " psurf = ",psurf
 
+! Compute pressures and altitudes of atmospheric levels
+! -----------------------------------------------------
+! Vertical Coordinates
+! """"""""""""""""""""
+hybrid = .true.
+write(*,*)'Hybrid coordinates?'
+call getin("hybrid",hybrid)
+write(*,*) " hybrid = ", hybrid
+
+! Reference pressures
+pa = 20.     ! Transition pressure (for hybrid coord.)
+preff = 610. ! Reference surface pressure
+
+call disvert_noterre
+! Now that disvert has been called, initialize module vertical_layers_mod
+call init_vertical_layers(nlayer,preff,scaleheight,ap,bp,aps,bps,presnivs,pseudoalt)
+
+plev = ap + psurf*bp
+play = aps + psurf*bps
+zlay = -200.*r*log(play/plev(1))/g
+
 ! Coefficient to control the surface pressure change
 ! To be defined in "callphys.def" so that the PEM can read it asa well
@@ -331 +352 @@
     error stop 'Please, specify a correct value!'
 endif
-
-! Reference pressures
-pa = 20.     ! Transition pressure (for hybrid coord.)
-preff = 610. ! Reference surface pressure
 
 ! Aerosol properties
@@ -512 +529 @@
 write(*,*) "init_testphys1d: initializing tracers"
 if (.not. therestart1D) then
-    call read_profile(nq,nlayer,qsurf(1,:,1),q)
+    call read_profile(nq,nlayer,qsurf(1,:,1),q,play)
     do iq = 1,nq
         qsurf(1,iq,:) = qsurf(1,iq,1)
@@ -631 +648 @@
     endif
 endif !(.not. therestartfi)
-
-! Compute pressures and altitudes of atmospheric levels
-! -----------------------------------------------------
-! Vertical Coordinates
-! """"""""""""""""""""
-hybrid = .true.
-write(*,*)'Hybrid coordinates?'
-call getin("hybrid",hybrid)
-write(*,*) " hybrid = ", hybrid
-
-call disvert_noterre
-! Now that disvert has been called, initialize module vertical_layers_mod
-call init_vertical_layers(nlayer,preff,scaleheight,ap,bp,aps,bps,presnivs,pseudoalt)
-
-plev = ap + psurf*bp
-play = aps + psurf*bps
-zlay = -200.*r*log(play/plev(1))/g
 
 ! Initialize temperature profile

trunk/LMDZ.MARS/libf/phymars/dyn1d/read_profile_mod.F90

@@ -19 +19 @@
 !          - q(co2) = 0.95
 !          - q(hdo) = q(h2o)*2*155.76e-6*5
+!          - q(dust_mass) and q(dust_number) initialized with Conrath dust
 !
 !     - igcm is not available at this part of the code
@@ -36 +37 @@
 !   4. Traitment for minor isotopologue if the major isotopologue input profile does not exist
 !=====================================================================================================================!
-  subroutine read_profile(nb_tracer, nb_layer, qsurf, q)
-
-  use infotrac, only: tname
-  use tracer_mod, only: igcm_co2, igcm_h2o_vap, igcm_h2o_ice, &
-                        igcm_dust_number, igcm_dust_mass,     &
-                        igcm_ccn_number, igcm_ccn_mass
+  subroutine read_profile(nb_tracer, nb_layer, qsurf, q, play)
+
+  use infotrac,       only: tname
+  use tracer_mod,     only: igcm_co2, igcm_h2o_vap, igcm_h2o_ice, &
+                            igcm_dust_number, igcm_dust_mass,     &
+                            igcm_ccn_number, igcm_ccn_mass
+  use aeropacity_mod, only: topdustref
+  use dust_param_mod, only: odpref
+  use comcstfi_h,     only: pi
+
   implicit none
 !---------------------------------------------------------------------------------------------------------------------!
@@ -50 +55 @@
   integer, intent(in) :: &
      nb_layer, & ! number of layer
-     nb_tracer  ! number of traceur read from traceur.def 
+     nb_tracer   ! number of traceur read from traceur.def 
+  real, dimension(nb_layer), intent(in) :: play ! Pressure at the middle of the layers (Pa)
 !---------------------------------------------------------------------------------------------------------------------!
 !  Output arguments:
@@ -61 +67 @@
 !--------
   integer :: &
-     iq,                 & ! loop over nb_tracer
-     ilayer,             & ! loop over nb_layer
-     ierr,               & ! open file iostat
-     indice_h2o_vap = 0, & ! indice of h2o_vap tracer 
-     indice_h2o_ice = 0, & ! indice of h2o_ice tracer
-     indice_hdo_vap = 0, & ! indice of hdo_vap tracer
-     indice_hdo_ice = 0    ! indice of hdo_ice tracer
+     iq,                     & ! loop over nb_tracer
+     ilayer,                 & ! loop over nb_layer
+     ierr,                   & ! open file iostat
+     indice_h2o_vap     = 0, & ! indice of h2o_vap tracer 
+     indice_h2o_ice     = 0, & ! indice of h2o_ice tracer
+     indice_hdo_vap     = 0, & ! indice of hdo_vap tracer
+     indice_hdo_ice     = 0, & ! indice of hdo_ice tracer
+     indice_dust_mass   = 0, & ! indice of hdo_vap tracer
+     indice_dust_number = 0    ! indice of hdo_ice tracer
 
   character(len=80), dimension(nb_tracer) :: &
@@ -73 +81 @@
 
   logical :: &
-     hdo_vap = .false., & ! used to compute hdo_vap profile if its input profile is missing (= .true)
-     hdo_ice = .false.    ! used to compute hdo_ice profile if its input profile is missing (= .true)
+     hdo_vap = .false., & ! used to compute hdo_vap profile if its input profile is missing (= .true.)
+     hdo_ice = .false., & ! used to compute hdo_ice profile if its input profile is missing (= .true.)
+     dust = .false.       ! used to compute dust mass & number profiles if their input profiles are missing (= .true.)
+
+  real :: zp, varian, r3n_q, ref_r0, r0_lift, mass2number_lift ! Dust variables
+  real, parameter ::    &
+     nueff_lift = 0.5,  &
+     reff_lift = 3.e-6, & ! Effective radius of lifted dust (m)
+     rho_dust = 2500      ! Mars dust density (kg.m-3)
 !=====================================================================================================================!
 !=== BEGIN                                                                                                            !
@@ -96 +111 @@
     else if (trim(name_tracer(iq)) == 'hdo_ice') then
       indice_hdo_ice = iq
+    else if (trim(name_tracer(iq)) == 'dust_mass') then
+      indice_dust_mass = iq
+    else if (trim(name_tracer(iq)) == 'dust_number') then
+      indice_dust_number = iq
     end if
   end do
@@ -126 +145 @@
 !---------------------------------------------------------------------------------------------------------------------!
     else
-      write(*,*)"File profile_"//trim(name_tracer(iq))//" not found!"
-      write(*,*)"q(traceur) and qsurf(traceur) are set to 0! Except:"
-      write(*,*)"    - q(co2) = 0.95"
-      write(*,*)"    - q(hdo) = q(h2o)*2*155.76e-6*5"
+      write(*,*)"Warning: file profile_"//trim(name_tracer(iq))//" not found!"
 !---------------------------------------------------------------------------------------------------------------------!
 ! 3.1.2.a. Some cases require a special initialization
@@ -135 +151 @@
       select case(trim(name_tracer(iq)))
         case("co2")
-          q(1:nb_layer, iq) = 0.95
+          q(1:nb_layer,iq) = 0.95
+          write(*,*)"Initilization: q(co2) = 0.95"
 
         case("hdo_vap")
           hdo_vap = .true.
+          write(*,*)"Initialization: q(hdo) = q(h2o)*2*155.76e-6*5"
 
         case("hdo_ice")
           hdo_ice = .true.
+          write(*,*)"Initialization: q(hdo) = q(h2o)*2*155.76e-6*5"
+
+        case("dust_mass","dust_number")
+          dust = .true.
+          write(*,*)"q(dust_mass) and q(dust_number) initialized with Conrath dust"
+
+        case default
+          write(*,*)"q("//trim(name_tracer(iq))//") and qsurf("//trim(name_tracer(iq))//") are set to 0!"
 
        end select
@@ -163 +189 @@
     end do
   end if
+
+  if (dust) then
+    varian = sqrt(log(1. + nueff_lift))
+    r3n_q = exp(-4.5*varian**2)*(3./4.)/(pi*rho_dust)
+    ref_r0 = exp(2.5*varian**2)
+    r0_lift = reff_lift/ref_r0
+    mass2number_lift = r3n_q/r0_lift**3
+    do ilayer = 1,nb_layer
+      zp = (odpref/play(ilayer))**(70./topdustref)
+      q(ilayer,indice_dust_mass) = max(exp(0.007*(1. - max(zp,1.))),1.e-3)
+      q(ilayer,indice_dust_number) = q(ilayer,indice_dust_mass)*mass2number_lift
+    enddo
+  endif
 !=====================================================================================================================!
 !=== END                                                                                                              !