source: trunk/LMDZ.COMMON/libf/evolution/surf_temp.F90 @ 4074

Last change on this file since 4074 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: 9.9 KB
Line 
1MODULE surf_temp
2!-----------------------------------------------------------------------
3! NAME
4!     surf_temp
5!
6! DESCRIPTION
7!     Surface temperature management.
8!
9! AUTHORS & DATE
10!     JB Clement, 2025
11!
12! NOTES
13!
14!-----------------------------------------------------------------------
15
16! DEPENDENCIES
17! ------------
18use numerics, only: dp, di, k4, minieps
19
20! DECLARATION
21! -----------
22implicit none
23
24! PARAMATERS
25! ----------
26real(dp), dimension(:,:), allocatable, protected :: tsurf_PCM ! Surface temperature in the PCM at the beginning [K]
27
28contains
29!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
30
31!=======================================================================
32SUBROUTINE ini_surf_temp()
33!-----------------------------------------------------------------------
34! NAME
35!     ini_surf_temp
36!
37! DESCRIPTION
38!     Initialize the parameters of module 'surf_temp'.
39!
40! AUTHORS & DATE
41!     JB Clement, 12/2025
42!
43! NOTES
44!
45!-----------------------------------------------------------------------
46
47! DEPENDENCIES
48! ------------
49use geometry, only: ngrid, nslope
50
51! DECLARATION
52! -----------
53implicit none
54
55! CODE
56! ----
57if (.not. allocated(tsurf_PCM)) allocate(tsurf_PCM(ngrid,nslope))
58tsurf_PCM(:,:) = 0._dp
59
60END SUBROUTINE ini_surf_temp
61!=======================================================================
62
63!=======================================================================
64SUBROUTINE end_surf_temp()
65!-----------------------------------------------------------------------
66! NAME
67!     end_surf_temp
68!
69! DESCRIPTION
70!     Deallocate surf_temp arrays.
71!
72! AUTHORS & DATE
73!     JB Clement, 12/2025
74!
75! NOTES
76!
77!-----------------------------------------------------------------------
78
79! DECLARATION
80! -----------
81implicit none
82
83! CODE
84! ----
85if (allocated(tsurf_PCM)) deallocate(tsurf_PCM)
86
87END SUBROUTINE end_surf_temp
88!=======================================================================
89
90!=======================================================================
91SUBROUTINE set_tsurf_PCM(tsurf_PCM_in)
92!-----------------------------------------------------------------------
93! NAME
94!     set_tsurf_PCM
95!
96! DESCRIPTION
97!     Setter for 'tsurf_PCM'.
98!
99! AUTHORS & DATE
100!     JB Clement, 12/2025
101!
102! NOTES
103!
104!-----------------------------------------------------------------------
105
106! DECLARATION
107! -----------
108implicit none
109
110! ARGUMENTS
111! ---------
112real(dp), dimension(:,:), intent(in) :: tsurf_PCM_in
113
114! CODE
115! ----
116tsurf_PCM(:,:) = tsurf_PCM_in(:,:)
117
118END SUBROUTINE set_tsurf_PCM
119!=======================================================================
120
121!=======================================================================
122SUBROUTINE adapt_tsurf2disappearedice(surfice,is_ice_ini,ice_disappeared,tsurf_avg)
123!-----------------------------------------------------------------------
124! NAME
125!     adapt_tsurf2disappearedice
126!
127! DESCRIPTION
128!     Update surface temperature where CO2 ice disappeared using nearby
129!     bare soil temperature.
130!
131! AUTHORS & DATE
132!     JB Clement, 2025
133!
134! NOTES
135!
136!-----------------------------------------------------------------------
137
138! DEPENDENCIES
139! ------------
140use geometry, only: ngrid, nslope, nlon, nlat, latitudes, vect2lonlat, lonlat2vect
141use display,  only: print_msg
142use utility,  only: int2str
143
144! DECLARATION
145! -----------
146implicit none
147
148! ARGUMENTS
149! ---------
150real(dp),    dimension(:,:), intent(in)    :: surfice
151logical(k4), dimension(:,:), intent(in)    :: is_ice_ini
152real(dp),    dimension(:,:), intent(inout) :: tsurf_avg
153logical(k4), dimension(:,:), intent(inout) :: ice_disappeared
154
155! LOCAL VARIABLES
156! ---------------
157integer(di)                              :: islope, i, j, k, radius, rmax, d_i, d_j, ii, jj
158logical(k4)                              :: found
159real(dp),    dimension(nlon,nlat,nslope) :: tsurf_ll, co2ice_ll, mask_co2ice_ini, mask_co2ice_disappeared
160real(dp),    dimension(nlon,nlat)        :: latitude_ll
161real(dp),    dimension(ngrid)            :: tmp
162integer(di), dimension(nslope - 1)       :: priority
163
164! CODE
165! ----
166! Check to escape the subroutine (not relevant in 1D)
167if (ngrid == 1) return
168
169call print_msg("> Adapting surface temperature where ice disappeared")
170! Convert from reduced grid to lon-lat grid
171call vect2lonlat(nlon,nlat,ngrid,latitudes,latitude_ll)
172do islope = 1,nslope
173    call vect2lonlat(nlon,nlat,ngrid,tsurf_avg(:,islope),tsurf_ll(:,:,islope))
174    call vect2lonlat(nlon,nlat,ngrid,surfice(:,islope),co2ice_ll(:,:,islope))
175    call vect2lonlat(nlon,nlat,ngrid,merge(1._dp,0._dp,is_ice_ini(:,islope)),mask_co2ice_ini(:,:,islope))
176    call vect2lonlat(nlon,nlat,ngrid,merge(1._dp,0._dp,ice_disappeared(:,islope)),mask_co2ice_disappeared(:,:,islope))
177end do
178
179! For each point where ice disappeared
180rmax = max(nlon,nlat)
181do j = 1,nlat
182    do i = 1,nlon
183        do islope = 1,nslope
184            if (mask_co2ice_ini(i,j,islope) > 0.5_dp .and. co2ice_ll(i,j,islope) < minieps .and. mask_co2ice_disappeared(i,j,islope) < 0.5_dp) then
185                found = .false.
186                mask_co2ice_disappeared(i,j,islope) = 1._dp
187                call get_slope_priority(latitude_ll(i,j),islope,priority)
188                do k = 1,nslope - 1
189                    if (mask_co2ice_ini(i,j,priority(k)) < 0.5_dp) then
190                        tsurf_ll(i,j,islope) = tsurf_ll(i,j,priority(k))
191                        found = .true.
192                        exit
193                    end if
194                end do
195
196                radius = 1
197                do while (.not. found .and. radius <= rmax) ! Only if no adjacent slopes holds bare soil
198                    do d_j = -radius,radius
199                        do d_i = -radius,radius
200                            if (abs(d_i) + abs(d_j) == radius) then
201                                ii = i + d_i
202                                jj = j + d_j
203                                ! Longitudinal periodicity
204                                if (ii < 1) then
205                                    ii = ii + nlon
206                                else if (ii > nlon) then
207                                    ii = ii - nlon
208                                end if
209                                ! Latitude boundaries
210                                if (jj >= 1 .and. jj <= nlat) then
211                                    call get_slope_priority(latitude_ll(ii,jj),islope,priority)
212                                    do k = 1,nslope - 1
213                                        if (mask_co2ice_ini(ii,jj,priority(k)) < 0.5_dp) then
214                                            tsurf_ll(i,j,islope) = tsurf_ll(ii,jj,priority(k))
215                                            found = .true.
216                                            exit
217                                        end if
218                                    end do
219                                end if
220                            end if
221                            if (found) exit
222                        end do
223                        if (found) exit
224                    end do
225                    radius = radius + 1
226                end do
227                if (.not. found) call print_msg("Warning: no bare soil found for ice disappeared at point ("//int2str(i)//','//int2str(j)//','//int2str(islope)//'!')
228            end if
229        end do
230    end do
231end do
232
233! Convert back from lon-lat grid to reduced grid
234do islope = 1,nslope
235    call lonlat2vect(nlon,nlat,ngrid,tsurf_ll(:,:,islope),tsurf_avg(:,islope))
236    call lonlat2vect(nlon,nlat,ngrid,mask_co2ice_disappeared(:,:,islope),tmp)
237    where (tmp > 0.5_dp) ice_disappeared(:,islope) = .true.
238end do
239
240END SUBROUTINE adapt_tsurf2disappearedice
241!=======================================================================
242
243!=======================================================================
244SUBROUTINE get_slope_priority(lat,islope,priority)
245!-----------------------------------------------------------------------
246! NAME
247!     get_slope_priority
248!
249! DESCRIPTION
250!     Determine slope priority based on latitudes (equator-ward favored).
251!
252! AUTHORS & DATE
253!     JB Clement, 2025
254!
255! NOTES
256!     Equator-ward slopes are most likely to hold no ice.
257!-----------------------------------------------------------------------
258
259! DEPENDENCIES
260! ------------
261use geometry, only: nslope
262
263! DECLARATION
264! -----------
265implicit none
266
267! ARGUMENTS
268! ---------
269real(dp),                  intent(in)  :: lat      ! Latitude [degrees]
270integer(di),               intent(in)  :: islope
271integer(di), dimension(:), intent(out) :: priority ! Priority ordering of slopes
272
273! LOCAL VARIABLES
274! ---------------
275integer(di) :: i, k
276
277! CODE
278! ----
279k = 1
280
281! Northern hemisphere
282if (lat > 0.) then
283    ! Equator-ward slopes
284    do i = islope - 1,1,-1
285        priority(k) = i
286        k = k + 1
287    end do
288    ! Pole-ward slopes
289    do i = islope + 1,nslope
290        priority(k) = i
291        k = k + 1
292    end do
293else ! Southern hemisphere
294    ! Equator-ward slopes
295    do i = islope + 1,nslope
296        priority(k) = i
297        k = k + 1
298    end do
299    ! Pole-ward slopes
300    do i = islope - 1,1,-1
301        priority(k) = i
302        k = k + 1
303    end do
304end if
305
306END SUBROUTINE get_slope_priority
307!=======================================================================
308
309!=======================================================================
310SUBROUTINE build4PCM_tsurf(tsurf_avg,tsurf_dev,tsurf4PCM)
311!-----------------------------------------------------------------------
312! NAME
313!     build4PCM_tsurf
314!
315! DESCRIPTION
316!     Reconstructs surface temperature for the PCM.
317!
318! AUTHORS & DATE
319!     JB Clement, 12/2025
320!
321! NOTES
322!
323!-----------------------------------------------------------------------
324
325! DEPENDENCIES
326! ------------
327use display, only: print_msg
328
329! DECLARATION
330! -----------
331implicit none
332
333! ARGUMENTS
334! ---------
335real(dp), dimension(:,:), intent(in)  :: tsurf_avg, tsurf_dev
336real(dp), dimension(:,:), intent(out) :: tsurf4PCM
337
338! CODE
339! ----
340call print_msg('> Building surface temperature for the PCM')
341tsurf4PCM(:,:) = tsurf_avg(:,:) + tsurf_dev(:,:)
342
343END SUBROUTINE build4PCM_tsurf
344!=======================================================================
345
346END MODULE surf_temp
Note: See TracBrowser for help on using the repository browser.