Changeset 3054

Timestamp:

Sep 27, 2023, 11:14:02 AM (14 months ago)

Author:

jbclement

Message:

Mars PCM 1D:
Correction of a bug: 'inertiedat(1,1)' was overwritten by 'inertieice'.
JBC

Location:

trunk/LMDZ.MARS

Files:

• changelog.txt (modified) (1 diff)
• libf/phymars/dyn1d/testphys1d.F90 (modified) (14 diffs)

trunk/LMDZ.MARS/changelog.txt

@@ -4214 +4214 @@
 == 27/09/2023 == JBC
 Upgrade of "testphys1d" to Fortran 90. Cleaning of the subroutine and minor optimizations of the code.
+Correction of a bug: 'inertiedat(1,1)' was overwritten by 'inertieice'.

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

@@ -35 +35 @@
 use iostart,                  only: open_startphy, get_var, close_startphy
 use write_output_mod,         only: write_output
-! Mostly for XIOS outputs: 
+! Mostly for XIOS outputs:
 use mod_const_mpi,            only: init_const_mpi, COMM_LMDZ
 use parallel_lmdz,            only: init_parallel
@@ -190 +190 @@
 
 if (startfiles_1D) then
-    inquire(file = 'start1D.txt',exist = therestart1D) 
+    inquire(file = 'start1D.txt',exist = therestart1D)
     if (.not. therestart1D) then
         write(*,*) 'There is no "start1D.txt" file!'
         write(*,*) 'Initialization is done with default values.'
     endif
-    inquire(file = 'startfi.nc',exist = therestartfi) 
+    inquire(file = 'startfi.nc',exist = therestartfi)
     if (.not. therestartfi) then
         write(*,*) 'There is no "startfi.nc" file!'
@@ -224 +224 @@
 obliquit = 25.2                  ! Obliquity (deg) ~25.2
 excentric = 0.0934               ! Eccentricity (0.0934)
- 
+
 ! Planetary Boundary Layer and Turbulence parameters
 ! --------------------------------------------------
@@ -230 +230 @@
 emin_turb = 1.e-6  ! minimal turbulent energy ~1.e-8
 lmixmin = 30       ! mixing length ~100
- 
+
 ! cap properties and surface emissivities
 ! ---------------------------------------
@@ -312 +312 @@
                 continu = .false.
             else !if (tnom_transp(iq) == tnom_0(ipere)) then
-                ipere = ipere + 1 
+                ipere = ipere + 1
                 if (ipere > nqtot) then
                     write(*,*) 'Le traceur',iq,'appele ',trim(tname(iq)),', est orpelin.'
@@ -406 +406 @@
 pa = 20.     ! transition pressure (for hybrid coord.)
 preff = 610. ! reference surface pressure
- 
+
 ! Aerosol properties
 ! ------------------
@@ -425 +425 @@
         write(*,*) 'Martian eccentricity (0<e<1)?'
         call getin('excentric ',excentric)
-        write(*,*)"excentric =",excentric 
+        write(*,*)"excentric =",excentric
         write(*,*) 'Solar longitude of perihelion (0<Ls<360)?'
         call getin('Lsperi',Lsperi )
-        write(*,*)"Lsperi=",Lsperi 
+        write(*,*)"Lsperi=",Lsperi
         Lsperi = Lsperi*pi/180.0 ! Put it in rad for peri_day
         periheli = halfaxe*(1 - excentric)
@@ -458 +458 @@
      endif
 endif !(.not. startfiles_1D )
- 
+
 ! Latitude/longitude
 ! ------------------
@@ -528 +528 @@
 def_slope(1) = -90 ! minimum slope angle
 def_slope(2) = 90  ! maximum slope angle
-subslope_dist(1,1) = 1 ! fraction of subslopes in mesh 
+subslope_dist(1,1) = 1 ! fraction of subslopes in mesh
 
 ! For the gravity wave scheme
@@ -611 +611 @@
         if (latitude(1) < 0) emis = emisice(2) ! southern hemisphere
     endif
-endif !(.not. startfiles_1D ) 
+endif !(.not. startfiles_1D )
 
 ! Compute pressures and altitudes of atmospheric levels
@@ -661 +661 @@
 
     ! Creating the new soil inertia table if there is subsurface ice:
-    if (ice_depth > 0) then 
+    if (ice_depth > 0) then
         iref = 1 ! ice/regolith boundary index
         if (ice_depth < layer(1)) then
             inertiedat(1,1) = sqrt(layer(1)/((ice_depth/inertiedat(1,1)**2) + ((layer(1) - ice_depth)/inertieice**2)))
-            inertiedat(1,:) = inertieice
+            inertiedat(1,2:) = inertieice
         else ! searching for the ice/regolith boundary:
-            do isoil = 1,nsoil 
+            do isoil = 1,nsoil
                 if ((ice_depth >= layer(isoil)) .and. (ice_depth < layer(isoil + 1))) then
                     iref = isoil + 1
@@ -680 +680 @@
             inertiedat(1,iref + 1:) = inertieice
         endif ! (ice_depth < layer(1))
-    else ! ice_depth <0 all is set to surface thermal inertia
+    else ! ice_depth < 0 all is set to surface thermal inertia
         inertiedat(1,:) = inertiedat(1,1) ! soil thermal inertia
     endif ! ice_depth > 0
@@ -888 +888 @@
 
 END PROGRAM testphys1d
- 
+
 !***********************************************************************
 !***********************************************************************
@@ -909 +909 @@
 
 END SUBROUTINE gr_fi_dyn
- 
+
 !***********************************************************************
 !***********************************************************************
-