Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 2354)
+++ /trunk/LMDZ.MARS/README	(revision 2355)
@@ -3087,2 +3087,6 @@
 - a .def file to use as input for the main program
 
+== 09/06/2020 == CM
+Improvement in the reading of input profiles for 1D simulation.
+ - add: dyn1d/read_profile_mod.F90
+
Index: /trunk/LMDZ.MARS/libf/phymars/dyn1d/read_profile_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/dyn1d/read_profile_mod.F90	(revision 2355)
+++ /trunk/LMDZ.MARS/libf/phymars/dyn1d/read_profile_mod.F90	(revision 2355)
@@ -0,0 +1,158 @@
+module read_profile_mod
+
+implicit none
+
+contains
+!=====================================================================================================================!
+! SUBROUTINE: read_profile ===========================================================================================!
+!=====================================================================================================================!
+! Author: Christophe Mathe
+! Date: 09/06/2020
+!---------------------------------------------------------------------------------------------------------------------!
+! Subject:
+!---------
+!   Read input profile of tracers listed in traceur.def
+!---------------------------------------------------------------------------------------------------------------------!
+! Comments:
+!----------
+!     - If no input profile found, q(traceur) and qsurf(traceur) set to 0, except:
+!          - q(co2) = 0.95
+!          - q(hdo) = q(h2o)*2*155.76e-6*5
+!
+!     - igcm is not available at this part of the code
+!
+!     - To ensure that major isotopologue input profile is read before compute minor isotopologue profile, the
+!       calculation is performed at the end of the subroutine
+!---------------------------------------------------------------------------------------------------------------------!
+! Algorithm:
+!-----------
+!   1. Initialization of q and qsurf to 0
+!   2. Get indices of some tracers
+!   3. Main
+!     3.1. Try to open the input profile
+!       3.1.1. Succeed
+!       3.1.2. Fail
+!         3.1.2.a. Some cases require a special initialization
+!   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
+
+  implicit none
+!---------------------------------------------------------------------------------------------------------------------!
+! VARIABLE DECLARATION
+!---------------------------------------------------------------------------------------------------------------------!
+!  Input arguments:
+!------------------
+  integer, intent(in) :: &
+     nb_layer, & ! number of layer
+     nb_tracer  ! number of traceur read from traceur.def 
+!---------------------------------------------------------------------------------------------------------------------!
+!  Output arguments:
+!-------------------
+  real, intent(out) :: &
+     qsurf(nb_tracer), &    ! kg/m2
+     q(nb_layer, nb_tracer) ! kg/kg of atmosphere
+!---------------------------------------------------------------------------------------------------------------------!
+!  Local:
+!--------
+  integer :: &
+     iq,             & ! loop over nb_tracer
+     ilayer,         & ! loop over nb_layer
+     ierr,           & ! open file iostat
+     indice_h2o_vap, & ! indice of h2o_vap tracer 
+     indice_h2o_ice, & ! indice of h2o_ice tracer
+     indice_hdo_vap, & ! indice of hdo_vap tracer
+     indice_hdo_ice    ! indice of hdo_ice tracer
+
+  character(len=80), dimension(nb_tracer) :: &
+     name_tracer ! array of all tracers already read in traceur.def
+
+  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)
+!=====================================================================================================================!
+!=== BEGIN                                                                                                            !
+!=====================================================================================================================!
+! 1. Initialization of q and qsurf to 0
+!---------------------------------------------------------------------------------------------------------------------!
+  q(1:nb_layer,1:nb_tracer) = 0. 
+  qsurf(1:nb_tracer) = 0. 
+  name_tracer(1:nb_tracer) = ""
+!---------------------------------------------------------------------------------------------------------------------!
+! 2. Get indices of some tracers: igcm is not available at this part of the code
+!---------------------------------------------------------------------------------------------------------------------!
+  do iq = 1, nb_tracer
+    write(name_tracer(iq),"(a)") tname(iq)
+    if (trim(name_tracer(iq)) == 'h2o_vap') then
+      indice_h2o_vap = iq
+    else if (trim(name_tracer(iq)) == 'h2o_ice') then
+      indice_h2o_ice = iq
+    else if (trim(name_tracer(iq)) == 'hdo_vap') then
+      indice_hdo_vap = iq
+    else if (trim(name_tracer(iq)) == 'hdo_ice') then
+      indice_hdo_ice = iq
+    end if
+  end do
+!---------------------------------------------------------------------------------------------------------------------!
+! 3. Main
+!---------------------------------------------------------------------------------------------------------------------!
+  do iq = 1, nb_tracer
+    write(*,*)"  tracer:",trim(name_tracer(iq))
+!---------------------------------------------------------------------------------------------------------------------!
+! 3.1. Try to open the input profile
+!---------------------------------------------------------------------------------------------------------------------!
+    open(91,file='profile_'//trim(name_tracer(iq)), status='old', form='formatted', iostat=ierr)
+!---------------------------------------------------------------------------------------------------------------------!
+! 3.1.1. Succeed: the input file exists
+!---------------------------------------------------------------------------------------------------------------------!
+    if (ierr.eq.0) then
+      read(91,*)qsurf(iq)
+      do ilayer = 1, nb_layer
+        read(91,*)q(ilayer,iq)
+      end do
+!---------------------------------------------------------------------------------------------------------------------!
+! 3.1.2. Fail: the input file does not exist
+!---------------------------------------------------------------------------------------------------------------------!
+    else
+      write(*,*)"File profile_"//trim(name_tracer(iq))//" not found!"
+!---------------------------------------------------------------------------------------------------------------------!
+! 3.1.2.a. Some cases require a special initialization
+!---------------------------------------------------------------------------------------------------------------------!
+      select case(trim(name_tracer(iq)))
+        case("co2")
+          q(1:nb_layer, iq) = 0.95
+
+        case("hdo_vap")
+          hdo_vap = .true.
+
+        case("hdo_ice")
+          hdo_ice = .true.
+
+       end select
+    end if
+    close(91)
+  end do ! of do iq = 1, nq
+!---------------------------------------------------------------------------------------------------------------------!
+! 4. Traitment for minor isotopologue if the major isotopologue input profile does not exist. At this part, ensure that
+!    the main isotopologue profile is already read before minors isotopologues
+!---------------------------------------------------------------------------------------------------------------------!
+  if (hdo_vap.eqv..true. .and. indice_h2o_vap.ne.0) then
+    do ilayer = 1, nb_layer
+      q(ilayer, indice_hdo_vap) = q(ilayer, indice_h2o_vap)*2*155.76e-6*5
+    end do
+  end if
+
+  if (hdo_ice.eqv..true. .and. indice_h2o_ice.ne.0) then
+    qsurf(indice_hdo_ice) = qsurf(indice_h2o_ice) * 2*155.76e-6*5
+    do ilayer = 1, nb_layer
+      q(ilayer, indice_hdo_ice) = q(ilayer, indice_h2o_ice) * 2*155.76e-6*5
+    end do
+  end if
+!=====================================================================================================================!
+!=== END                                                                                                              !
+!=====================================================================================================================!
+  end subroutine read_profile
+end module read_profile_mod
+
Index: /trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 2354)
+++ /trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 2355)
@@ -32,4 +32,5 @@
       USE phys_state_var_init_mod, ONLY: phys_state_var_init
       USE physiq_mod, ONLY: physiq
