Changeset 3979

Timestamp:

Nov 28, 2025, 5:34:32 PM (2 days ago)

Author:

jbclement

Message:

PEM:
Addition of the periodicity to search along the longitudes to find the nearest bare soil from the place where ice disappeared + Searching with slope priority according to equator-ward orientation to try gaining efficiency + Warning if the search is unsuccessful.
JBC

Location:

trunk/LMDZ.COMMON/libf/evolution

Files:

• changelog.txt (modified) (1 diff)
• pem.F90 (modified) (1 diff)
• surf_temp.F90 (modified) (6 diffs)

trunk/LMDZ.COMMON/libf/evolution/changelog.txt

@@ -802 +802 @@
 - Deletion of the reshaping tool "reshape_XIOS_output" used to convert XIOS outputs onto the PCM grid. Thus, the PEM is now able to read directly the format of XIOS outputs.
 - Addition of subroutines to convert data between a lon x lat array and a vector.
+
+== 28/11/2025 == JBC
+Addition of the periodicity to search along the longitudes to find the nearest bare soil from the place where ice disappeared + Searching with slope priority according to equator-ward orientation to try gaining efficiency + Warning if the search is unsuccessful.

trunk/LMDZ.COMMON/libf/evolution/pem.F90

@@ -743 +743 @@
 !------------------------
 ! II_d.1 Update Tsurf
-    call update_tsurf_nearest_baresoil(ngrid,nslope,iim,jjm_value,tsurf_avg,co2_ice,is_co2ice_ini,co2ice_disappeared)
+    call update_tsurf_nearest_baresoil(ngrid,nslope,iim,jjm - 1,latitude,tsurf_avg,co2_ice,is_co2ice_ini,co2ice_disappeared)
 
     if (soil_pem) then

trunk/LMDZ.COMMON/libf/evolution/surf_temp.F90

@@ -7 +7 @@
 !=======================================================================
 
-SUBROUTINE update_tsurf_nearest_baresoil(ngrid,nslope,iim_input,jjm_input,tsurf_avg,co2_ice,is_co2ice_ini,co2ice_disappeared)
+SUBROUTINE update_tsurf_nearest_baresoil(ngrid,nslope,nlon,nlat,latitude,tsurf_avg,co2_ice,is_co2ice_ini,co2ice_disappeared)
+
+use grid_conversion, only: vect2lonlat, lonlat2vect
 
 implicit none
 
 ! Inputs:
-integer,                          intent(in) :: iim_input, jjm_input, nslope, ngrid
+integer,                          intent(in) :: nlon, nlat, nslope, ngrid
 real,    dimension(ngrid,nslope), intent(in) :: co2_ice
+real,    dimension(ngrid),        intent(in) :: latitude
 logical, dimension(ngrid,nslope), intent(in) :: is_co2ice_ini
 ! Outputs:
@@ -19 +22 @@
 logical, dimension(ngrid,nslope), intent(inout) :: co2ice_disappeared
 ! Local variables:
-real, parameter                                     :: eps = 1.e-10
-integer                                             :: islope, i, j, k, radius, rmax, di, dj, ii, jj
-logical                                             :: found
-real, dimension(iim_input + 1,jjm_input + 1,nslope) :: tsurf_ll, co2ice_ll, mask_co2ice_ini, co2ice_disappeared_ll
-real, dimension(ngrid)                              :: tmp
+real, parameter                   :: eps = 1.e-10
+integer                           :: islope, i, j, k, radius, rmax, di, dj, ii, jj
+logical                           :: found
+real, dimension(nlon,nlat,nslope) :: tsurf_ll, co2ice_ll, mask_co2ice_ini, co2ice_disappeared_ll
+real, dimension(nlon,nlat)        :: latitude_ll
+real, dimension(ngrid)            :: tmp
+integer, dimension(nslope - 1)    :: priority
 
-! Check to escape the subroutine
+! Check to escape the subroutine (not relevant in 1D)
 if (ngrid == 1) return
 
 write(*,*) "> Updating surface temperature where ice disappeared"
 ! Convert from reduced grid to lon-lat grid
-#ifndef CPP_1D
+call vect2lonlat(nlon,nlat,ngrid,latitude,latitude_ll)
 do islope = 1,nslope
-    call gr_fi_dyn(1,ngrid,iim_input + 1,jjm_input + 1,tsurf_avg(:,islope),tsurf_ll(:,:,islope))
-    call gr_fi_dyn(1,ngrid,iim_input + 1,jjm_input + 1,co2_ice(:,islope),co2ice_ll(:,:,islope))
-    call gr_fi_dyn(1,ngrid,iim_input + 1,jjm_input + 1,merge(1.,0.,is_co2ice_ini(:,islope)),mask_co2ice_ini(:,:,islope))
-    call gr_fi_dyn(1,ngrid,iim_input + 1,jjm_input + 1,merge(1.,0.,co2ice_disappeared(:,islope)),co2ice_disappeared_ll(:,:,islope))
+    call vect2lonlat(nlon,nlat,ngrid,tsurf_avg(:,islope),tsurf_ll(:,:,islope))
+    call vect2lonlat(nlon,nlat,ngrid,co2_ice(:,islope),co2ice_ll(:,:,islope))
+    call vect2lonlat(nlon,nlat,ngrid,merge(1.,0.,is_co2ice_ini(:,islope)),mask_co2ice_ini(:,:,islope))
+    call vect2lonlat(nlon,nlat,ngrid,merge(1.,0.,co2ice_disappeared(:,islope)),co2ice_disappeared_ll(:,:,islope))
 enddo
