Changeset 3941 for trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/newstart.F

Timestamp:

Oct 29, 2025, 9:10:13 PM (5 weeks ago)

Author:

tbertrand

Message:

PLUTO PCM :
Some options in newstart for Triton simulations
TB

File:

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

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

@@ -564 +564 @@
 
           write(*,*) "Enter value of albedo of the bare ground:"
+          write(*,*) "Might not count"
           read(*,*) alb(1,1)
           alb(:,:)=alb(1,1)
 
           write(*,*) "Enter value of thermal inertia of soil:"
+          write(*,*) "Might not count"
           read(*,*) surfith(1,1)
           surfith(:,:)=surfith(1,1)
@@ -599 +601 @@
         ! copy soil thermal inertia
         ithfi(:,:)=inertiedat(:,:)
+        ! copy topo
+        CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phisold_newgrid,phisfi)
+        ! Initialize n2frac
+        n2fracfi(:)=1.
 
         ierr= NF_CLOSE(nid)
 
       else if (choix_1.eq.1) then
-         !do nothing, start and startfi have already been read
+        !do nothing, start and startfi have already been read
+        ! Initialize n2frac
+        n2fracfi(:)=1.
       else
         CALL exit(1)
@@ -1891 +1899 @@
  533      read(*,*,iostat=ierr) val3
           if(ierr.ne.0) goto 533
-!          write(*,*) 'iip1,jjp1=',iip1,jjp1,ngridmx
-
-          DO ig=1,ngridmx
+
+          if (val1.le.0.) then ! Case asym in S. hemisphere
+
+           DO ig=1,ngridmx
              ! Latitude of the sinusoid:
              val11=val1+val2*cos(lonfi(ig)*1.57079633*2/pi-val3*pi/180.)
-             ! If above line and ice: remove ice
+             ! If we are above the sinus line: we remove any ice
              IF ( (latfi(ig)*180./pi.ge.val11) .and.
      &             (latfi(ig)*180./pi.le.val1+val2) .and.
      &             (qsurf(ig,igcm_n2).gt.0.) ) then
-               ! Look for same longitude point where no ice
+               ! Looking for point at same longitude but northward with no ice
+               ! to apply its surface and soil temperatures
                val5=1.
                jref=ig
@@ -1906 +1916 @@
                ! ... iip1 ... x (jjp1-2)     32 x 23
                ! 1
-               do while (val5.ge.1..and.jref.gt.iip1)
-                  ! northward point
+               do while (val5.ge.1..and.jref.gt.iip1) ! will exit if qsurf=0 
+                  ! We take the northward point
                   jref=jref-iip1+1
-                  ! ice at the point
+                  ! We check of ice is present at that point. If not we exit
                   val5=qsurf(jref,igcm_n2)
-!                 write(*,*) jref,
-!     &                   latfi(jref)*180./pi,lonfi(jref)*180/pi,val5
                enddo
-               if (val5.ge.1) write(*,*) 'NO POINT FOUND WITH NO ICE'
+               if (val5.ge.1) write(*,*) 'PB NO POINT FOUND WITH NO ICE'
 
                ! Copy point in the selected area
@@ -1920 +1928 @@
                qsurf(ig,igcm_n2)=qsurf(jref,igcm_n2)
                qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice)
-               qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice)
+               !qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice)
                DO l=1,nsoilmx
                   tsoil(ig,l)=tsoil(jref,l)
                ENDDO
              ENDIF
-             ! If below line and no ice: add ice
+             ! If we are below the sinus line and there no ice: we add ice
              IF ( (latfi(ig)*180./pi.le.val11) .and.
-     &             (latfi(ig)*180./pi.le.val1+val2) .and.
+     &             (latfi(ig)*180./pi.ge.val1-val2) .and.
      &             (qsurf(ig,igcm_n2).eq.0.) ) then
-               ! Look for same longitude point where ice
                val5=1.
                jref=ig
@@ -1935 +1942 @@
                   ! southward point
                   jref=jref+iip1-1
-                  ! ice at the point
+                  ! We check of ice is present at that point. If yes we exit
                   val5=qsurf(jref,igcm_n2)
-                  write(*,*) jref,
-     &                   latfi(jref)*180./pi,lonfi(jref)*180/pi,val5
                enddo
-               if (val5.le.1) write(*,*) 'NO POINT FOUND WITH ICE'
+               if (val5.le.1) write(*,*) 'PB NO POINT FOUND WITH ICE'
 
                ! Copy point in the selected area
@@ -1946 +1951 @@
                qsurf(ig,igcm_n2)=qsurf(jref,igcm_n2)
                qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice)
-               qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice)
+               !qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice)
                DO l=1,nsoilmx
                   tsoil(ig,l)=tsoil(jref,l)
@@ -1952 +1957 @@
              ENDIF
 
-          ENDDO
+           ENDDO
+          
+          else ! Case N. hemisphere
+
+           DO ig=1,ngridmx
+             ! Latitude of the sinusoid:
+             val11=val1+val2*cos(lonfi(ig)*1.57079633*2/pi-val3*pi/180.)
+             ! If we are below the sinus line: we remove any ice
+             IF ( (latfi(ig)*180./pi.le.val11) .and.
+     &             (latfi(ig)*180./pi.ge.val1-val2) .and.
+     &             (qsurf(ig,igcm_n2).gt.0.) ) then
+               ! Looking for point at same longitude but southward with no ice
+               ! to apply its surface and soil temperatures
+               val5=1.
+               jref=ig
+               ! 1
+               ! ... iip1 ... x (jjp1-2)     32 x 23
+               ! 1
+               do while (val5.ge.1..and.jref.lt.ngridmx-iip1) ! will exit if qsurf=0 
+                  ! We take the southward point
+                  jref=jref+iip1-1
+                  ! We check of ice is present at that point. If not we exit
+                  val5=qsurf(jref,igcm_n2)
+               enddo
+               if (val5.ge.1) write(*,*) 'PB NO POINT FOUND WITH NO ICE'
+
+               ! Copy point in the selected area
+               tsurf(ig)=tsurf(jref)
+               qsurf(ig,igcm_n2)=qsurf(jref,igcm_n2)
+               qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice)
+               !qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice)
+               DO l=1,nsoilmx
+                  tsoil(ig,l)=tsoil(jref,l)
+               ENDDO
+             ENDIF
+             ! If we are above the sinus line and there no ice: we add ice
+             IF ( (latfi(ig)*180./pi.ge.val11) .and.
+     &             (latfi(ig)*180./pi.le.val1+val2) .and.
+     &             (qsurf(ig,igcm_n2).eq.0.) ) then
+               val5=1.
+               jref=ig
+               do while (val5.le.1.and.jref.gt.iip1)
+                  ! northward point
+                  jref=jref-iip1-1
+                  ! We check of ice is present at that point. If yes we exit
+                  val5=qsurf(jref,igcm_n2)
+               enddo
+               if (val5.le.1) write(*,*) 'PB NO POINT FOUND WITH ICE'
+
+               ! Copy point in the selected area
+               tsurf(ig)=tsurf(jref)
+               qsurf(ig,igcm_n2)=qsurf(jref,igcm_n2)
+               qsurf(ig,igcm_ch4_ice)=qsurf(jref,igcm_ch4_ice)
+               !qsurf(ig,igcm_co_ice)=qsurf(jref,igcm_co_ice)
+               DO l=1,nsoilmx
+                  tsoil(ig,l)=tsoil(jref,l)
+               ENDDO
+             ENDIF
+
+           ENDDO
+          
+          endif ! Case N. hemisphere
 
 c       source_ch4 : adding source ch4