Changeset 2857 for trunk/LMDZ.COMMON/libf/evolution

Timestamp:

Dec 29, 2022, 10:58:39 AM (2 years ago)

Author:

llange

Message:

PEM
Fixing bug when updating tracer mass mixing ratio (ccn_number and dust_number were previously set to 1)
Tiny edits in evol_h2o_ice and criterion_ice_stop_mod_water
LL

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• criterion_ice_stop_mod_water_slope.F90 (modified) (1 diff)
• evol_h2o_ice_s_mod_slope.F90 (modified) (3 diffs)
• pem.F90 (modified) (11 diffs)

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

@@ -62 +62 @@
     STOPPING=.FALSE.
   endif
-
 END SUBROUTINE criterion_ice_stop_water_slope
 

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

@@ -26 +26 @@
 !   OUTPUT
   REAL, INTENT(INOUT) ::  qsurf(ngrid,nslope)                ! physical point field : Previous and actual density of water ice
-  LOGICAL :: STOPPING
+  LOGICAL, INTENT(INOUT) :: STOPPING
   REAL, intent(inout) ::  tendencies_h2o_ice_phys(ngrid,nslope) ! physical point field : Evolution of perenial ice over one year
 
@@ -80 +80 @@
     print *, "Tendencies on ice increasing=", pos_tend
     print *, "This can be due to the absence of water ice in the PCM run!!"
-      call criterion_ice_stop_water_slope(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
@@ -101 +101 @@
     enddo
   enddo
+  
 
   real_coefficient=negative_part/pos_tend
+  if(pos_tend.eq.0) real_coefficient = 0.
 
   do i=1,ngrid

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

@@ -843 +843 @@
 !------------------------
      year_iter=0
-     print *, "Max_iter_pem", Max_iter_pem
-     print *, 'year_iter_max',year_iter_max
      do while (year_iter.LT.year_iter_max)
 
@@ -852 +850 @@
        do islope=1,nslope
            global_ave_press_new=global_ave_press_new-g*cell_area(i)*tendencies_co2_ice_phys_slope(i,islope)*subslope_dist(i,islope)/cos(pi*def_slope_mean(islope)/180.)/Total_surface
-       enddo
+      enddo
      enddo
        print *, 'Global average pressure old time step',global_ave_press_old
@@ -932 +930 @@
      call evol_h2o_ice_s_slope(qsurf_slope(:,igcm_h2o_ice,:),tendencies_h2o_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_water,nslope)
 
+
       print *, "Evolution of co2 ice"
      call evol_co2_ice_s_slope(co2ice_slope,tendencies_co2_ice_phys_slope,iim,jjm,ngrid,cell_area,STOPPING_1_co2,nslope)
@@ -972 +971 @@
       DO ig = 1,ngrid
         DO islope = 1,nslope
-          if(initial_co2_ice_slope(ig,islope).gt.0.5 .and. co2ice_slope(ig,islope).LT. 1E-5) THEN !co2ice disappeared, look for closest point without co2ice
+          if(initial_co2_ice_slope(ig,islope).gt.0.5 .and. co2ice_slope(ig,islope).LT. 1E-10) THEN !co2ice disappeared, look for closest point without co2ice
             if(latitude_deg(ig).gt.0) then
               do ig_loop=ig,ngrid
                 DO islope_loop=islope,iflat,-1
-                  if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-5) then
+                  if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then
                     tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
                     bool_sublim=1
@@ -989 +988 @@
               do ig_loop=ig,1,-1
                 DO islope_loop=islope,iflat
-                  if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-5) then
+                  if(initial_co2_ice_slope(ig_loop,islope_loop).lt.0.5 .and. co2ice_slope(ig_loop,islope_loop).LT. 1E-10) then
                     tsurf_ave_phys(ig,islope)=tsurf_ave_phys(ig_loop,islope_loop)
                     bool_sublim=1
@@ -1001 +1000 @@
             endif
             initial_co2_ice_slope(ig,islope)=0
-            if ((co2ice_slope(ig,islope).lt.1e-5).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold))  then   
+            if ((co2ice_slope(ig,islope).lt.1e-10).and. (qsurf_slope(ig,igcm_h2o_ice,islope).gt.frost_albedo_threshold))  then   
               albedo_slope(ig,1,islope) = albedo_h2o_frost
               albedo_slope(ig,2,islope) = albedo_h2o_frost
@@ -1009 +1008 @@
               emiss_slope(ig,islope) = emissiv         
             endif
-          elseif( (co2ice_slope(ig,islope).GT. 1E-3).and.(tendencies_co2_ice_phys_slope(ig,islope).gt.1e-5) )THEN !Put tsurf as tcond co2
+          elseif( (co2ice_slope(ig,islope).GT. 1E-3).and.(tendencies_co2_ice_phys_slope(ig,islope).gt.1e-10) )THEN !Put tsurf as tcond co2
             ave=0.
             do t=1,timelen
@@ -1065 +1064 @@
           if(isnan(Tsoil_locslope(ig,isoil))) then
              write(*,*)'Tsoil=',ig,isoil,Tsoil_locslope(ig,isoil)
+             stop 
           endif
          enddo
@@ -1148 +1148 @@
         print *, "STOPPING because tendencies on co2 ice=0, see message above", STOPPING_1_co2
       endif
+
+
 
       if (STOPPING_water .or. STOPPING_1_water .or. STOPPING_co2 .or. STOPPING_1_co2)  then
@@ -1260 +1262 @@
      enddo
 
+     write(*,*) 'rapport ps',global_ave_press_new/global_ave_press_GCM
+
 ! III_a.5 Tracer (for start)
      allocate(zplev_new(ngrid,nlayer+1))
@@ -1294 +1298 @@
        DO ig=1,ngrid
          DO l=1,llm-1
-           if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number" .or. noms(nnq).NE."ccn_number") ) then
+           if(q(ig,l,nnq).GT.1 .and. (noms(nnq).NE."dust_number") .and. (noms(nnq).NE."ccn_number") ) then
              extra_mass=(q(ig,l,nnq)-1)*(zplev_new(ig,l)-zplev_new(ig,l+1))
              q(ig,l,nnq)=1.
              q(ig,l+1,nnq)=q(ig,l+1,nnq)+extra_mass*(zplev_new(ig,l+1)-zplev_new(ig,l+2))
+             write(*,*) 'extra ',noms(nnq),extra_mass, noms(nnq).NE."dust_number",noms(nnq).NE."ccn_number"
            endif
            if(q(ig,l,nnq).LT.0) then