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

Timestamp:

Apr 25, 2025, 6:26:55 PM (3 months ago)

Author:

jbclement

Message:

1D Mars PCM:

• The value for atmospheric water profile is now imposed as a uniform mixing ratio and not.
• When atmospheric water profile is monitored, ice value is not put to 0 anymore.
• The user can now prescribe a non-uniform atmospheric water profile with a file "profile_def_h2o_vap".

JBC

Location:

trunk/LMDZ.MARS/libf/phymars/dyn1d

Files:

• init_testphys1d_mod.F90 (modified) (8 diffs)
• testphys1d.F90 (modified) (3 diffs)

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

@@ -7 +7 @@
 SUBROUTINE init_testphys1d(start1Dname,startfiname,therestart1D,therestartfi,ngrid,nlayer,odpref,               &
                            nq,q,time,psurf,u,v,temp,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w,  &
-                           play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
+                           play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau,q_def_h2o_vap)
 
 use ioipsl_getincom,          only: getin ! To use 'getin'
@@ -71 +71 @@
 real,                                intent(out) :: ptif, pks
 real,                                intent(out) :: dttestphys                   ! testphys1d timestep
-real, dimension(:),     allocatable, intent(out) :: zqsat                        ! useful for (atm_wat_profile=2)
+real, dimension(:),     allocatable, intent(out) :: zqsat
 real, dimension(:,:,:), allocatable, intent(out) :: dq, dqdyn                    ! Physical and dynamical tandencies
 integer,                             intent(out) :: day0                         ! initial (sol ; =0 at Ls=0) and final date
@@ -81 +81 @@
 real, dimension(1),                  intent(out) :: latitude, longitude, cell_area
 real,                                intent(out) :: atm_wat_profile, atm_wat_tau ! Force atmospheric water profiles
+real, dimension(nlayer),             intent(out) :: q_def_h2o_vap                ! User-defined atmospheric water profile
 
 !=======================================================================
@@ -94 +95 @@
 real, dimension(:,:),     allocatable :: psdyn
 
-! RV & JBC: Use of starting files for 1D
+! Use of starting files for 1D
 logical                :: found
 character(30)          :: header
@@ -778 +779 @@
 ! Check if the atmospheric water profile is specified
 ! ---------------------------------------------------
-! Adding an option to force atmospheric water values JN
+! Adding an option to force atmospheric water values
 atm_wat_profile = -1. ! Default: free atmospheric water profile
+q_def_h2o_vap = 0. ! Default
 if (water) then
-    write(*,*)'Force atmospheric water vapor profile (uniform column value in kg.m-2)?'
+    write(*,*)'Force atmospheric water vapor profile (mixing ratio in kg/kg)?'
     call getin('atm_wat_profile',atm_wat_profile)
     write(*,*) 'atm_wat_profile = ', atm_wat_profile
@@ -789 +791 @@
     else if (abs(atm_wat_profile) < 1.e-15) then ! if == 0.
         write(*,*) 'Dry atmospheric water vapor profile'
-    else if (0. < atm_wat_profile) then
-        write(*,*) 'Prescribed atmospheric water vapor profile'
+    else if (0. < atm_wat_profile .and. atm_wat_profile <= 1.) then ! if 0. <= and <= 1.
+        write(*,*) 'Prescribed uniform atmospheric water vapor profile'
         write(*,*) 'Unless it reaches saturation (maximal value)'
+    else if (abs(atm_wat_profile + 2.) < 1.e-15) then ! if == -2.
+        write(*,*) 'Prescribed user-defined atmospheric water vapor profile'
+        write(*,*) 'Unless it reaches saturation (maximal value)'
+        open(10,file = 'profile_def_h2o_vap',status = 'old',form = 'formatted',action = 'read',iostat = ierr)
+        if (ierr == 0) then
+            do ilayer = 1,nlayer
+                read(10,*) q_def_h2o_vap(ilayer)
+            enddo
+        else
+            error stop 'File "profile_def_h2o_vap" was not found!'
+        endif 
+        close(10)
     else
-        error stop 'Water vapor profile value not correct!'
+        error stop 'Water vapor profile value is not correct!'
     endif
 endif
