Changeset 3910 for trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto

Timestamp:

Sep 3, 2025, 11:19:36 PM (5 months ago)

Author:

tbertrand

Message:

Pluto PCM:

• adding option to use N2 ice fractional maps (n2frac) read in the startfi.nc
• adding option in newstart to correct (tsurf and tsoil) for too warm or too cold N2-free (correct_t_non2) or N2-rich (correct_t_n2) patches

TB

File:

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

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

@@ -112 +112 @@
       REAL field_input(ngridmx)
       REAL,ALLOCATABLE :: field_inputs(:,:)
+      REAL, ALLOCATABLE :: n2fracfi(:)
+      REAL, ALLOCATABLE :: n2fracdat(:,:)
 
       INTEGER i,j,l,isoil,ig,idum,k
@@ -155 +157 @@
       real val, val1, val2, val3, val4, dist, ref ! to store temporary variables
       real val5, val6, val7, val8, val9, val10,val11, val12 ! to store temporary variables
+      real :: val_lat1, val_lat2, val_lon1, val_lon2
+      real :: phi_min, phi_max, albomin, albomax !temp variables
+      real :: fact, addch4
       real :: iceith=2000 ! thermal inertia of subterranean ice
       integer :: iref,jref,iref1,jref1,iref2,jref2
@@ -233 +238 @@
       allocate(qsurf(ngridmx,nqtot))
       allocate(qsurf_tmp(ngridmx,nqtot))
+!allocate n2frac
+      allocate(n2fracfi(ngridmx))
+      allocate(n2fracdat(iip1, jjp1))
 
 ! get value of nsoilmx and allocate corresponding arrays
@@ -671 +679 @@
         write(*,*) 'bladed: add ch4 on bladed terrains'
         write(*,*) 'cpbladed: add extension bladed terrains'
+        write(*,*) 'modch4_topo: add/remove ch4 (alb+topo)' 
         write(*,*) 'slope: add slope over all long. (for triton)'
         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(*,*) 'correct_t_non2: correct surf temp of no-n2 patches'
+        write(*,*) 'correct_t_n2: correct surface temp of warm n2'
         write(*,*) 'bugch4: correct bug of warm ch4 due to HighRes'
         write(*,*) 'flatnosp : topo flat except SP (topo controled)'
@@ -689 +699 @@
         write(*,*) 'albedomap: read in an albedomap albedo.nc'
         write(*,*) 'mod_distrib_ice: modify ice distrib. from albedo'
-
+        write(*,*) 'use_n2frac: read a n2frac map'
         write(*,*)
         write(*,*) 'Please note that emis and albedo are set in physiq'
@@ -1315 +1325 @@
           ENDDO
 
+c       modch4_topo : adding/removing ch4 ice with albedo + altitude filter
+c       ----------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'modch4_topo') then
+
+          write(*,*) 'Adding or removing ch4 (alb+top)'
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+ 1200      read(*,*,iostat=ierr) val_lat1,val_lat2
+          if(ierr.ne.0) goto 1200
+
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+ 1201      read(*,*,iostat=ierr) val_lon1,val_lon2
+          if(ierr.ne.0) goto 1201
+
+          write(*,*) 'Choice topo window (phi_min phi_max)'
+ 1202      read(*,*,iostat=ierr) phi_min,phi_max
+          if(ierr.ne.0) goto 1202
+
+          write(*,*) 'Choice albedo window (albomin albomax)'
+ 1203      read(*,*,iostat=ierr) albomin,albomax
+          if(ierr.ne.0) goto 1203
+
+          write(*,*) 'Choice multiplicative factor for ch4'
+ 1204      read(*,*,iostat=ierr) fact
+          if(ierr.ne.0) goto 1204
+
+          write(*,*) 'Choice additional ch4'
+ 1205     read(*,*,iostat=ierr) addch4
+          if(ierr.ne.0) goto 1205
+
+          ! Loop over gridpoints
+          DO ig=1,ngridmx
+             IF (   (latfi(ig)*180./pi.ge.val_lat1) .and.
+     &              (latfi(ig)*180./pi.le.val_lat2) .and.
+     &              (lonfi(ig)*180./pi.ge.val_lon1) .and.
+     &              (lonfi(ig)*180./pi.le.val_lon2) .and.
+     &              (zmea(ig).ge.phi_min)         .and.
+     &              (zmea(ig).le.phi_max)         .and.
+     &              (albedodat(ig).ge.albomin)      .and.
+     &              (albedodat(ig).le.albomax) ) then
+
+                qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)*fact
+                qsurf(ig,igcm_ch4_ice)=qsurf(ig,igcm_ch4_ice)+addch4
+
+             ENDIF
+          ENDDO
+
 c       cpbladed : add extension bladed terrains
 c       ----------------------------------------------------------
@@ -2756 +2812 @@
 !        END DO
 
-c       bugn2_new correct bug warm n2 due to change of resolution
+c       correct_t_non2 correct tsurf for N2-free surfaces
 c       -----------------------------------------------------
-        else if (modif(1:len_trim(modif)) .eq. 'bugn2_new') then
+        else if (modif(1:len_trim(modif)) .eq. 'correct_t_non2') then
                 
           write(*,*) 'Where do you want to apply the change?'
