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

Timestamp:

Apr 28, 2025, 3:21:28 PM (3 months ago)

Author:

jbclement

Message:

Mars PCM:
Simplification and clarification for the prescription of the atmospheric water profile in 1D: the boolean 'prescribed_h2ovap' activates the option, the file "profile_prescribed_h2ovap" defines the prescribed profile and the value 'h2ovap_relax_time', if positive, activates the relaxation and gives the constant of time.
JBC

Location:

trunk/LMDZ.MARS/libf/phymars/dyn1d

Files:

• init_testphys1d_mod.F90 (modified) (3 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,q_def_h2o_vap)
+                           play,plev,latitude,longitude,cell_area,prescribed_h2ovap,h2ovap_relax_time,q_prescribed_h2o_vap)
 
 use ioipsl_getincom,          only: getin ! To use 'getin'
@@ -81 +81 @@
 real, dimension(nlayer + 1),         intent(out) :: plev                         ! intermediate pressure levels (pa)
 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
+logical,                             intent(out) :: prescribed_h2ovap            ! Prescribe atmospheric water profile
+real,                                intent(out) :: h2ovap_relax_time            ! Relaxtion time towards the atmospheric water profile
+real, dimension(nlayer),             intent(out) :: q_prescribed_h2o_vap         ! User-defined atmospheric water profile
 
 !=======================================================================
@@ -778 +779 @@
 if (.not. therestartfi .and. paleoclimate) perennial_co2ice = 0.
 
-! Check if the atmospheric water profile is specified
-! ---------------------------------------------------
-! Adding an option to force atmospheric water values
-atm_wat_profile = -1. ! Default: free atmospheric water profile
-q_def_h2o_vap = 0. ! Default
+! Prescription of atmospheric water profile
+! -----------------------------------------
+prescribed_h2ovap = .false. ! Default: free atmospheric water profile
+h2ovap_relax_time = -1. ! Default: no time relaxation
+q_prescribed_h2o_vap = 0. ! Default
 if (water) then
-    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
-    if (abs(atm_wat_profile + 1.) < 1.e-15) then ! if == -1.
+    write(*,*)'Force atmospheric water profile (mixing ratio in kg/kg)?'
+    call getin('prescribed_h2ovap',prescribed_h2ovap)
+    write(*,*) 'prescribed_h2ovap = ', prescribed_h2ovap
+    if (prescribed_h2ovap) then
+        write(*,*) 'Atmospheric water profile is prescribed'
+        write(*,*) 'Unless it reaches saturation (maximal value)'
+        open(10,file = 'profile_prescribed_h2o_vap',status = 'old',form = 'formatted',action = 'read',iostat = ierr)
+        if (ierr == 0) then
+            do ilayer = 1,nlayer
+                read(10,*) q_prescribed_h2o_vap(ilayer)
+            enddo
+        else
+            error stop 'File "profile_prescribed_h2o_vap" was not found!'
+        endif 
+        close(10)
+
+        write(*,*) 'Relax atmospheric water profile with time constant (s)?'
+        
+        call getin('h2ovap_relax_time',h2ovap_relax_time)
+        write(*,*) 'h2ovap_relax_time = ', h2ovap_relax_time
+        if (h2ovap_relax_time < 0.) then
+            write(*,*) 'Atmospheric water vapor profile is not relaxed (forcing).'
+        else
+            write(*,*) 'Atmospheric water profile is relaxed towards the profile in "profile_prescribed_h2o_vap"'
+            write(*,*) 'Unless it reaches saturation (maximal value)'
+        endif
+    else
         write(*,*) 'Free atmospheric water vapor profile'
         write(*,*) 'Total water is conserved in the column'
-    else if (abs(atm_wat_profile) < 1.e-15) then ! if == 0.
-        write(*,*) 'Dry 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 is not correct!'
     endif
 endif
 
-! Check if the atmospheric water profile relaxation is specified
-! --------------------------------------------------------------
-! Adding an option to relax atmospheric water values
-atm_wat_tau = -1. ! Default: no time relaxation
-if (water) then
-    write(*,*) 'Relax atmospheric water vapor profile with time constant (s)?'
-    call getin('atm_wat_tau',atm_wat_tau)
-    write(*,*) 'atm_wat_tau = ', atm_wat_tau
-    if (atm_wat_tau < 0.) then
-        write(*,*) 'Atmospheric water vapor profile is not relaxed (forcing).'
-    else
-        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)'
-        else
-            write(*,*) 'Reference atmospheric water vapor profile not known!'
-            error stop 'Please, specify atm_wat_profile'
-        endif
-    endif
-endif
-
 END SUBROUTINE init_testphys1d
 

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

@@ -84 +84 @@
 logical :: startfiles_1D, therestart1D, therestartfi, there
 
-! Force atmospheric water profiles
-real                            :: atm_wat_profile, atm_wat_tau
-real, dimension(nlayer)         :: q_def_h2o_vap
+! Prescription of atmospheric water profile
+logical                         :: prescribed_h2ovap
+real                            :: h2ovap_relax_time
+real, dimension(nlayer)         :: q_prescribed_h2o_vap
 real, dimension(:), allocatable :: zqsat
 !=======================================================================
@@ -150 +151 @@
 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,q_def_h2o_vap)
+                     play,plev,latitude,longitude,cell_area,prescribed_h2ovap,h2ovap_relax_time,q_prescribed_h2o_vap)
 
 ! Write a "startfi" file
@@ -196 +197 @@
     enddo
 
-    ! Modify atmospheric water if asked
-    ! Added "atm_wat_profile" flag
-    ! ----------------------------
+    ! Prescription of atmospheric water if asked
+    ! ------------------------------------------
     if (water) then
         call watersat(nlayer,temp,play,zqsat)
-        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)
-            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)
+        if (prescribed_h2ovap) then ! If atmospheric water profile is prescribed
+            if (h2ovap_relax_time < 0.) then ! Forcing
+                q(1,:,igcm_h2o_vap) = min(zqsat(:),q_prescribed_h2o_vap(:))
+            else ! Relaxation
+                q(1,:,igcm_h2o_vap) = q_prescribed_h2o_vap(:) + (q(1,:,igcm_h2o_vap) - q_prescribed_h2o_vap(:))*dexp(-dttestphys/h2ovap_relax_time)
             endif
         endif