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

Last change on this file since 3321 was 3206, checked in by jbclement, 10 months ago

PEM:
Cleanings of unused variables/arguments and bad type conversions.
JBC
File size: 6.3 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,surf_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
29!   arguments:
30!   ----------
31!   INPUT
32integer,                       intent(in) :: ngrid     ! # of physical grid points
33real, dimension(ngrid),        intent(in) :: cell_area ! Area of the cells
34real, dimension(ngrid,nslope), intent(in) :: h2o_ice   ! Actual density of h2o ice
35real,                          intent(in) :: surf_ini  ! Initial surface of h2o ice that was sublimating
36!   OUTPUT
37integer, intent(inout) :: stopPEM ! Stopping criterion code
38!   local:
39!   ------
40integer :: i, islope ! Loop
41real    :: surf_now  ! Current surface of h2o ice
42
43!=======================================================================
44! Computation of the present surface of h2o ice
45surf_now = 0.
46do i = 1,ngrid
47    do islope = 1,nslope
48        if (h2o_ice(i,islope) > 0.) surf_now = surf_now + cell_area(i)*subslope_dist(i,islope)
49    enddo
50enddo
51
52! Check of the criterion
53if (surf_now < surf_ini*(1. - h2o_ice_crit)) then
54    stopPEM = 1
55    write(*,*) "Reason of stopping: the surface of h2o ice sublimating reaches the threshold"
56    write(*,*) "surf_now < surf_ini*(1. - h2o_ice_crit)", surf_now < surf_ini*(1. - h2o_ice_crit)
57    write(*,*) "Current surface of h2o ice sublimating =", surf_now
58    write(*,*) "Initial surface of h2o ice sublimating =", surf_ini
59    write(*,*) "Percentage of change accepted =", h2o_ice_crit*100
60else if (surf_now > surf_ini*(1. + h2o_ice_crit)) then
61    stopPEM = 1
62    write(*,*) "Reason of stopping: the surface of h2o ice sublimating reaches the threshold"
63    write(*,*) "surf_now > surf_ini*(1. + h2o_ice_crit)", surf_now > surf_ini*(1. + h2o_ice_crit)
64    write(*,*) "Current surface of h2o ice sublimating =", surf_now
65    write(*,*) "Initial surface of h2o ice sublimating =", surf_ini
66    write(*,*) "Percentage of change accepted =", h2o_ice_crit*100
67endif
68
69if (abs(surf_ini) < 1.e-5) stopPEM = 0
70
71END SUBROUTINE stopping_crit_h2o_ice
72
73!=======================================================================
74
75SUBROUTINE stopping_crit_co2(cell_area,surf_ini,co2_ice,stopPEM,ngrid,global_avg_press_PCM,global_avg_press_new,nslope)
76
77use time_evol_mod, only: co2_ice_crit, ps_criterion
78use comslope_mod,  only: subslope_dist
79
80implicit none
81
82!=======================================================================
83!
84! Routine to check if the co2 and pressure criteria to stop the PEM are reached
85!
86!=======================================================================
87
88!   arguments:
89!   ----------
90!   INPUT
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 h2o ice
94real,                          intent(in) :: surf_ini             ! Initial surface of co2 ice that was sublimating
95real,                          intent(in) :: global_avg_press_PCM ! Planet average pressure from the PCM start files
96real,                          intent(in) :: global_avg_press_new ! Planet average pressure from the PEM computations
97!   OUTPUT
98integer, intent(inout) :: stopPEM ! Stopping criterion code
99
100!   local:
101!   ------
102integer :: i, islope ! Loop
103real    :: surf_now  ! Current surface of co2 ice
104
105!=======================================================================
106! Computation of the present surface of co2 ice
107surf_now = 0.
108do i = 1,ngrid
109    do islope = 1,nslope
110        if (co2_ice(i,islope) > 0.) surf_now = surf_now + cell_area(i)*subslope_dist(i,islope)
111    enddo
112enddo
113
114! Check of the criterion
115if (surf_now < surf_ini*(1. - co2_ice_crit)) then
116    stopPEM = 3
117    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
118    write(*,*) "surf_now < surf_ini*(1. - co2_ice_crit)", surf_now < surf_ini*(1. - co2_ice_crit)
119    write(*,*) "Current surface of co2 ice sublimating =", surf_now
120    write(*,*) "Initial surface of co2 ice sublimating =", surf_ini
121    write(*,*) "Percentage of change accepted =", co2_ice_crit*100.
122else if (surf_now > surf_ini*(1. + co2_ice_crit)) then
123    stopPEM = 3
124    write(*,*) "Reason of stopping: the surface of co2 ice sublimating reaches the threshold"
125    write(*,*) "surf_now > surf_ini*(1. + co2_ice_crit)", surf_now > surf_ini*(1. + co2_ice_crit)
126    write(*,*) "Current surface of co2 ice sublimating =", surf_now
127    write(*,*) "Initial surface of co2 ice sublimating =", surf_ini
128    write(*,*) "Percentage of change accepted =", co2_ice_crit*100.
129endif
130
131if (abs(surf_ini) < 1.e-5) stopPEM = 0
132
133if (global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)) then
134    stopPEM = 4
135    write(*,*) "Reason of stopping: the global pressure reaches the threshold"
136    write(*,*) "global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)", global_avg_press_new < global_avg_press_PCM*(1. - ps_criterion)
137    write(*,*) "Current global pressure =", global_avg_press_new
138    write(*,*) "PCM global pressure =", global_avg_press_PCM
139    write(*,*) "Percentage of change accepted =", ps_criterion*100.
140else if (global_avg_press_new > global_avg_press_PCM*(1. + ps_criterion)) then
141    stopPEM = 4
142    write(*,*) "Reason of stopping: the global pressure reaches the threshold"
143    write(*,*) "global_avg_press_new > global_avg_press_PCM*(1. + ps_criterion)", global_avg_press_new > global_avg_press_PCM*(1. + ps_criterion)
144    write(*,*) "Current global pressure =", global_avg_press_new
145    write(*,*) "PCM global pressure =", global_avg_press_PCM
146    write(*,*) "Percentage of change accepted =", ps_criterion*100.
147endif
148
149END SUBROUTINE stopping_crit_co2
150
151END MODULE
Note: See TracBrowser for help on using the repository browser.