@@ -2770 +2826 @@
  566      read(*,*,iostat=ierr) val5,val6
           if(ierr.ne.0) goto 566
+          write(*,*) 'Choice of temperature threshold ?'
+ 567      read(*,*,iostat=ierr) val7
+          if(ierr.ne.0) goto 567
 
           ! let's find where we need to apply the correction
           ierr=1
-          do while (ierr.ne.0)
+          do while (ierr.ne.0) ! do a loop until no change is made 
            ierr=0
            do ig=1,ngridmx
-            IF (qsurf(ig,igcm_n2).gt.0.) then
+            IF (qsurf(ig,igcm_n2).lt.1.e-6.and.tsurf(ig).le.val7) then
 
              IF (   (latfi(ig)*180./pi.ge.val)  .and.
@@ -2789 +2848 @@
                  ignew=0 
                  DO compt=1,4
-                    if (qsurf(edge(compt),igcm_n2).le.1.e-6) then
+                    if (qsurf(edge(compt),igcm_n2).le.1.e-6.and.
+     &                                       tsurf(ig).gt.val7) then
                        ignew=edge(compt)
                     endif
@@ -2797 +2857 @@
                    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       correct_t_n2 : correct tsurf of warm n2 patches
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'correct_t_n2') then
+                
+          write(*,*) 'Where do you want to apply the change?'
+          write(*,*) 'Choice band : lat1 and lat2 ?'
+ 580      read(*,*,iostat=ierr) val,val2
+          if(ierr.ne.0) goto 580
+          write(*,*) 'Choice band : lon1 and lon2 ?'
+ 581      read(*,*,iostat=ierr) val3,val4
+          if(ierr.ne.0) goto 581
+          write(*,*) 'Choice of topography range ?'
+ 582      read(*,*,iostat=ierr) val5,val6
+          if(ierr.ne.0) goto 582
+          write(*,*) 'Choice of temperature threshold ?'
+ 583      read(*,*,iostat=ierr) val7
+          if(ierr.ne.0) goto 583
+          write(*,*) 'Do you want to remove (0) or keep (1) N2 ice ?'
+ 584      read(*,*,iostat=ierr) val8
+          if(ierr.ne.0) goto 584
+
+          ! 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..and.tsurf(ig).gt.val7) 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).ge.1.e-6) then
+                       ignew=edge(compt)
+                    endif
+                 ENDDO
+                 if (ignew.gt.0) then  ! Applying correction
+                   ierr=1
+                   write(*,*) 'ig=',ig,' replaced with ',ignew
+                   if (val8.eq.0) then
+                      qsurf(ig,igcm_n2)=0.
+                   endif
                    tsurf(ig)=tsurf(ignew)
                    DO l=1,nsoilmx
@@ -3709 +3833 @@
           ENDDO
 
+c       read an n2frac map
+c       -----------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'use_n2frac') then
+
+          ! Get field 2D
+          fichnom = 'n2frac.nc'
+          ierr = NF_OPEN(fichnom, NF_NOWRITE, nid_fi_input)
+          IF (ierr .NE. NF_NOERR) THEN
+            write(6,*) 'Problem opening n2frac file:', fichnom
+            write(6,*) 'ierr = ', ierr
+            CALL ABORT
+          ENDIF
+
+          ierr = NF_INQ_VARID(nid_fi_input, "n2frac", nvarid_input)
+          IF (ierr .NE. NF_NOERR) THEN
+            PRINT*, "Could not find asked variable in n2frac.nc"
+            CALL abort
+          ENDIF
+#ifdef NC_DOUBLE
+          ierr = NF_GET_VAR_DOUBLE(nid_fi_input,nvarid_input,
+     &                                 field_input)
+#else
+          ierr = NF_GET_VAR_REAL(nid_fi_input,nvarid_input,field_input)
+#endif
+          IF (ierr .NE. NF_NOERR) THEN
+            PRINT*, "Could not get asked variable in n2frac.nc"
+            CALL abort
+          ELSE
+            PRINT*, "Got variable in n2frac.nc"
+          ENDIF
+
+          ! Transfer the data to the new variables
+          DO ig = 1, ngridmx
+            n2fracfi(ig) = field_input(ig)
+          ENDDO
+          print*, "Transferred data to n2fracfi"
+          ! Confirmation message
+          print*, "n2frac data has been successfully implemented."
+
 c       slope : add slope on all longitudes
 c       -----------------------------------------------------
@@ -4146 +4309 @@
       call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
      &                dtphys,real(day_ini),
-     &                tsurf,tsoil,ithfi,emis,albfi,q2,qsurf)
+     &                tsurf,tsoil,ithfi,emis,albfi,q2,qsurf,n2fracfi)
 !     &                cloudfrac,totalfrac,hice,
 !     &                rnat,pctsrf_sic,tslab,tsea_ice,sea_ice)
@@ -4204 +4367 @@
 
       ! Selection of the adjacent index depending on the grid point
-      !! poles (now well implemented yet) 
-      IF (ig.eq.1) then 
+      !! 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