Changeset 3070

Timestamp:

Oct 5, 2023, 11:29:45 AM (14 months ago)

Author:

jbclement

Message:

PEM:
Correction of a bug: the variable 'g' was not correctly initialized in 1D. A little of code cleaning.
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• compute_tendencies_slope.F90 (modified) (2 diffs)
• pem.F90 (modified) (17 diffs)

trunk/LMDZ.COMMON/libf/evolution/changelog.txt

@@ -86 +86 @@
 == 03/10/2023 == JBC
 Following the commits r3066 and r3067, the PEM initialization in 1D has been adapted.
+
+== 05/10/2023 == JBC
+Correction of a bug: the variable 'g' was not correctly initialized in 1D. A little of code cleaning.

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

@@ -1 @@
-!
-! $Id $
-!
-SUBROUTINE compute_tendencies_slope(ngrid,nslope,min_ice_Y1,&
-     min_ice_Y2,tendencies_ice)
+SUBROUTINE compute_tendencies_slope(ngrid,nslope,min_ice_Y1,min_ice_Y2,tendencies_ice)
 
-      IMPLICIT NONE
+implicit none
 
 !=======================================================================
@@ -15 +11 @@
 !   arguments:
 !   ----------
-
 !   INPUT
-
-     INTEGER, intent(in) :: ngrid, nslope                           ! # of grid points along longitude/latitude/ total
-     REAL, intent(in) , dimension(ngrid,nslope):: min_ice_Y1       ! LON x LAT field : minimum of water ice at each point for the first year
-     REAL, intent(in) , dimension(ngrid,nslope):: min_ice_Y2       ! LON x LAT field : minimum of water ice at each point for the second year
+integer,                       intent(in) :: ngrid, nslope ! # of grid points along longitude/latitude/ total
+real, dimension(ngrid,nslope), intent(in) :: min_ice_Y1    ! LON x LAT field : minimum of water ice at each point for the first year
+real, dimension(ngrid,nslope), intent(in) :: min_ice_Y2    ! LON x LAT field : minimum of water ice at each point for the second year
 
 !   OUTPUT
-     REAL, intent(out) , dimension(ngrid,nslope)   :: tendencies_ice            ! physical point field : difference between the minima = evolution of perenial ice
-
-!   local:
-!   ------
-     INTEGER :: ig,islope                                                           ! loop variable
+real, dimension(ngrid,nslope), intent(out) :: tendencies_ice ! physical point field : difference between the minima = evolution of perenial ice
 
 !=======================================================================
 
 !  We compute the difference
-
-  DO ig=1,ngrid
-    DO islope = 1, nslope
-      tendencies_ice(ig,islope)=min_ice_Y2(ig,islope)-min_ice_Y1(ig,islope)
-    enddo
-  ENDDO
+tendencies_ice = min_ice_Y2 - min_ice_Y1
 
 !  If the difference is too small; there is no evolution
-  DO ig=1,ngrid
-    DO islope = 1, nslope
-      if(abs(tendencies_ice(ig,islope)).LT.1.0E-10) then
-        tendencies_ice(ig,islope)=0.
-      endif
-    enddo
-  ENDDO
+where (abs(tendencies_ice) < 1.e-10) tendencies_ice = 0.
 
 END SUBROUTINE compute_tendencies_slope

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

@@ -34 +34 @@
 use logic_mod,                  only: iflag_phys
 use mod_const_mpi,              only: COMM_LMDZ
-use comconst_mod,               only: rad, g, cpp, pi
 use infotrac
 use geometry_mod,               only: latitude_deg
@@ -73 +72 @@
     use mod_phys_lmdz_para, only: is_parallel, is_sequential, is_mpi_root, is_omp_root, is_master
     use planete_h,          only: aphelie, periheli, year_day, peri_day, obliquit
-    use comcstfi_h,         only: r, mugaz
+    use comcstfi_h,         only: r, mugaz, g
     use paleoclimate_mod,   only: albedo_perenialco2
 #else
@@ -88 +87 @@
     use iniphysiq_mod,    only: iniphysiq
     use control_mod,      only: iphysiq, day_step, nsplit_phys
+    use comconst_mod,     only: rad, g, cpp, pi
 #else
     use time_phylmdz_mod, only: iphysiq, day_step
