Changeset 2794 for trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90

Timestamp:

Sep 9, 2022, 4:02:51 PM (2 years ago)

Author:

llange

Message:

MARS PEM:

• Add a PEMETAT0 that read "startfi_pem.nc"
• Add the soil in the model: soil temperature, thermal properties, ice table
• Add a routine that compute CO2 + H2O adsorption
• Minor corrections in PEM.F90

LL

File:

• trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90 (modified) (9 diffs)

trunk/LMDZ.COMMON/libf/evolution/evol_h2o_ice_s_mod_slope.F90

@@ -6 +6 @@
 
   USE temps_mod_evol, ONLY: dt_pem
-  use comslope_mod,   ONLY: subslope_dist
+          use comslope_mod, ONLY: subslope_dist
 
       IMPLICIT NONE
 
 !=======================================================================
-!  
-!  Routine that compute the evolution of the water ice.
-!  
-!  We suppose that the water exchange uniquely form one point to another 
-!  (water ice disappear from a sublimating point and increase at another 
-!  accumulation point). Therfor the total positive and negative tendencies 
-!  must be equal to conserve the total amount of water
-!  
+!
+!  Routine that compute the evolution of the water ice
+!
 !=======================================================================
 
@@ -26 +21 @@
 !   INPUT
 
-  INTEGER, intent(in) :: iim_input,jjm_input, ngrid,nslope      ! # of grid points along longitude/latitude/ total
-  REAL, intent(in) ::  cell_area(ngrid)                         ! Area of each cell
+  INTEGER, intent(in) :: iim_input,jjm_input, ngrid,nslope   ! # of grid points along longitude/latitude/ total
+!  REAL, intent(in) ::  tendencies_h2o_ice_phys(ngrid) ! physical point field : Evolution of perenial ice over one year
+  REAL, intent(in) ::  cell_area(ngrid)
 
 !   OUTPUT
-  REAL, INTENT(INOUT) ::  qsurf(ngrid,nslope)                   ! physical point field : Previous and actual density of water ice
-  LOGICAL :: STOPPING                                           ! Boolean to alert that there is no sublimating point anymore
+  REAL, INTENT(INOUT) ::  qsurf(ngrid,nslope)                ! physical point field : Previous and actual density of water ice
+  LOGICAL :: STOPPING
   REAL, intent(inout) ::  tendencies_h2o_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
 
@@ -39 +35 @@
 
   INTEGER :: i,j,ig0,islope                                  ! loop variable
-  REAL :: pos_tend, neg_tend, real_coefficient               ! The integral of positive/negative tendencies overall the planet. Ratio
-  REAL :: negative_part                                      ! After applying the tendencies, some point can be negative, we sum all of them in this variable
-  REAL :: new_tendencies(ngrid,nslope)                       ! Tendencies after taking into account the non physical negative part above
+!  REAL :: not_budget, budget
+  REAL :: pos_tend, neg_tend, real_coefficient,negative_part
+  REAL ::  new_tendencies(ngrid,nslope)
 
 
 !=======================================================================
 
-! Integral of the positive and negative tendencies
+!  budget=sum(qsurf(:))
+
   pos_tend=0.
   neg_tend=0.
@@ -62 +59 @@
   enddo
 
-! These two quantities must be equal to conserve water, therfor we apply a ratio to them
+!  print *, "pos_tend", pos_tend
+!  print *, "neg_tend", neg_tend
 
   if(neg_tend.GT.pos_tend .and. pos_tend.GT.0) then
@@ -68 +66 @@
        do islope=1,nslope
        if(tendencies_h2o_ice_phys(i,islope).LT.0) then
+!          print *, "pos_tend/neg_tend", pos_tend/neg_tend
           new_tendencies(i,islope)=tendencies_h2o_ice_phys(i,islope)*(pos_tend/neg_tend)
        else
@@ -75 +74 @@
      enddo
   elseif(neg_tend.LT.pos_tend .and. neg_tend.GT.0) then
+!          print *, "neg_tend/pos_tend", neg_tend/pos_tend
      do i=1,ngrid
        do islope=1,nslope
@@ -85 +85 @@
      enddo
   elseif(pos_tend.EQ.0 .OR. neg_tend.EQ.0) then 
-      call criterion_ice_stop_water_slope(cell_area,1,qsurf(:,:)*0,STOPPING,ngrid,cell_area)
+
+!      call criterion_ice_stop(cell_area,1,qsurf*0.,STOPPING,ngrid,cell_area)
+      call criterion_ice_stop_water_slope(cell_area,1,qsurf(:,:)*0.,STOPPING,ngrid,cell_area)
       do i=1,ngrid
          do islope=1,nslope
@@ -93 +95 @@
   endif
 
+ 
+  
+  negative_part = 0.
+
 
 ! Evolution of the water ice for each physical point
   do i=1,ngrid
     do islope=1, nslope
+!      qsurf(i)=qsurf(i)+tendencies_h2o_ice_phys(i)*dt_pem
       qsurf(i,islope)=qsurf(i,islope)+new_tendencies(i,islope)*dt_pem
-      if (qsurf(i,islope).lt.0) then         ! If we have negative number, we put them to zero add sum them
+!      budget=budget+tendencies_h2o_ice_phys(i)*dt_pem
+      if (qsurf(i,islope).lt.0) then
+!        not_budget=not_budget+qsurf(i)
+!       print *, "NNqsurf(i,islope)", qsurf(i,islope)
+!       print *, "NNnew_tendencies(i,islope)", new_tendencies(i,islope)
+!       print *, "NNtendencies_h2o_ice_phys(i,islope)", tendencies_h2o_ice_phys(i,islope)
         negative_part=negative_part-qsurf(i,islope)*cell_area(i)*subslope_dist(i,islope)
         qsurf(i,islope)=0.
         tendencies_h2o_ice_phys(i,islope)=0.
+!        print *, "NNineg", i
+      endif
+      if(qsurf(i,islope).NE.qsurf(i,islope)) then
+!          print *, "qsurf(i,islope)",qsurf(i,islope)
+!          print *, "new_tendencies",new_tendencies(i,islope)
+!          print *, "tendencies_h2o_ice_phys",tendencies_h2o_ice_phys(i,islope)
+!          print *, "i", i
+!          print *,"islope",islope
       endif
     enddo
   enddo
 
-! This negative part must be put somewhere else to conserve water, therfor we compute new tendencies
+!  print *, "negative_part", negative_part
   real_coefficient=negative_part/pos_tend
+!  print *, "real_coefficient", real_coefficient
 
   do i=1,ngrid
@@ -118 +139 @@
 
 
+
+! Conservation of water ice
+!  qsurf(:)=qsurf(:)*budget/sum(qsurf(:))
+
+
 END SUBROUTINE evol_h2o_ice_s_slope