source: trunk/LMDZ.COMMON/libf/evolution/stopping_crit_mod.F90 @ 3574

Last change on this file since 3574 was 3571, checked in by jbclement, 10 months ago

PEM:

  • New way to manage the pressure: now the PEM manages only the average pressure and keeps the pressure deviation with the instantaneous pressure from the start to reconstruct the pressure at the end ('ps_avg = ps_start + ps_dev'). As a consequence, everything related to pressure in the PEM is modified accordingly.
  • Surface temperatures management is now simpler. It follows the strategy for the pressure (and soil temperature) described above.
  • Soil temperatures are now adapted to match the surface temperature changes occured during the PEM by modifying the soil temperature deviation at the end.
  • Few simplifications/optimizations: notably, the two PCM years are now read in one go in 'read_data_PCM_mod.F90' and only the needed variables are extracted.
  • Deletion of unused variables and unnecessary intermediate variables (memory saving and loop deletion in some cases).
  • Renaming of variables and subroutines to make everything clearer. In particular, the suffixes: '_avg' = average, '_start' = PCM start file, '_dev' = deviation, '_ini' or '0' = initial, '_dyn' = dynamical grid, '_timeseries' = daily average of last PCM year.
  • Cosmetic cleanings for readability.

JBC
File size: 7.0 KB
Line 
1MODULE stopping_crit_mod
2
3implicit none
4
5!=======================================================================
6contains
7!=======================================================================
8
9!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
10!!!
11!!! Purpose: Criterions to check if the PEM needs to call the PCM
12!!! Author: RV & LL, 02/2023
13!!!
14!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15
16SUBROUTINE stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,is_h2oice_sublim_ini,h2o_ice,stopPEM,ngrid)
17
18use time_evol_mod, only: h2o_ice_crit
19use comslope_mod,  only: subslope_dist, nslope
20
21implicit none
22
23!=======================================================================
24!
25! Routine to check if the h2o ice criterion to stop the PEM is reached
26!
27!=======================================================================
28! Inputs
29!-------
30integer,                          intent(in) :: ngrid                ! # of physical grid points
31real,    dimension(ngrid),        intent(in) :: cell_area            ! Area of the cells
32real,    dimension(ngrid,nslope), intent(in) :: h2o_ice              ! Actual density of h2o ice
33real,                             intent(in) :: h2oice_ini_surf      ! Initial surface of sublimating h2o ice
34logical, dimension(ngrid,nslope), intent(in) :: is_h2oice_sublim_ini ! Grid points where h2o ice was initially sublimating
35! Outputs
36!--------
37integer, intent(inout) :: stopPEM ! Stopping criterion code
38! Locals
39! ------
40integer :: i, islope       ! Loop
41real    :: h2oice_now_surf ! Current surface of h2o ice
42
43!=======================================================================
44if (stopPEM > 0) return
45
46! Computation of the present surface of h2o ice still sublimating
47h2oice_now_surf = 0.
48do i = 1,ngrid
49    do islope = 1,nslope
50        if (is_h2oice_sublim_ini(i,islope) .and. h2o_ice(i,islope) > 0.) h2oice_now_surf = h2oice_now_surf + cell_area(i)*subslope_dist(i,islope)
51    enddo
52enddo
53
54! Check of the criterion
55if (h2oice_now_surf < h2oice_ini_surf*(1. - h2o_ice_crit)) then
56    stopPEM = 1
57    write(*,*) "Reason of stopping: the surface of h2o ice sublimating reaches the threshold"
58    write(*,*) "h2oice_now_surf < h2oice_ini_surf*(1. - h2o_ice_crit)", h2oice_now_surf < h2oice_ini_surf*(1. - h2o_ice_crit)
59    write(*,*) "Initial surface of h2o ice sublimating =", h2oice_ini_surf
60    write(*,*) "Current surface of h2o ice sublimating =", h2oice_now_surf
61    write(*,*) "Percentage of change accepted =", h2o_ice_crit*100
62else if (h2oice_now_surf > h2oice_ini_surf*(1. + h2o_ice_crit)) then
63    stopPEM = 1
64    write(*,*) "Reason of stopping: the surface of h2o ice sublimating reaches the threshold"
65    write(*,*) "h2oice_now_surf > h2oice_ini_surf*(1. + h2o_ice_crit)", h2oice_now_surf > h2oice_ini_surf*(1. + h2o_ice_crit)
66    write(*,*) "Initial surface of h2o ice sublimating =", h2oice_ini_surf
67    write(*,*) "Current surface of h2o ice sublimating =", h2oice_now_surf
68    write(*,*) "Percentage of change accepted =", h2o_ice_crit*100
69endif
70
71if (abs(h2oice_ini_surf) < 1.e-5) stopPEM = 0
72
73END SUBROUTINE stopping_crit_h2o_ice
74
75!=======================================================================
76
77SUBROUTINE stopping_crit_co2(cell_area,co2ice_sublim_surf_ini,is_co2ice_sublim_ini,co2_ice,stopPEM,ngrid,ps_avg_global_ini,ps_avg_global,nslope)
78
79use time_evol_mod, only: co2_ice_crit, ps_criterion
80use comslope_mod,  only: subslope_dist
81
82implicit none
83
84!=======================================================================
85!
86! Routine to check if the co2 and pressure criteria to stop the PEM are reached
87!
88!=======================================================================
89! Inputs
90!-------
91integer,                          intent(in) :: ngrid, nslope          ! # of grid physical grid points
92real,    dimension(ngrid),        intent(in) :: cell_area              ! Area of the cells
93real,    dimension(ngrid,nslope), intent(in) :: co2_ice                ! Actual density of co2 ice
94real,                             intent(in) :: co2ice_sublim_surf_ini ! Initial surface of sublimatingco2 ice
95logical, dimension(ngrid,nslope), intent(in) :: is_co2ice_sublim_ini   ! Grid points where co2 ice was initially sublimating
96real,                             intent(in) :: ps_avg_global_ini      ! Planet average pressure from the PCM start files
97real,                             intent(in) :: ps_avg_global          ! Planet average pressure from the PEM computations
98! Outputs
99!--------
100integer, intent(inout) :: stopPEM ! Stopping criterion code
101
102! Locals
103! ------
104integer :: i, islope       ! Loop
105real    :: co2ice_now_surf ! Current surface of co2 ice
106
107!=======================================================================
108if (stopPEM > 0) return
109
110! Computation of the present surface of co2 ice still sublimating
111co2ice_now_surf = 0.
112do i = 1,ngrid
113    do islope = 1,nslope
114        if (is_co2ice_sublim_ini(i,islope) .and. co2_ice(i,islope) > 0.) co2ice_now_surf = co2ice_now_surf + cell_area(i)*subslope_dist(i,islope)
115    enddo
116enddo
117
118! Check of the criterion
119if (co2ice_now_surf < co2ice_sublim_surf_ini*(1. - co2_ice_crit)) then
120    stopPEM = 3
121    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
122    write(*,*) "co2ice_now_surf < co2ice_sublim_surf_ini*(1. - co2_ice_crit)", co2ice_now_surf < co2ice_sublim_surf_ini*(1. - co2_ice_crit)
123    write(*,*) "Initial surface of co2 ice sublimating =", co2ice_sublim_surf_ini
124    write(*,*) "Current surface of co2 ice sublimating =", co2ice_now_surf
125    write(*,*) "Percentage of change accepted =", co2_ice_crit*100.
126else if (co2ice_now_surf > co2ice_sublim_surf_ini*(1. + co2_ice_crit)) then
127    stopPEM = 3
128    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
129    write(*,*) "co2ice_now_surf > co2ice_sublim_surf_ini*(1. + co2_ice_crit)", co2ice_now_surf > co2ice_sublim_surf_ini*(1. + co2_ice_crit)
130    write(*,*) "Current surface of co2 ice sublimating =", co2ice_now_surf
131    write(*,*) "Initial surface of co2 ice sublimating =", co2ice_sublim_surf_ini
132    write(*,*) "Percentage of change accepted =", co2_ice_crit*100.
133endif
134
135if (abs(co2ice_sublim_surf_ini) < 1.e-5) stopPEM = 0
136
137if (ps_avg_global < ps_avg_global_ini*(1. - ps_criterion)) then
138    stopPEM = 4
139    write(*,*) "Reason of stopping: the global pressure reaches the threshold"
140    write(*,*) "ps_avg_global < ps_avg_global_ini*(1. - ps_criterion)", ps_avg_global < ps_avg_global_ini*(1. - ps_criterion)
141    write(*,*) "Initial global pressure =", ps_avg_global_ini
142    write(*,*) "Current global pressure =", ps_avg_global
143    write(*,*) "Percentage of change accepted =", ps_criterion*100.
144else if (ps_avg_global > ps_avg_global_ini*(1. + ps_criterion)) then
145    stopPEM = 4
146    write(*,*) "Reason of stopping: the global pressure reaches the threshold"
147    write(*,*) "ps_avg_global > ps_avg_global_ini*(1. + ps_criterion)", ps_avg_global > ps_avg_global_ini*(1. + ps_criterion)
148    write(*,*) "Initial global pressure =", ps_avg_global_ini
149    write(*,*) "Current global pressure =", ps_avg_global
150    write(*,*) "Percentage of change accepted =", ps_criterion*100.
151endif
152
153END SUBROUTINE stopping_crit_co2
154
155END MODULE
Note: See TracBrowser for help on using the repository browser.