+    use comcstfi_h,       only: pi, rad, g, cpp
 #endif
 
@@ -339 +340 @@
 ! In the gcm, these values are given to the physic by the dynamic.
 ! Here we simply read them in the startfi_evol.nc file
-status =nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
+status = nf90_open(FILE_NAME, NF90_NOWRITE, ncid)
 
 status = nf90_inq_varid(ncid,"longitude",lonvarid)
@@ -362 +363 @@
                   watercap,perenial_co2ice,def_slope,def_slope_mean,subslope_dist)
 
-! Remove unphysical values of surface tracer
-    do i = 1,ngrid
-        do nnq = 1,nqtot
-            do islope = 1,nslope
-                if (qsurf(i,nnq,islope) < 0) qsurf(i,nnq,islope) = 0.
-            enddo
-        enddo
-    enddo
+    ! Remove unphysical values of surface tracer
+    where (qsurf < 0.) qsurf = 0.
 
     call surfini(ngrid,qsurf)
@@ -416 +411 @@
     endif
 
-! Remove unphysical values of surface tracer
-    do i = 1,ngrid
-        do nnq = 1,nqtot
-            qsurf(i,nnq,1) = qsurf_read_generic(i,nnq)
-            if (qsurf(i,nnq,1) < 0) qsurf(i,nnq,1) = 0.
-        enddo
-    enddo
+    ! Remove unphysical values of surface tracer
+    qsurf(:,:,1) = qsurf_read_generic(:,:)
+    where (qsurf < 0.) qsurf = 0.
 #endif
 
@@ -517 +508 @@
 if (soil_pem) then
     call soil_settings_PEM(ngrid,nslope,nsoilmx_PEM,nsoilmx,inertiesoil,TI_PEM)
-    do l=1,nsoilmx
-        tsoil_PEM(:,l,:)=tsoil_ave(:,l,:)
-        tsoil_phys_PEM_timeseries(:,l,:,:)=tsoil_GCM_timeseries(:,l,:,:)
-        watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,l,:,:)
-    enddo
-    do l=nsoilmx+1,nsoilmx_PEM
-        tsoil_PEM(:,l,:)=tsoil_ave(:,nsoilmx,:)
-        watersoil_density_PEM_timeseries(:,l,:,:)=watersoil_density_timeseries(:,nsoilmx,:,:)
-    enddo
-    watersoil_density_PEM_ave(:,:,:) = SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
+    do l = 1,nsoilmx
+        tsoil_PEM(:,l,:) = tsoil_ave(:,l,:)
+        tsoil_phys_PEM_timeseries(:,l,:,:) = tsoil_GCM_timeseries(:,l,:,:)
+        watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,l,:,:)
+    enddo
+    do l = nsoilmx + 1,nsoilmx_PEM
+        tsoil_PEM(:,l,:) = tsoil_ave(:,nsoilmx,:)
+        watersoil_density_PEM_timeseries(:,l,:,:) = watersoil_density_timeseries(:,nsoilmx,:,:)
+    enddo
+    watersoil_density_PEM_ave(:,:,:) = sum(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
 endif !soil_pem
 deallocate(tsoil_ave,tsoil_GCM_timeseries)
@@ -559 +550 @@
 Total_surface = 0.
 do i = 1,ngrid
-    Total_surface = Total_surface+cell_area(i)
+    Total_surface = Total_surface + cell_area(i)
     do islope = 1,nslope
         if (tendencies_co2_ice(i,islope) < 0) then
             initial_co2_ice_sublim(i,islope) = 1.
-            ini_surf_co2 = ini_surf_co2+cell_area(i)*subslope_dist(i,islope)
+            ini_surf_co2 = ini_surf_co2 + cell_area(i)*subslope_dist(i,islope)
         else
             initial_co2_ice_sublim(i,islope) = 0.
@@ -586 +577 @@
 allocate(zplev_gcm(ngrid,nlayer + 1))
 
-do l = 1,nlayer + 1
-    do ig = 1,ngrid
-        zplev_gcm(ig,l) = ap(l) + bp(l)*ps_start_GCM(ig)
-    enddo
+do ig = 1,ngrid
+    zplev_gcm(ig,:) = ap(:) + bp(:)*ps_start_GCM(ig)
 enddo
 
@@ -623 +612 @@
     totmassh2o_adsorbded = 0.
     do ig = 1,ngrid
-        do islope = 1, nslope
+        do islope = 1,nslope
             do l = 1,nsoilmx_PEM - 1
                 totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
@@ -665 +654 @@
 ! II.a.1. Compute updated global pressure
     write(*,*) "Recomputing the new pressure..."
-
     do i = 1,ngrid
         do islope = 1,nslope
@@ -682 +670 @@
 
 ! II.a.2. Old pressure levels for the timeseries, this value is deleted when unused and recreated each time (big memory consuption)
-    allocate(zplev_old_timeseries(ngrid,nlayer+1,timelen))
+    allocate(zplev_old_timeseries(ngrid,nlayer + 1,timelen))
     write(*,*) "Recomputing the old pressure levels timeserie adapted to the old pressure..."
     do l = 1,nlayer + 1
@@ -721 +709 @@
         do t = 1,timelen
             q_co2_PEM_phys(ig,t) = q_co2_PEM_phys(ig,t)*(zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t))/ &
