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

Timestamp:

May 11, 2023, 5:40:03 PM (20 months ago)

Author:

romain.vande

Message:

Mars PEM :

Adapt PEM to the subslope PCM configuration, it is now fully compatible.

Create a PEM folder in deftank that contains:

run_pem1: a bash file that runs chained simulation of PEM as well as running a parameterizable number of PCM simulation in between.

It also takes care of reshaping XIOS output as well as renaming outputs etc… in the spirit of run_month1.

It is written for Irene machine and the header needs to be adapted for other machines.

run_PEM.def: A text file that shows the possible parameters to choose before a PEM simulation.

It should be included at the end of run.def just like callphys.def

ob_ex_lsp.asc: An ascii file containing the obliquity, eccentricity, ls_peri data from Laskar in Martian year.
README: A txt file explaining the content of the folder

Adapt field_def_physics_mars.xml to consider the case with 7 subslopes in the PCM.
Change context_lmdz_physics.xml to be able to output the file needed by the PEM.

Correct a few other minor bugs.

RV

Location:

trunk/LMDZ.MARS/libf/phymars

Files:

• call_dayperi.F (modified) (1 diff)
• iniorbit.F (modified) (1 diff)
• physiq_mod.F (modified) (6 diffs)
• rocketduststorm_mod.F90 (modified) (1 diff)

trunk/LMDZ.MARS/libf/phymars/call_dayperi.F

@@ -17 +17 @@
 c   Lsperi      solar longitude (Ls) of perohelion (rad)
 c   e_elips       Excentricity
+c   real year_day      ! number of sols per Mars yar
 c
 c   output
 c   ------
 c   dayperi       Martian date at perihelion (sol)
-c   real year_day      ! number of sols per Mars yar
+
 c-----------------------------------------------------------------------
 

trunk/LMDZ.MARS/libf/phymars/iniorbit.F

@@ -68 +68 @@
 
       timeperi=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.))
+      if(timeperi.lt.0) timeperi=-timeperi
       PRINT*,'iniorbit: Perihelion solar long. Ls (deg)=',
      &       360.-timeperi*180./pi

trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

@@ -545 +545 @@
       logical :: write_restart
 
+! Variable for ice table
+      REAL :: rhowater_surf(ngrid,nslope)
+      REAL :: rhowater_surf_sat(ngrid,nslope)
+      REAL :: rhowater_soil(ngrid,nsoilmx,nslope)
+      REAL,PARAMETER  :: alph_clap = -6143.7
+      REAL,PARAMETER :: beta_clap = 28.9074
+      REAL :: pvap_surf(ngrid)
+      REAL,PARAMETER :: m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)   
+      REAL,PARAMETER :: m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)   
+      REAL :: ztmp1,ztmp2
+
 c=======================================================================
       pdq(:,:,:) = 0.
@@ -789 +800 @@
       call update_xios_timestep     
 #endif 
+
+
 
 
@@ -1498 +1511 @@
              dwatercap(ig,:)=dwatercap(ig,:)+dwatercap_dif(ig,:)
           ENDDO
+
           call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,zdqsdif,
      &   albedo_meshavg,emis_meshavg,tsurf_meshavg,zdqsdif_meshavg_tmp)
@@ -1559 +1573 @@
            ENDDO
          ENDDO
+
          IF (turb_resolved) THEN
             write(*,*) 'Turbulent-resolving mode !' 
@@ -2160 +2175 @@
       endif ! of if (callthermos)
 #endif
+
 c-----------------------------------------------------------------------
 c   11. Carbon dioxide condensation-sublimation:
@@ -3822 +3838 @@
            endif !not.scavenging
         ENDIF ! of IF (water)
+
+! Output needed by the PEM
+          DO ig = 1,ngrid
+            ztmp1 =(1/m_co2 - 1/m_noco2)
+            ztmp2=1/m_noco2
+            pvap_surf(ig) = 1/(ztmp1*zq(ig,1,igcm_co2)+ztmp2)
+     &      * zq(ig,1,igcm_h2o_vap)/(mmol(igcm_h2o_vap)*1.e-3)*ps(ig)
+
+            DO islope = 1,nslope
+             rhowater_surf_sat(ig,islope)  =
+     &         exp(alph_clap/tsurf(ig,islope)+beta_clap)
+     &         / tsurf(ig,islope)
+ 
+             if(qsurf(ig,igcm_h2o_ice,islope).gt.(1.e-4)) then
+               rhowater_surf(ig,islope) =
+     &         exp(alph_clap/tsurf(ig,islope)+beta_clap)
+     &         / tsurf(ig,islope)
+             else
+              rhowater_surf(ig,islope) = pvap_surf(ig)
+     &         / tsurf(ig,islope)
+             endif
+             DO isoil = 1,nsoilmx
+             rhowater_soil(ig,isoil,islope) =
+     &         exp(alph_clap/tsoil(ig,isoil,islope)+beta_clap)
+     &         / tsoil(ig,isoil,islope)
+             ENDDO
+          ENDDO
+        ENDDO
+
+      DO islope = 1,nslope
+        write(str2(1:2),'(i2.2)') islope
+        CALL send_xios_field("Waterdensity_soil_slope"//str2,
+     &     rhowater_soil(:,:,islope))
+        CALL send_xios_field("Waterdensity_surface"//str2,
+     &     rhowater_surf(:,islope))
+      ENDDO
+
+      CALL send_xios_field("h2o_layer1",zq(:,1,igcm_h2o_vap))
+      CALL send_xios_field("co2_layer1",zq(:,1,igcm_co2))
+
 !PREVIOUSLY IN 1D ONLY
 

trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90

@@ -532 +532 @@
                           'lapse rate in the rocket dust storm', &
                                     'K/m',lapserate(:,:))
-        call write_output('rds_deltahr', &
-                          'extra heating rate in the rocket dust storm', &
-                                    'K/s',deltahr(:,:))
+!        call write_output('rds_deltahr', &
+!                          'extra heating rate in the rocket dust storm', &
+!                                    'K/s',deltahr(:,:))
 !        call write_output('scheme','which scheme',&
 !                                                   ' ',real(scheme(:)))