Changeset 3021 for trunk/LMDZ.MARS/libf

Timestamp:

Jul 31, 2023, 5:49:13 PM (16 months ago)

Author:

jbclement

Message:

Mars PCM 1D:
Addition in 1D of relaxation towards a profile for atmospheric water. It needs the flag "atm_wat_tau" (real):

• atm_wat_profile < 0. -> no relaxation (default);
• atm_wat_profile >= 0. -> relaxation towards the value "atm_wat_profile" with relaxation time "atm_wat_tau".

File:

• trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F (modified) (7 diffs)

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

@@ -142 +142 @@
 c   JN : Force atmospheric water profiles
       REAL atm_wat_profile
+      REAL atm_wat_tau
       REAL,ALLOCATABLE :: zqsat(:) ! useful for (atm_wat_profile=2)
 
@@ -813 +814 @@
 c ---------------------------------------------------
       ! Adding an option to force atmospheric water values JN
-      atm_wat_profile=-1 ! Default: free atm wat profile
-      print *,'Force atmospheric water vapor profile ?'
+      atm_wat_profile = -1. ! Default: free atm wat profile
+      print *,'Force atmospheric water vapor profile?'
       call getin("atm_wat_profile",atm_wat_profile)
       write(*,*) "atm_wat_profile = ", atm_wat_profile
-      if (atm_wat_profile.EQ.-1) then
+      if (abs(atm_wat_profile + 1.) < 1.e-15) then ! if == -1.
         write(*,*) "Free atmospheric water vapor profile"
-        write(*,*) "Total water is conserved on the column"
-      else if (atm_wat_profile.EQ.0) then
+        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 ((atm_wat_profile.GT.0).and.(atm_wat_profile.LE.1)) then
-        write(*,*) "Prescribed atmospheric water vapor MMR profile"
+      else if (0. < atm_wat_profile .and. atm_wat_profile <= 1.) then
+        write(*,*) "Prescribed atmospheric water vapor profile"
         write(*,*) "Unless it reaches saturation (maximal value)"
-        write(*,*) "MMR chosen : ", atm_wat_profile
-      endif 
+      else
+        write(*,*) 'Water vapor profile value not correct!'
+        stop
+      endif
+
+! Check if the atmospheric water profile relaxation is specified
+! --------------------------------------------------------------
+      ! Adding an option to relax atmospheric water values JBC
+      atm_wat_tau = -1. ! Default: no time relaxation
+      print*, 'Relax atmospheric water vapor profile?'
+      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"
+      else
+        if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) 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!'
+       write(*,*) 'Please, specify atm_wat_profile'
+       stop
+        endif
+      endif
 
 c  Write a "startfi" file
@@ -881 +905 @@
       ENDDO
 
-!       Force atmospheric water if asked
-!       Added "atm_wat_profile" flag (JN)
-!       -------------------------------------
+!       Modify atmospheric water if asked
+!       Added "atm_wat_profile" flag (JN + JBC)
+!       ---------------------------------------
         call watersat(nlayer,temp,play,zqsat)
         if (0. <= atm_wat_profile .and. atm_wat_profile <= 1.) then
-          q(:,igcm_h2o_vap) = min(zqsat(:),atm_wat_profile)
-!          write(*,*) "prescribed atm_wat_profile (wet if >0, dry if =0)"
-          q(:,igcm_h2o_ice) = 0.
-!          write(*,*) "atm_wat_profile: reset h2o ice"
+        ! If atmospheric water is monitored
+          if (atm_wat_tau < 0.) then ! Prescribed atm_wat_profile
+            ! Wet if >0, Dry if =0
+            q(:,igcm_h2o_vap) = min(zqsat(:),atm_wat_profile*g/psurf)
+            q(:,igcm_h2o_ice) = 0. ! reset h2o ice
+          else ! Relaxation towards the value atm_wat_profile with relaxation time atm_wat_tau
+        q(:,igcm_h2o_vap) = atm_wat_profile*g/psurf + (q(:,igcm_h2o_vap)
+     &              - atm_wat_profile*g/psurf)*dexp(-dtphys/atm_wat_tau)
+            q(:,igcm_h2o_vap) = min(zqsat(:),q(:,igcm_h2o_vap))
+            q(:,igcm_h2o_ice) = 0. ! reset h2o ice
+          endif
         endif
-!        CALL write_output('qsat',
-!     &                         'qsat',
-!     &                         'kg/kg',zqsat)
-!        CALL write_output('qvap',
-!     &                         'qvap',
-!     &                         'kg/kg',q(:,igcm_h2o_vap))
-
 
 !      write(*,*) "testphys1d avant q", q(1,:)
@@ -912 +936 @@
 !      write(*,*) "testphys1d apres q", q(1,:)
 
-
 c       wind increment : specific for 1D
 c       --------------------------------
- 
 c       The physics compute the tendencies on u and v,
 c       here we just add Coriolos effect
@@ -931 +953 @@
           ENDDO
 c       end if
-c      
-c
+
 c       Compute time for next time step
 c       ---------------------------------------
@@ -953 +974 @@
 c       compute pressure for next time step 
 c       ----------------------------------------------------------
-
            psurf=psurf+dtphys*dpsurf(1)   ! surface pressure change
            DO ilevel=1,nlevel
@@ -968 +988 @@
           ENDDO
         ENDDO
-
       ENDDO   ! of idt=1,ndt ! end of time stepping loop
-      
+
       ! update the profiles files at the end of the run
       write_prof=.false.