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

Last change on this file since 3991 was 3991, checked in by jbclement, 6 days ago

PEM:
Apply documentation template everywhere: standardized headers format with short description, separators between functions/subroutines, normalized code sections, aligned dependencies/arguments/variables declaration.
JBC
File size: 7.1 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! DECLARATION
17! -----------
18implicit none
19
20contains
21!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
22
23!=======================================================================
24SUBROUTINE update_tsurf_nearest_baresoil(ngrid,nslope,nlon,nlat,latitude,tsurf_avg,co2_ice,is_co2ice_ini,co2ice_disappeared)
25!-----------------------------------------------------------------------
26! NAME
27!     update_tsurf_nearest_baresoil
28!
29! DESCRIPTION
30!     Update surface temperature where CO2 ice disappeared using nearby
31!     bare soil temperature.
32!
33! AUTHORS & DATE
34!     JB Clement, 2025
35!
36! NOTES
37!
38!-----------------------------------------------------------------------
39
40! DEPENDENCIES
41! ------------
42use grid_conversion, only: vect2lonlat, lonlat2vect
43
44! DECLARATION
45! -----------
46implicit none
47
48! ARGUMENTS
49! ---------
50integer,                          intent(in)    :: nlon, nlat, nslope, ngrid ! Grid dimensions
51real,    dimension(ngrid,nslope), intent(in)    :: co2_ice                   ! CO2 ice density
52real,    dimension(ngrid),        intent(in)    :: latitude                  ! Latitude
53logical, dimension(ngrid,nslope), intent(in)    :: is_co2ice_ini             ! Initial CO2 ice flag
54real,    dimension(ngrid,nslope), intent(inout) :: tsurf_avg                 ! Average surface temperature
55logical, dimension(ngrid,nslope), intent(inout) :: co2ice_disappeared        ! Ice disappeared flag
56
57! LOCAL VARIABLES
58! ---------------
59real, parameter                      :: eps = 1.e-10
60integer                              :: islope, i, j, k, radius, rmax, di, dj, ii, jj
61logical                              :: found
62real,    dimension(nlon,nlat,nslope) :: tsurf_ll, co2ice_ll, mask_co2ice_ini, co2ice_disappeared_ll
63real,    dimension(nlon,nlat)        :: latitude_ll
64real,    dimension(ngrid)            :: tmp
65integer, dimension(nslope - 1)       :: priority
66
67! CODE
68! ----
69! Check to escape the subroutine (not relevant in 1D)
70if (ngrid == 1) return
71
72write(*,*) "> Updating surface temperature where ice disappeared"
73! Convert from reduced grid to lon-lat grid
74call vect2lonlat(nlon,nlat,ngrid,latitude,latitude_ll)
75do islope = 1,nslope
76    call vect2lonlat(nlon,nlat,ngrid,tsurf_avg(:,islope),tsurf_ll(:,:,islope))
77    call vect2lonlat(nlon,nlat,ngrid,co2_ice(:,islope),co2ice_ll(:,:,islope))
78    call vect2lonlat(nlon,nlat,ngrid,merge(1.,0.,is_co2ice_ini(:,islope)),mask_co2ice_ini(:,:,islope))
79    call vect2lonlat(nlon,nlat,ngrid,merge(1.,0.,co2ice_disappeared(:,islope)),co2ice_disappeared_ll(:,:,islope))
80enddo
81
82! For each point where ice disappeared
83rmax = max(nlon,nlat)
84do j = 1,nlat
85    do i = 1,nlon
86        do islope = 1,nslope
87            if (mask_co2ice_ini(i,j,islope) > 0.5 .and. co2ice_ll(i,j,islope) < eps .and. co2ice_disappeared_ll(i,j,islope) < 0.5) then
88                found = .false.
89                co2ice_disappeared_ll(i,j,islope) = 1.
90                call get_slope_priority(latitude_ll(i,j),nslope,islope,priority)
91                do k = 1,nslope - 1
92                    if (mask_co2ice_ini(i,j,priority(k)) < 0.5) then
93                        tsurf_ll(i,j,islope) = tsurf_ll(i,j,priority(k))
94                        found = .true.
95                        exit
96                    endif
97                enddo
98
99                radius = 1
100                do while (.not. found .and. radius <= rmax) ! Only if no adjacent slopes holds bare soil
101                    do dj = -radius,radius
102                        do di = -radius,radius
103                            if (abs(di) + abs(dj) == radius) then
104                                ii = i + di
105                                jj = j + dj
106                                ! Longitudinal periodicity
107                                if (ii < 1) then
108                                    ii = ii + nlon
109                                else if (ii > nlon) then
110                                    ii = ii - nlon
111                                endif
112                                ! Latitude boundaries
113                                if (jj >= 1 .and. jj <= nlat) then
114                                    call get_slope_priority(latitude_ll(ii,jj),nslope,islope,priority)
115                                    do k = 1,nslope - 1
116                                        if (mask_co2ice_ini(ii,jj,priority(k)) < 0.5) then
117                                            tsurf_ll(i,j,islope) = tsurf_ll(ii,jj,priority(k))
118                                            found = .true.
119                                            exit
120                                        endif
121                                    enddo
122                                endif
123                            endif
124                            if (found) exit
125                        enddo
126                        if (found) exit
127                    enddo
128                    radius = radius + 1
129                enddo
130                if (.not. found) write(*,*) "WARNING: no bare soil found for ice disappeared on  point:",i,j,islope
131            endif
132        enddo
133    enddo
134enddo
135
136! Convert back from lon-lat grid to reduced grid
137do islope = 1,nslope
138    call lonlat2vect(nlon,nlat,ngrid,tsurf_ll(:,:,islope),tsurf_avg(:,islope))
139    call lonlat2vect(nlon,nlat,ngrid,co2ice_disappeared_ll(:,:,islope),tmp)
140    where (tmp > 0.5) co2ice_disappeared(:,islope) = .true.
141enddo
142
143END SUBROUTINE update_tsurf_nearest_baresoil
144!=======================================================================
145
146!=======================================================================
147SUBROUTINE get_slope_priority(lat,nslope,islope,priority)
148!-----------------------------------------------------------------------
149! NAME
150!     get_slope_priority
151!
152! DESCRIPTION
153!     Determine slope priority based on latitude (equator-ward favored).
154!
155! AUTHORS & DATE
156!     JB Clement, 2025
157!
158! NOTES
159!     Equator-ward slopes are most likely to hold no ice.
160!-----------------------------------------------------------------------
161
162! DECLARATION
163! -----------
164implicit none
165
166! ARGUMENTS
167! ---------
168real,                           intent(in)  :: lat      ! Latitude [degrees]
169integer,                        intent(in)  :: nslope, islope
170integer, dimension(nslope - 1), intent(out) :: priority ! Priority ordering of slopes
171
172! LOCAL VARIABLES
173! ---------------
174integer :: i, k
175
176! CODE
177! ----
178k = 1
179
180! Northern hemisphere
181if (lat > 0.) then
182    ! Equator-ward slopes
183    do i = islope - 1,1,-1
184        priority(k) = i
185        k = k + 1
186    enddo
187    ! Pole-ward slopes
188    do i = islope + 1,nslope
189        priority(k) = i
190        k = k + 1
191    enddo
192else ! Southern hemisphere
193    ! Equator-ward slopes
194    do i = islope + 1,nslope
195        priority(k) = i
196        k = k + 1
197    enddo
198    ! Pole-ward slopes
199    do i = islope - 1,1,-1
200        priority(k) = i
201        k = k + 1
202    enddo
203endif
204
205END SUBROUTINE get_slope_priority
206!=======================================================================
207
208END MODULE surf_temp
Note: See TracBrowser for help on using the repository browser.