-                                   (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))  &
-                                + ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))  &
+                                   (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))   &
+                                + ((zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))   &
                                 -  (zplev_old_timeseries(ig,l,t) - zplev_old_timeseries(ig,l + 1,t)))/ &
                                    (zplev_new_timeseries(ig,l,t) - zplev_new_timeseries(ig,l + 1,t))
@@ -881 +869 @@
             enddo
         enddo
-        tsoil_PEM(:,:,:) = SUM(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
-        watersoil_density_PEM_ave(:,:,:)= SUM(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
+        tsoil_PEM(:,:,:) = sum(tsoil_phys_PEM_timeseries(:,:,:,:),4)/timelen
+        watersoil_density_PEM_ave(:,:,:) = sum(watersoil_density_PEM_timeseries(:,:,:,:),4)/timelen
 
         write(*,*) "Update of soil temperature done"
@@ -893 +881 @@
         call computeice_table_equilibrium(ngrid,nslope,nsoilmx_PEM,watercaptag,watersurf_density_ave,watersoil_density_PEM_ave,TI_PEM(:,1,:),porefillingice_depth,porefillingice_thickness)
 
-        call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsoil_PEM, &
-                                          delta_h2o_icetablesublim)        ! Mass of H2O exchange between the ssi and the atmosphere
+        call compute_massh2o_exchange_ssi(ngrid,nslope,nsoilmx_PEM,porefillingice_thickness_prev_iter,porefillingice_thickness,porefillingice_depth,tsoil_PEM,delta_h2o_icetablesublim) ! Mass of H2O exchange between the ssi and the atmosphere
 
         write(*,*) "Update soil propreties"
 
 ! II_d.4 Update the soil thermal properties
-        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new, &
-                                           porefillingice_depth,porefillingice_thickness,TI_PEM)
+        call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,tendencies_h2o_ice,qsurf(:,igcm_h2o_ice,:),global_ave_press_new,porefillingice_depth,porefillingice_thickness,TI_PEM)
 
 ! II_d.5 Update the mass of the regolith adsorbded
@@ -914 +900 @@
                     do l = 1,nsoilmx_PEM - 1
                         totmassco2_adsorbded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
-                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
-                        cell_area(ig)
+                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
                         totmassh2o_adsorbded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l+1) - layer_PEM(l))* &
-                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.) * &
-                        cell_area(ig)
+                        subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig)
                     enddo
                 enddo
@@ -1077 +1061 @@
             do ig = 1,ngrid
                 q(ig,l,nnq) = q(ig,l,nnq)*(zplev_gcm(ig,l) - zplev_gcm(ig,l + 1))/(zplev_new(ig,l) - zplev_new(ig,l + 1)) &
-                              + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_gcm(ig,l) - zplev_gcm(ig,l + 1)))/      &
-                              (zplev_new(ig,l) - zplev_new(ig,l + 1))
+                              + ((zplev_new(ig,l) - zplev_new(ig,l + 1)) - (zplev_gcm(ig,l) - zplev_gcm(ig,l + 1)))/(zplev_new(ig,l) - zplev_new(ig,l + 1))
             enddo
             q(:,llm,nnq) = q(:,llm - 1,nnq)