Changeset 3862 for trunk/LMDZ.PLUTO

Timestamp:

Jul 23, 2025, 9:32:26 AM (12 days ago)

Author:

tbertrand

Message:

PLUTO PCM:
Newstart.F : new option (bugn2_new) to correct for warm N2 patches that form when interpolating on the grid
TB

File:

• trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F (modified) (4 diffs)

trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F

@@ -158 +158 @@
       integer :: iref,jref,iref1,jref1,iref2,jref2
       integer :: igref,igref1,igref2,igref3
+      integer :: edge(4) ! index of the adjecent points, N,S,E,W
+      integer :: ignew,compt
 
       INTEGER :: itau
@@ -672 +674 @@
         write(*,*) 'digsp: specific to dig SP'
         write(*,*) 'bugn2: correct bug of warm n2 due to HighRes'
+        write(*,*) 'bugn2_new: correct bug of warm n2 due to HighRes'
         write(*,*) 'bugch4: correct bug of warm ch4 due to HighRes'
         write(*,*) 'flatnosp : topo flat except SP (topo controled)'
@@ -2753 +2756 @@
 !        END DO
 
+c       bugn2_new correct bug warm n2 due to change of resolution
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'bugn2_new') then
+                
+          write(*,*) 'Where do you want to apply the change?'
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+ 564      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 564
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+ 565      read(*,*,iostat=ierr) val3,val4
+          if(ierr.ne.0) goto 565
+          write(*,*) 'Choice of topography range ?'
+ 566      read(*,*,iostat=ierr) val5,val6
+          if(ierr.ne.0) goto 566
+
+          ! let's find where we need to apply the correction
+          ierr=1
+          do while (ierr.ne.0)
+           ierr=0
+           do ig=1,ngridmx
+            IF (qsurf(ig,igcm_n2).gt.0.) then
+
+             IF (   (latfi(ig)*180./pi.ge.val)  .and.
+     &              (latfi(ig)*180./pi.le.val2)  .and.
+     &              (lonfi(ig)*180./pi.ge.val3)  .and.
+     &              (lonfi(ig)*180./pi.le.val4)  ) then
+
+              IF (   (phisfi(ig).ge.val5)  .and.
+     &               (phisfi(ig).le.val6)  ) then
+                 ! Get the index of the adjecent points, N,S,E,W
+                 call get_next_point(ig,edge) 
+                 ignew=0 
+                 DO compt=1,4
+                    if (qsurf(edge(compt),igcm_n2).le.1.e-6) then
+                       ignew=edge(compt)
+                    endif
+                 ENDDO
+                 if (ignew.gt.0) then  ! Applying correction
+                   ierr=1
+                   write(*,*) 'ig=',ig,' replaced with ',ignew
+                   qsurf(ig,igcm_n2)=0.
+                   tsurf(ig)=tsurf(ignew)
+                   DO l=1,nsoilmx
+                       tsoil(ig,l)=tsoil(ignew,l)
+                   ENDDO
+                 else
+                   write(*,*) 'No solution for ig=',ig
+                   write(*,*) '  edge=',edge
+                 endif
+
+              ENDIF
+             ENDIF
+            ENDIF
+           end do
+          end do
+
 c       flatnosp : flat topo outside a specific terrain (SP)
 c       -----------------------------------------------------
@@ -4124 +4183 @@
       end
 
+
+      SUBROUTINE get_next_point(ig,edge)
+
+      IMPLICIT NONE
+
+      include "dimensions.h"
+c=======================================================================
+c  Auteur :  T. Bertrand
+c   Objet:
+c   Get next point where arr=0
+c=======================================================================
+c   Arguments:
+c   ----------
+      INTEGER ig
+      INTEGER :: edge(4) ! index of the adjecent points, N,S,E,W
+
+c   Local:
+c   ------
+      integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)
+
+      ! Selection of the adjacent index depending on the grid point
+      !! poles (now well implemented yet) 
+      IF (ig.eq.1) then 
+         edge = (/2, 2+int(iim/4),2+2*int(iim/4),iim+1 /)
+      ELSEIF (ig.eq.ngridmx) then
+         edge =(/ngridmx-iim,ngridmx-int(iim/4),
+     &                        ngridmx-2*int(iim/4),ngridmx-1 /)
+      !! 9 other cases:  edges East, west, or in the center , at edges north south or in the center
+      ELSEIF (mod(ig,iim).eq.2) then ! Edge east :  N,S,W,E
+          IF (ig.eq.2) then
+             edge = (/1, ig+iim,ig-1+iim,ig+1 /)
+          ELSEIF (ig.eq.ngridmx-iim) then
+             edge = (/ig-iim, ngridmx,ig-1+iim,ig+1 /)
+          ELSE
+             edge = (/ig-iim, ig+iim,ig-1+iim,ig+1 /)
+          ENDIF
+      ELSEIF (mod(ig,iim).eq.1) then ! Edge west
+          IF (ig.eq.iim+1) then
+             edge = (/1,ig+iim,ig-1,ig+1-iim /)
+          ELSEIF (ig.eq.ngridmx-1) then
+             edge = (/ig-iim,ngridmx,ig-1,ig+1-iim /)
+          ELSE
+             edge = (/ig-iim,ig+iim,ig-1,ig+1-iim /)
+          ENDIF
+      ELSE
+          IF ((ig.gt.2).and.(ig.lt.iim+1)) then
+             edge = (/1,ig+iim,ig-1,ig+1 /)
+          ELSEIF ((ig.gt.ngridmx-iim).and.(ig.lt.ngridmx-1)) then
+             edge = (/ig-iim,ngridmx,ig-1,ig+1 /)
+          ELSE
+             edge = (/ig-iim,ig+iim,ig-1,ig+1 /)
+          ENDIF
+      ENDIF
+
+      end subroutine get_next_point