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

Last change on this file since 3989 was 3989, checked in by jbclement, 5 days ago

PEM:
Massive structural refactor of the PEM codebase for improved readability, consistency and maintainability. The goal is to modernize, standardize and consolidate the code while removing legacy complexity. In detail, this change:

  • Performs large-scale cleanup removing unused code, obsolete routines, duplicated functionality and deprecated initialization logic;
  • Removes "*_mod" wrappers;
  • Replaces mixed naming conventions with clearer variable names, domain-based file/module names and purpose-based routine names;
  • Adds native reading/writing capabilities to the PEM removing the cumbersome dependency on Mars PCM subroutines;
  • Centralizes soil/ice/orbital/PEM configuration logic into dedicated modules;
  • Simplifies routines and updates calls/interfaces to match the new structure.

This refactor significantly clarifies the codebase and provides a cleaner foundation for forthcoming developments.
JBC
File size: 5.7 KB
Line 
1MODULE surf_temp
2
3implicit none
4
5!=======================================================================
6contains
7!=======================================================================
8
9!=======================================================================
10SUBROUTINE update_tsurf_nearest_baresoil(ngrid,nslope,nlon,nlat,latitude,tsurf_avg,co2_ice,is_co2ice_ini,co2ice_disappeared)
11
12use grid_conversion, only: vect2lonlat, lonlat2vect
13
14implicit none
15
16! Inputs:
17integer,                          intent(in) :: nlon, nlat, nslope, ngrid
18real,    dimension(ngrid,nslope), intent(in) :: co2_ice
19real,    dimension(ngrid),        intent(in) :: latitude
20logical, dimension(ngrid,nslope), intent(in) :: is_co2ice_ini
21! Outputs:
22real,    dimension(ngrid,nslope), intent(inout) :: tsurf_avg
23logical, dimension(ngrid,nslope), intent(inout) :: co2ice_disappeared
24! Local variables:
25real, parameter                   :: eps = 1.e-10
26integer                           :: islope, i, j, k, radius, rmax, di, dj, ii, jj
27logical                           :: found
28real, dimension(nlon,nlat,nslope) :: tsurf_ll, co2ice_ll, mask_co2ice_ini, co2ice_disappeared_ll
29real, dimension(nlon,nlat)        :: latitude_ll
30real, dimension(ngrid)            :: tmp
31integer, dimension(nslope - 1)    :: priority
32
33! Check to escape the subroutine (not relevant in 1D)
34if (ngrid == 1) return
35
36write(*,*) "> Updating surface temperature where ice disappeared"
37! Convert from reduced grid to lon-lat grid
38call vect2lonlat(nlon,nlat,ngrid,latitude,latitude_ll)
39do islope = 1,nslope
40    call vect2lonlat(nlon,nlat,ngrid,tsurf_avg(:,islope),tsurf_ll(:,:,islope))
41    call vect2lonlat(nlon,nlat,ngrid,co2_ice(:,islope),co2ice_ll(:,:,islope))
42    call vect2lonlat(nlon,nlat,ngrid,merge(1.,0.,is_co2ice_ini(:,islope)),mask_co2ice_ini(:,:,islope))
43    call vect2lonlat(nlon,nlat,ngrid,merge(1.,0.,co2ice_disappeared(:,islope)),co2ice_disappeared_ll(:,:,islope))
44enddo
45
46! For each point where ice disappeared
47rmax = max(nlon,nlat)
48do j = 1,nlat
49    do i = 1,nlon
50        do islope = 1,nslope
51            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
52                found = .false.
53                co2ice_disappeared_ll(i,j,islope) = 1.
54                call get_slope_priority(latitude_ll(i,j),nslope,islope,priority)
55                do k = 1,nslope - 1
56                    if (mask_co2ice_ini(i,j,priority(k)) < 0.5) then
57                        tsurf_ll(i,j,islope) = tsurf_ll(i,j,priority(k))
58                        found = .true.
59                        exit
60                    endif
61                enddo
62
63                radius = 1
64                do while (.not. found .and. radius <= rmax) ! Only if no adjacent slopes holds bare soil
65                    do dj = -radius,radius
66                        do di = -radius,radius
67                            if (abs(di) + abs(dj) == radius) then
68                                ii = i + di
69                                jj = j + dj
70                                ! Longitudinal periodicity
71                                if (ii < 1) then
72                                    ii = ii + nlon
73                                else if (ii > nlon) then
74                                    ii = ii - nlon
75                                endif
76                                ! Latitude boundaries
77                                if (jj >= 1 .and. jj <= nlat) then
78                                    call get_slope_priority(latitude_ll(ii,jj),nslope,islope,priority)
79                                    do k = 1,nslope - 1
80                                        if (mask_co2ice_ini(ii,jj,priority(k)) < 0.5) then
81                                            tsurf_ll(i,j,islope) = tsurf_ll(ii,jj,priority(k))
82                                            found = .true.
83                                            exit
84                                        endif
85                                    enddo
86                                endif
87                            endif
88                            if (found) exit
89                        enddo
90                        if (found) exit
91                    enddo
92                    radius = radius + 1
93                enddo
94                if (.not. found) write(*,*) "WARNING: no bare soil found for ice disappeared on  point:",i,j,islope
95            endif
96        enddo
97    enddo
98enddo
99
100! Convert back from lon-lat grid to reduced grid
101do islope = 1,nslope
102    call lonlat2vect(nlon,nlat,ngrid,tsurf_ll(:,:,islope),tsurf_avg(:,islope))
103    call lonlat2vect(nlon,nlat,ngrid,co2ice_disappeared_ll(:,:,islope),tmp)
104    where (tmp > 0.5) co2ice_disappeared(:,islope) = .true.
105enddo
106
107END SUBROUTINE update_tsurf_nearest_baresoil
108!=======================================================================
109
110!=======================================================================
111SUBROUTINE get_slope_priority(lat,nslope,islope,priority)
112! Priority given to equator-ward slope which are most likely to hold no ice
113
114implicit none
115
116! Inputs:
117real,    intent(in) :: lat
118integer, intent(in) :: nslope, islope
119! Outputs:
120integer, dimension(nslope - 1), intent(out) :: priority
121! Locals:
122integer :: i, k
123
124! Code
125!-----
126k = 1
127
128! Northern hemisphere
129if (lat > 0.) then
130    ! Equator-ward slopes
131    do i = islope - 1,1,-1
132        priority(k) = i
133        k = k + 1
134    enddo
135    ! Pole-ward slopes
136    do i = islope + 1,nslope
137        priority(k) = i
138        k = k + 1
139    enddo
140else ! Southern hemisphere
141    ! Equator-ward slopes
142    do i = islope + 1,nslope
143        priority(k) = i
144        k = k + 1
145    enddo
146    ! Pole-ward slopes
147    do i = islope - 1,1,-1
148        priority(k) = i
149        k = k + 1
150    enddo
151endif
152
153END SUBROUTINE get_slope_priority
154!=======================================================================
155
156END MODULE surf_temp
Note: See TracBrowser for help on using the repository browser.