@@ -799 +813 @@
 ! Check if the atmospheric water profile relaxation is specified
 ! --------------------------------------------------------------
-! Adding an option to relax atmospheric water values JBC
+! Adding an option to relax atmospheric water values
 atm_wat_tau = -1. ! Default: no time relaxation
 if (water) then
@@ -808 +822 @@
         write(*,*) 'Atmospheric water vapor profile is not relaxed (forcing).'
     else
-        if (0. <= atm_wat_profile) then
+        if (0. <= atm_wat_profile .and. atm_wat_profile <= 1. .or. abs(atm_wat_profile + 2.) < 1.e-15) then
             write(*,*) 'Relaxed atmospheric water vapor profile towards ', atm_wat_profile
             write(*,*) 'Unless it reaches saturation (maximal value)'

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

@@ -81 +81 @@
 real, dimension(1)      :: latitude, longitude, cell_area
 
-! RV & JBC: Use of starting files for 1D
+! Use of starting files for 1D
 logical :: startfiles_1D, therestart1D, therestartfi, there
 
-! JN & JBC: Force atmospheric water profiles
+! Force atmospheric water profiles
 real                            :: atm_wat_profile, atm_wat_tau
-real, dimension(:), allocatable :: zqsat ! useful for (atm_wat_profile=2)
+real, dimension(nlayer)         :: q_def_h2o_vap
+real, dimension(:), allocatable :: zqsat
 !=======================================================================
 
@@ -149 +150 @@
 call init_testphys1d('start1D.txt','startfi.nc',therestart1D,therestartfi,ngrid,nlayer,odpref,nq,q, &
                      time,psurf,u,v,temp,ndt,ptif,pks,dttestphys,zqsat,dq,dqdyn,day0,day,gru,grv,w, &
-                     play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau)
+                     play,plev,latitude,longitude,cell_area,atm_wat_profile,atm_wat_tau,q_def_h2o_vap)
 
 ! Write a "startfi" file
@@ -196 +197 @@
 
     ! Modify atmospheric water if asked
-    ! Added "atm_wat_profile" flag (JN + JBC)
-    ! ---------------------------------------
+    ! Added "atm_wat_profile" flag
+    ! ----------------------------
     if (water) then
         call watersat(nlayer,temp,play,zqsat)
-        if (0. <= atm_wat_profile) then ! If atmospheric water is monitored
+        if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) then ! If atmospheric water is monitored with a uniform value
             if (atm_wat_tau < 0.) then ! Case for prescribed atm_wat_profile: wet if >0, dry if =0
-                q(1,:,igcm_h2o_vap) = min(zqsat,atm_wat_profile*g/psurf)
-                q(1,:,igcm_h2o_ice) = 0. ! reset h2o ice
-            else !  Case for relaxation towards the value atm_wat_profile with relaxation time atm_wat_tau
-                q(1,:,igcm_h2o_vap) = atm_wat_profile*g/psurf + (q(1,:,igcm_h2o_vap) - atm_wat_profile*g/psurf)*dexp(-dttestphys/atm_wat_tau)
-                q(1,:,igcm_h2o_ice) = 0. ! reset h2o ice
+                q(1,:,igcm_h2o_vap) = min(zqsat,atm_wat_profile)
+            else ! Case for relaxation towards the value atm_wat_profile with relaxation time atm_wat_tau
+                q(1,:,igcm_h2o_vap) = atm_wat_profile + (q(1,:,igcm_h2o_vap) - atm_wat_profile)*dexp(-dttestphys/atm_wat_tau)
+            endif
+        else if (abs(atm_wat_profile + 2.) < 1.e-15) then ! If atmospheric water is monitored with a user-defined water profile (with file "profile_def_h2o_vap")
+            if (atm_wat_tau < 0.) then ! Case for prescribed atm_wat_profile: wet if >0, dry if =0
+                q(1,:,igcm_h2o_vap) = min(zqsat(:),q_def_h2o_vap(:))
+            else ! Case for relaxation towards the value atm_wat_profile with relaxation time atm_wat_tau
+                q(1,:,igcm_h2o_vap) = q_def_h2o_vap(:) + (q(1,:,igcm_h2o_vap) - q_def_h2o_vap(:))*dexp(-dttestphys/atm_wat_tau)
             endif
         endif