source: trunk/LMDZ.COMMON/libf/evolution/surface.F90 @ 4076

Last change on this file since 4076 was 4074, checked in by jbclement, 11 days ago

PEM:

  • Correct management of H2O ice tendency in 1D when there is not enough ice anymore.
  • Clean initialization of allocatable module arrays (especially needed when no slope)
  • One more renaming for consistency + few small updates thoughout the code.

JBC
File size: 11.1 KB
Line 
1MODULE surface
2!-----------------------------------------------------------------------
3! NAME
4!     surface
5!
6! DESCRIPTION
7!     Contains global parameters used for the surface.
8!
9! AUTHORS & DATE
10!     JB Clement, 12/2025
11!
12! NOTES
13!
14!-----------------------------------------------------------------------
15
16! DEPENDENCIES
17! ------------
18use numerics, only: dp, di, k4
19
20! DECLARATION
21! -----------
22implicit none
23
24! PARAMETERS
25! ----------
26real(dp), dimension(:),   allocatable, protected :: albedodat_PCM  ! Albedo of bare ground
27real(dp), dimension(:,:), allocatable, protected :: albedo_PCM     ! Surface albedo_PCM
28real(dp), dimension(:,:), allocatable, protected :: emissivity_PCM ! Thermal IR surface emissivity_PCM
29
30contains
31!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32
33!=======================================================================
34SUBROUTINE ini_surface()
35!-----------------------------------------------------------------------
36! NAME
37!     ini_surface
38!
39! DESCRIPTION
40!     Initialize the parameters of module 'surface'.
41!
42! AUTHORS & DATE
43!     JB Clement, 12/2025
44!
45! NOTES
46!
47!-----------------------------------------------------------------------
48
49! DEPENDENCIES
50! ------------
51use geometry, only: ngrid, nslope
52
53! DECLARATION
54! -----------
55implicit none
56
57! CODE
58! ----
59if (.not. allocated(albedo_PCM)) allocate(albedo_PCM(ngrid,nslope))
60if (.not. allocated(albedodat_PCM)) allocate(albedodat_PCM(ngrid))
61if (.not. allocated(emissivity_PCM)) allocate(emissivity_PCM(ngrid,nslope))
62albedo_PCM(:,:) = 0._dp
63albedodat_PCM(:) = 0._dp
64emissivity_PCM(:,:) = 1._dp
65
66END SUBROUTINE ini_surface
67!=======================================================================
68
69!=======================================================================
70SUBROUTINE end_surface()
71!-----------------------------------------------------------------------
72! NAME
73!     end_surface
74!
75! DESCRIPTION
76!     Deallocate surface arrays.
77!
78! AUTHORS & DATE
79!     JB Clement, 12/2025
80!
81! NOTES
82!
83!-----------------------------------------------------------------------
84
85! DECLARATION
86! -----------
87implicit none
88
89! CODE
90! ----
91if (allocated(albedo_PCM)) deallocate(albedo_PCM)
92if (allocated(albedodat_PCM)) deallocate(albedodat_PCM)
93if (allocated(emissivity_PCM)) deallocate(emissivity_PCM )
94
95END SUBROUTINE end_surface
96!=======================================================================
97
98!=======================================================================
99SUBROUTINE set_albedodat_PCM(albedodat_PCM_in)
100!-----------------------------------------------------------------------
101! NAME
102!     set_albedodat_PCM
103!
104! DESCRIPTION
105!     Setter for 'albedodat_PCM'.
106!
107! AUTHORS & DATE
108!     JB Clement, 12/2025
109!
110! NOTES
111!
112!-----------------------------------------------------------------------
113
114! DECLARATION
115! -----------
116implicit none
117
118! ARGUMENTS
119! ---------
120real(dp), dimension(:), intent(in) :: albedodat_PCM_in
121
122! CODE
123! ----
124albedodat_PCM(:) = albedodat_PCM_in(:)
125
126END SUBROUTINE set_albedodat_PCM
127!=======================================================================
128
129!=======================================================================
130SUBROUTINE set_albedo_PCM(albedo_PCM_in)
131!-----------------------------------------------------------------------
132! NAME
133!     set_albedo_PCM
134!
135! DESCRIPTION
136!     Setter for 'albedo_PCM'.
137!
138! AUTHORS & DATE
139!     JB Clement, 12/2025
140!
141! NOTES
142!
143!-----------------------------------------------------------------------
144
145! DECLARATION
146! -----------
147implicit none
148
149! ARGUMENTS
150! ---------
151real(dp), dimension(:,:), intent(in) :: albedo_PCM_in
152
153! CODE
154! ----
155albedo_PCM(:,:) = albedo_PCM_in(:,:)
156
157END SUBROUTINE set_albedo_PCM
158!=======================================================================
159
160!=======================================================================
161SUBROUTINE set_emissivity_PCM(emissivity_PCM_in)
162!-----------------------------------------------------------------------
163! NAME
164!     set_emissivity_PCM
165!
166! DESCRIPTION
167!     Setter for 'emissivity_PCM'.
168!
169! AUTHORS & DATE
170!     JB Clement, 12/2025
171!
172! NOTES
173!
174!-----------------------------------------------------------------------
175
176! DECLARATION
177! -----------
178implicit none
179
180! ARGUMENTS
181! ---------
182real(dp), dimension(:,:), intent(in) :: emissivity_PCM_in
183
184! CODE
185! ----
186emissivity_PCM(:,:) = emissivity_PCM_in(:,:)
187
188END SUBROUTINE set_emissivity_PCM
189!=======================================================================
190
191!=======================================================================
192SUBROUTINE build4PCM_surf_rad_prop(h2o_ice,co2_ice,albedo4PCM,emissivity4PCM)
193!-----------------------------------------------------------------------
194! NAME
195!     build4PCM_surf_rad_prop
196!
197! DESCRIPTION
198!     Reconstructs albedo and emissivity for the PCM.
199!
200! AUTHORS & DATE
201!     JB Clement, 12/2025
202!     C. Metz, 01/2026
203!
204! NOTES
205!
206!-----------------------------------------------------------------------
207
208! DEPENDENCIES
209! ------------
210use geometry, only: ngrid, nslope, latitudes
211use frost,    only: h2o_frost4PCM, co2_frost4PCM
212use display,  only: print_msg
213
214! DECLARATION
215! -----------
216implicit none
217
218! ARGUMENTS
219! ---------
220real(dp), dimension(:,:), intent(in)  :: h2o_ice, co2_ice
221real(dp), dimension(:,:), intent(out) :: albedo4PCM, emissivity4PCM
222
223! LOCAL VARIABLES
224! ---------------
225integer(di)            :: i, islope, icap
226real(dp)               :: albedo_h2oice, albedo_h2ofrost, h2ofrost_albedo_threshold, emissivity_baresoil
227real(dp), dimension(2) :: albedo_co2frost, emissivity_co2ice, albedo_co2ice
228
229! CODE
230! ----
231! Fetch parameters from "callphys.def" and "startfi.nc"
232call read_albedo_emis(albedo_h2oice,albedo_h2ofrost,h2ofrost_albedo_threshold, &
233                      albedo_co2ice,albedo_co2frost,emissivity_co2ice,emissivity_baresoil)
234
235! Reconstruction Loop
236call print_msg('> Building albedo and emmissivity for the PCM')
237do i = 1,ngrid
238    ! Determine hemisphere: 1 = Northern, 2 = Southern
239    if (latitudes(i) < 0._dp) then
240        icap = 2
241    else
242        icap = 1
243    end if
244
245    do islope = 1,nslope
246        ! Bare ground (default)
247        albedo4PCM(i,islope) = albedodat_PCM(i)
248        emissivity4PCM(i,islope) = emissivity_baresoil
249
250        ! H2O ice
251        if (h2o_ice(i,islope) > 0._dp) then
252            albedo4PCM(i,islope) = albedo_h2oice
253            emissivity4PCM(i,islope) = 1._dp
254        end if
255
256        ! CO2 ice (dominant over H2O ice)
257        if (co2_ice(i,islope) > 0._dp) then
258            albedo4PCM(i,islope) = albedo_co2ice(icap)
259            emissivity4PCM(i,islope) = albedo_co2ice(icap)
260        end if
261
262        ! H2O frost (subject to threshold)
263        if (h2o_frost4PCM(i,islope) > h2ofrost_albedo_threshold) then
264            albedo4PCM(i,islope) = albedo_h2ofrost
265            emissivity4PCM(i,islope) = 1._dp
266        end if
267
268        ! CO2 frost (final priority)
269        if (co2_frost4PCM(i,islope) > 0.) then
270            albedo4PCM(i,islope) = albedo_co2frost(icap)
271            emissivity4PCM(i,islope) = emissivity_co2ice(icap)
272        end if
273    end do
274end do
275
276END SUBROUTINE build4PCM_surf_rad_prop
277!=======================================================================
278
279!=======================================================================
280SUBROUTINE read_albedo_emis(albedo_h2oice,albedo_h2ofrost,h2ofrost_albedo_threshold, &
281                            albedo_co2ice,albedo_co2frost,emissivity_co2ice,emissivity_baresoil)
282!-----------------------------------------------------------------------
283! NAME
284!     read_albedo_emis
285!
286! DESCRIPTION
287!     Reads albedo/emissivity parameters from "callphys.def" and
288!     "startfi.nc" ('controle').
289!
290! AUTHORS & DATE
291!     C. Metz, 01/2026
292!
293! NOTES
294!
295!-----------------------------------------------------------------------
296
297! DEPENDENCIES
298! ------------
299#ifdef CPP_IOIPSL
300use IOIPSL,          only: getin
301#else
302use ioipsl_getincom, only: getin
303#endif
304use config,          only: callphys_name
305use io_netcdf,       only: open_nc, close_nc, startfi_name, get_dim_nc, get_var_nc
306use stoppage,        only: stop_clean
307
308! DECLARATION
309! -----------
310implicit none
311
312! ARGUMENTS
313! ---------
314real(dp),               intent(out) :: albedo_h2oice, albedo_h2ofrost, h2ofrost_albedo_threshold, emissivity_baresoil
315real(dp), dimension(2), intent(out) :: albedo_co2frost, emissivity_co2ice, albedo_co2ice
316
317! LOCAL VARIABLES
318! ---------------
319logical(k4)                         :: here, cst_h2oice_albedo
320integer(di)                         :: i, nindex
321real(dp), dimension(:), allocatable :: controle
322
323! CODE
324! ----
325inquire(file = callphys_name,exist = here)
326if (.not. here) call stop_clean(__FILE__,__LINE__,'cannot find "'//callphys_name//'"!', 1)
327
328! Read albedos of H2O frost and perennial ice from "callphys.def"
329albedo_h2oice = 0.35_dp ! Default
330call getin('albedo_h2o_cap',albedo_h2oice)
331if (albedo_h2oice < 0._dp .or. albedo_h2oice > 1._dp) call stop_clean(__FILE__,__LINE__,'''albedo_h2o_cap'' in "'//callphys_name//'" is out of bounds [0,1]!',1)
332
333albedo_h2ofrost = 0.35_dp ! Default
334call getin('albedo_h2ofrost',albedo_h2ofrost)
335if (albedo_h2ofrost < 0._dp .or. albedo_h2ofrost > 1._dp) call stop_clean(__FILE__,__LINE__,'''albedo_h2ofrost'' in "'//callphys_name//'" is out of bounds [0,1]!',1)
336
337h2ofrost_albedo_threshold = 0.005_dp ! Default [kg/m2]
338call getin('frost_albedo_threshold',h2ofrost_albedo_threshold)
339
340cst_h2oice_albedo = .false. ! Default
341call getin('cst_cap_albedo',cst_h2oice_albedo)
342if (cst_h2oice_albedo) albedo_h2ofrost = albedo_h2oice ! If true, then we don't account for water frost albedo effect
343
344! Read albedos of CO2 perennial ice from "callphys.def"
345albedo_co2ice(1) = 0.6_dp ! Default, north hemisphere
346albedo_co2ice(2) = 0.85_dp ! Default, south hemisphere
347call getin('albedo_co2ice_north',albedo_co2ice(1))
348call getin('albedo_co2ice_south',albedo_co2ice(2))
349if (albedo_co2ice(1) < 0._dp .or. albedo_co2ice(1) > 1._dp) call stop_clean(__FILE__,__LINE__,'''albedo_co2ice_north'' in "'//callphys_name//'" is out of bounds [0,1]!',1)
350if (albedo_co2ice(2) < 0._dp .or. albedo_co2ice(2) > 1._dp) call stop_clean(__FILE__,__LINE__,'''albedo_co2ice_south'' in "'//callphys_name//'" is out of bounds [0,1]!',1)
351
352! Read albedos of CO2 frost from "callphys.def"
353call open_nc(startfi_name,'read')
354call get_dim_nc('index',nindex)
355allocate(controle(nindex))
356call get_var_nc('controle',controle)
357call close_nc(startfi_name)
358albedo_co2frost(1)   = controle(22)
359albedo_co2frost(2)   = controle(23)
360emissivity_co2ice(1) = controle(24)
361emissivity_co2ice(2) = controle(25)
362emissivity_baresoil  = controle(26)
363deallocate(controle)
364do i = 1,2
365    if (albedo_co2frost(i) < 0._dp .or. albedo_co2frost(i) > 1._dp) call stop_clean(__FILE__,__LINE__,'''albedo_co2frost'' from ''controle(22:23)'' in "'//startfi_name//'" is out of bounds [0,1]!',1)
366    if (emissivity_co2ice(i) < 0._dp .or. emissivity_co2ice(i) > 1._dp) call stop_clean(__FILE__,__LINE__,'''emissivity_co2ice'' from ''controle(24:25)'' in "'//startfi_name//'" is out of bounds [0,1]!',1)
367end do
368if (emissivity_baresoil < 0._dp .or. emissivity_baresoil > 1._dp) call stop_clean(__FILE__,__LINE__,'''emissivity_baresoil'' from ''controle(26)'' in "'//startfi_name//'" is out of bounds [0,1]!',1)
369
370END SUBROUTINE read_albedo_emis
371!=======================================================================
372
373END MODULE surface
Note: See TracBrowser for help on using the repository browser.