-#endif
 
 ! For each point where ice disappeared
-rmax = max(iim_input + 1,jjm_input + 1)
-do j = 1,jjm_input + 1
-    do i = 1,iim_input + 1
+rmax = max(nlon,nlat)
+do j = 1,nlat
+    do i = 1,nlon
         do islope = 1,nslope
             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
                 found = .false.
                 co2ice_disappeared_ll(i,j,islope) = 1.
-                do k = 1,nslope
-                    if (k /= islope .and. mask_co2ice_ini(i,j,k) < 0.5) then
-                        tsurf_ll(i,j,islope) = tsurf_ll(i,j,k)
+                call get_slope_priority(latitude_ll(i,j),nslope,islope,priority)
+                do k = 1,nslope - 1
+                    if (mask_co2ice_ini(i,j,priority(k)) < 0.5) then
+                        tsurf_ll(i,j,islope) = tsurf_ll(i,j,priority(k))
                         found = .true.
                         exit
@@ -56 +61 @@
 
                 radius = 1
-                do while (.not. found .and. radius <= rmax) ! only if no adjacent slopes holds bare soil
+                do while (.not. found .and. radius <= rmax) ! Only if no adjacent slopes holds bare soil
                     do dj = -radius,radius
                         do di = -radius,radius
@@ -62 +67 @@
                                 ii = i + di
                                 jj = j + dj
-                                if (ii >= 1 .and. ii <= iim_input + 1 .and. jj >= 1 .and. jj <= jjm_input + 1) then
-                                    do k = 1,nslope
-                                        if (mask_co2ice_ini(ii,jj,k) < 0.5) then
-                                            tsurf_ll(i,j,islope) = tsurf_ll(i,j,k)
+                                ! Longitudinal periodicity
+                                if (ii < 1) then
+                                    ii = ii + nlon
+                                else if (ii > nlon) then
+                                    ii = ii - nlon
+                                endif
+                                ! Latitude boundaries
+                                if (jj >= 1 .and. jj <= nlat) then
+                                    call get_slope_priority(latitude_ll(ii,jj),nslope,islope,priority)
+                                    do k = 1,nslope - 1
+                                        if (mask_co2ice_ini(ii,jj,priority(k)) < 0.5) then
+                                            tsurf_ll(i,j,islope) = tsurf_ll(ii,jj,priority(k))
                                             found = .true.
                                             exit
@@ -78 +91 @@
                     radius = radius + 1
                 enddo
+                if (.not. found) write(*,*) "WARNING: no bare soil found for ice disappeared on  point:",i,j,islope
             endif
-
         enddo
     enddo
@@ -85 +98 @@
 
 ! Convert back from lon-lat grid to reduced grid
-#ifndef CPP_1D
 do islope = 1,nslope
-    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,tsurf_ll(:,:,islope),tsurf_avg(:,islope))
-    call gr_dyn_fi(1,iim_input + 1,jjm_input + 1,ngrid,co2ice_disappeared_ll(:,:,islope),tmp)
+    call lonlat2vect(nlon,nlat,ngrid,tsurf_ll(:,:,islope),tsurf_avg(:,islope))
+    call lonlat2vect(nlon,nlat,ngrid,co2ice_disappeared_ll(:,:,islope),tmp)
     where (tmp > 0.5) co2ice_disappeared(:,islope) = .true.
 enddo
-#endif
 
 END SUBROUTINE update_tsurf_nearest_baresoil
 
+!=======================================================================
+SUBROUTINE get_slope_priority(lat,nslope,islope,priority)
+! Priority given to equator-ward slope which are most likely to hold no ice
+
+implicit none
+
+! Inputs:
+real,    intent(in) :: lat
+integer, intent(in) :: nslope, islope
+! Outputs:
+integer, dimension(nslope - 1), intent(out) :: priority
+! Locals:
+integer :: i, k
+
+! Code
+!-----
+k = 1
+
+! Northern hemisphere
+if (lat > 0.) then
+    ! Equator-ward slopes
+    do i = islope - 1,1,-1
+        priority(k) = i
+        k = k + 1
+    enddo
+    ! Pole-ward slopes
+    do i = islope + 1,nslope
+        priority(k) = i
+        k = k + 1
+    enddo
+else ! Southern hemisphere
+    ! Equator-ward slopes
+    do i = islope + 1,nslope
+        priority(k) = i
+        k = k + 1
+    enddo
+    ! Pole-ward slopes
+    do i = islope - 1,1,-1
+        priority(k) = i
+        k = k + 1
+    enddo
+endif
+
+END SUBROUTINE get_slope_priority
+
 END MODULE surf_temp