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

Last change on this file since 3532 was 3432, checked in by jbclement, 3 months ago

PEM:
Small corrections for the launching script (run numbers must be integer) and the PEM stopping criteria algorithm related to r3430.
JBC

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