+      USE read_profile_mod, ONLY: read_profile
       IMPLICIT NONE
 
@@ -117,6 +118,4 @@
       integer :: nq=1 ! number of tracers
       real :: latitude(1), longitude(1), cell_area(1)
-      integer iqh2ovap
-      integer iqh2oice
 
       character*2 str2
@@ -307,186 +306,5 @@
         ! initialize tracers here:
         write(*,*) "testphys1d: initializing tracers"
-        q(:,:)=0 ! default, set everything to zero
-        qsurf(:)=0
-        ! "smarter" initialization of some tracers
-        ! (get values from "profile_*" files, if these are available)
-        do iq=1,nq
-          txt=""
-          write(txt,"(a)") tname(iq)
-          write(*,*)"  tracer:",trim(txt)
-          ! CO2
-          if (txt.eq."co2") then
-            q(:,iq)=0.95   ! kg /kg of atmosphere
-            qsurf(iq)=0. ! kg/m2 (not used for CO2)
-            ! even better, look for a "profile_co2" input file
-            open(91,file='profile_co2',status='old',
-     &       form='formatted',iostat=ierr)
-            if (ierr.eq.0) then
-              read(91,*) qsurf(iq)
-              do ilayer=1,nlayer
-                read(91,*) q(ilayer,iq)
-              enddo
-            endif
-            close(91)
-          endif ! of if (txt.eq."co2")
-          ! Allow for an initial profile of argon
-          ! Can also be used to introduce a decaying tracer
-          ! in the 1D (TBD) to study thermals
-          if (txt.eq."ar") then
-            !look for a "profile_ar" input file
-            open(91,file='profile_ar',status='old',
-     &       form='formatted',iostat=ierr)
-            if (ierr.eq.0) then
-              read(91,*) qsurf(iq)
-              do ilayer=1,nlayer
-                read(91,*) q(ilayer,iq)
-              enddo
-            else
-              write(*,*) "No profile_ar file!"
-            endif
-            close(91)
-          endif ! of if (txt.eq."ar")
-
-          ! WATER VAPOUR
-          if (txt.eq."h2o_vap") then
-             iqh2ovap=iq !remember index for water vap
-            !look for a "profile_h2o_vap" input file
-            open(91,file='profile_h2o_vap',status='old',
-     &       form='formatted',iostat=ierr)
-            if (ierr.eq.0) then
-              read(91,*) qsurf(iq)
-              do ilayer=1,nlayer
-                read(91,*) q(ilayer,iq)
-              enddo
-            else
-              write(*,*) "No profile_h2o_vap file!"
-            endif
-            close(91)
-          endif ! of if (txt.eq."h2o_ice")
-          ! WATER ICE
-          if (txt.eq."h2o_ice") then
-                iqh2oice = iq !remember index for water ice
-            !look for a "profile_h2o_vap" input file
-            open(91,file='profile_h2o_ice',status='old',
-     &       form='formatted',iostat=ierr)
-            if (ierr.eq.0) then
-              read(91,*) qsurf(iq)
-              do ilayer=1,nlayer
-                read(91,*) q(ilayer,iq)
-              enddo
-            else
-              write(*,*) "No profile_h2o_ice file!"
-            endif
-            close(91)
-          endif ! of if (txt.eq."h2o_ice")
-
-          ! HDO VAPOUR
-          if (txt.eq."hdo_vap") then
-            !look for a "profile_hdo_vap" input file
-            open(91,file='profile_hdo_vap',status='old',
-     &       form='formatted',iostat=ierr)
-            if (ierr.eq.0) then
-              read(91,*) qsurf(iq)
-              do ilayer=1,nlayer
-                read(91,*) q(ilayer,iq)
-              enddo
-            else
-              write(*,*) "No profile_hdo_vap file!"
-              do ilayer=1,nlayer
-                q(ilayer,iq) =  q(ilayer,iqh2ovap)*2*155.76e-6*5 
-              enddo
-            endif
-            close(91)
-          endif ! of if (txt.eq."hdo_ice")
-          ! HDO ICE
-          if (txt.eq."hdo_ice") then
-            !look for a "profile_hdo_vap" input file
-            open(91,file='profile_hdo_ice',status='old',
-     &       form='formatted',iostat=ierr)
-            if (ierr.eq.0) then
-              read(91,*) qsurf(iq)
-              do ilayer=1,nlayer
-                read(91,*) q(ilayer,iq)
-              enddo
-            else
-              write(*,*) "No profile_hdo_ice file!"
-                qsurf(iq) = qsurf(iqh2oice) * 2*155.76e-6*5
-               do ilayer=1,nlayer
-                q(ilayer,iq) = q(ilayer,iqh2oice) * 2*155.76e-6*5 
-              enddo
-            endif
-            close(91)
-          endif ! of if (txt.eq."hdo_ice")
-
-
-          ! DUST
-          !if (txt(1:4).eq."dust") then
-          !  q(:,iq)=0.4    ! kg/kg of atmosphere
-          !  qsurf(iq)=100 ! kg/m2
-          !endif
-          ! DUST MMR
-          if (txt.eq."dust_mass") then
-            !look for a "profile_dust_mass" input file
-            open(91,file='profile_dust_mass',status='old',
-     &       form='formatted',iostat=ierr)
-            if (ierr.eq.0) then
-              read(91,*) qsurf(iq)
-              do ilayer=1,nlayer
-                read(91,*) q(ilayer,iq)
-!                write(*,*) "l=",ilayer," q(ilayer,iq)=",q(ilayer,iq)
-              enddo
-            else
-              write(*,*) "No profile_dust_mass file!"
-            endif
-            close(91)
-          endif ! of if (txt.eq."dust_mass")
-          ! DUST NUMBER
-          if (txt.eq."dust_number") then
-            !look for a "profile_dust_number" input file
-            open(91,file='profile_dust_number',status='old',
-     &       form='formatted',iostat=ierr)
-            if (ierr.eq.0) then
-              read(91,*) qsurf(iq)
-              do ilayer=1,nlayer
-                read(91,*) q(ilayer,iq)
-              enddo
-            else
-              write(*,*) "No profile_dust_number file!"
-            endif
-            close(91)
-          endif ! of if (txt.eq."dust_number")
-          ! NB: some more initializations (chemistry) is done later
-          ! CCN MASS
-          if (txt.eq."ccn_mass") then
-            !look for a "profile_ccn_mass" input file
-            open(91,file='profile_ccn_mass',status='old',
-     &       form='formatted',iostat=ierr)
-            if (ierr.eq.0) then
-              read(91,*) qsurf(iq)
-              do ilayer=1,nlayer
-                read(91,*) q(ilayer,iq)
-              enddo
-            else
-              write(*,*) "No profile_ccn_mass file!"
-            endif
-            close(91)
-          endif ! of if (txt.eq."ccn_mass")
-          ! CCN NUMBER
-          if (txt.eq."ccn_number") then
-            !look for a "profile_ccn_number" input file
-            open(91,file='profile_ccn_number',status='old',
-     &       form='formatted',iostat=ierr)
-            if (ierr.eq.0) then
-              read(91,*) qsurf(iq)
-              do ilayer=1,nlayer
-                read(91,*) q(ilayer,iq)
-              enddo
-            else
-              write(*,*) "No profile_ccn_number file!"
-            endif
-            close(91)
-          endif ! of if (txt.eq."ccn_number")
-        enddo ! of do iq=1,nq
-
+        call read_profile(nq, nlayer, qsurf, q)
       else
       ! we still need to set (dummy) tracer number and names for physdem1
