Changeset 3228 for trunk/LMDZ.PLUTO

Timestamp:

Feb 20, 2024, 4:46:35 PM (9 months ago)

Author:

tbertrand

Message:

Pluto GCM:
Cleaning code and adding extra options for newstart.F

Location:

trunk/LMDZ.PLUTO/libf

Files:

• dynphy_lonlat/phypluto/newstart.F (modified) (3 diffs)
• phypluto/comsoil_h.F90 (modified) (1 diff)
• phypluto/condense_n2.F90 (modified) (5 diffs)
• phypluto/inifis_mod.F90 (modified) (1 diff)
• phypluto/physiq_mod.F90 (modified) (8 diffs)
• phypluto/soil.F (modified) (1 diff)
• phypluto/surfini.F (modified) (1 diff)
• phypluto/turbdiff_mod.F90 (modified) (1 diff)

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

@@ -120 +120 @@
       real mugaz ! molar mass of the atmosphere
 
-      integer ierr 
+      integer ierr,iref 
 
 c Variables on the new grid along scalar points 
@@ -585 +585 @@
       write(*,*) 'qs=x : give a uniform value to a surface tracer'
       write(*,*) 'q=profile    : specify a profile for a tracer'
+      write(*,*) 'subsoil_all : set seasonal subsurface thermal inertia'
+      write(*,*) 'diurnal_TI : set diurnal subsurface thermal inertia'
 
         write(*,*)
@@ -824 +826 @@
 
 
+
+c       subsoil_all : initialize subsurface thermal inertia
+c       --------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'subsoil_all') then
+
+          write(*,*) 'New value for subsoil thermal inertia  ?'
+ 104      read(*,*,iostat=ierr) ith_bb
+          if(ierr.ne.0) goto 104
+          write(*,*) 'thermal inertia (new value):',ith_bb
+
+          write(*,*)'At which depth (in m.) does the ice layer begin?'
+          write(*,*)'(here, the deepest soil layer extends down to:'
+     &              ,layer(1),' - ',layer(nsoilmx),')'
+          write(*,*)'write 0 for uniform value for all subsurf levels?'
+          ierr=1
+          do while (ierr.ne.0)
+           read(*,*,iostat=ierr) val2
+           write(*,*)'val2 in subsoil_all:',val2,'ierr=',ierr
+           if (ierr.eq.0) then ! got a value, but do a sanity check
+             if(val2.gt.layer(nsoilmx)) then
+               write(*,*)'Depth should be less than ',layer(nsoilmx)
+               ierr=1
+             endif
+             if(val2.lt.layer(1)) then
+              if(val2.eq.0) then
+               write(*,*)'Thermal inertia set for all subsurface layers'
+               ierr=0
+              else
+               write(*,*)'Depth should be more than ',layer(1)
+               ierr=1
+              endif
+             endif
+           endif
+          enddo ! of do while
+
+          ! find the reference index iref the depth corresponds to
+          if(val2.eq.0) then
+           iref=1
+           write(*,*)'Level selected is first level: ',layer(iref),' m'
+          else
+           do isoil=1,nsoilmx-1
+            if ((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
+     &       then
+             iref=isoil+1
+             write(*,*)'Level selected : ',layer(isoil+1),' m'
+             exit
+            endif
+           enddo
+          endif
+
+          DO ig=1,ngridmx
+             DO j=iref,nsoilmx
+                   ithfi(ig,j)=ith_bb
+             ENDDO
+          ENDDO
+
+          CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
+
+c       diurnal_TI : choice of thermal inertia values (global)
+c       ----------------------------------------------------------------
+        else if (modif(1:len_trim(modif)) .eq. 'diurnal_TI') then
+
+         write(*,*) 'New value for diurnal thermal inertia  ?'
+ 106     read(*,*,iostat=ierr) ith_bb
+         if(ierr.ne.0) goto 106
+         write(*,*) 'Diurnal thermal inertia (new value):',ith_bb
+
+         write(*,*)'At which depth (in m.) does the ice layer ends?'
+         write(*,*)'(currently, the soil layer 1 and nsoil are:'
+     &              ,layer(1),' - ',layer(nsoilmx),')'
+         ierr=1
+         do while (ierr.ne.0)
+          read(*,*,iostat=ierr) val2
+          write(*,*)'val2 in diurnal_TI:',val2,'ierr=',ierr
+          if (ierr.eq.0) then ! got a value, but do a sanity check
+            if(val2.gt.layer(nsoilmx)) then
+              write(*,*)'Depth should be less than ',layer(nsoilmx)
+              ierr=1
+            endif
+            if(val2.lt.layer(1)) then
+              write(*,*)'Depth should be more than ',layer(1)
+              ierr=1
+            endif
+          endif
+         enddo ! of do while
+
+           ! find the reference index iref the depth corresponds to
+         do isoil=1,nsoilmx-1
+            !write(*,*)'isoil, ',isoil,val2
+            !write(*,*)'lay(i),lay(i+1):',layer(isoil),layer(isoil+1),' m'
+            if ((val2.gt.layer(isoil)).and.(val2.lt.layer(isoil+1)))
+     &       then
+             iref=isoil+1
+             write(*,*)'Level selected : ',layer(isoil+1),' m'
+             exit
+            endif
+         enddo
+
+         DO ig=1,ngridmx
+             DO j=1,iref
+                   ithfi(ig,j)=ith_bb
+             ENDDO
+         ENDDO
+
+         CALL gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,ithfi,ith)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
         endif
              

trunk/LMDZ.PLUTO/libf/phypluto/comsoil_h.F90

@@ -5 +5 @@
 !integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case 
 !integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
-  integer,save :: nsoilmx = 18 ! default, but may be set in callphys.def
+  integer,save :: nsoilmx = 24 ! default, but may be set in callphys.def
 ! Full soil layer depths are set as: layer(k)=lay1_soil*alpha_soil**(k-1) , k=1,nsoil
 ! Mid soil layer depths are set as: mlayer(k)=lay1_soil*alpha_soil**(k-1/2) , k=0,nsoil-1 

trunk/LMDZ.PLUTO/libf/phypluto/condense_n2.F90

@@ -130 +130 @@
   REAL globzplevnew
 
-  REAL vmrn2(klon)
+!   REAL vmrn2(klon)
 !   SAVE vmrn2
   REAL stephan
@@ -172 +172 @@
      ! calculate global mean surface pressure for the fast mode
      IF (fast) THEN
+        IF (.not. ALLOCATED(kp)) ALLOCATE(kp(klon))
         DO ig=1,klon
            kp(ig) = exp(-phisfi(ig)/(r*38.))
@@ -178 +179 @@
      ENDIF
 
-     vmrn2(:) = 1.
-     IF (ch4lag) then
-        DO ig=1,klon
-           if (latitude(ig)*180./pi.ge.latlag) then
-              vmrn2(ig) = vmrlag
-           endif
-        ENDDO
-     ENDIF
-     IF (no_n2frost) then
-        DO ig=1,klon
-           if (picen2(ig).eq.0.) then
-              vmrn2(ig) = 1.e-15
-           endif
-        ENDDO
-     ENDIF
+     !vmrn2(:) = 1.
+     !IF (ch4lag) then
+     !   DO ig=1,klon
+     !      if (latitude(ig)*180./pi.ge.latlag) then
+     !         vmrn2(ig) = vmrlag
+     !      endif
+     !   ENDDO
+     !ENDIF
+     !IF (no_n2frost) then
+     !   DO ig=1,klon
+     !      if (picen2(ig).eq.0.) then
+     !         vmrn2(ig) = 1.e-15
+     !      endif
+     !   ENDDO
+     !ENDIF
      firstcall=.false.
   ENDIF
@@ -343 +344 @@
    DO ig=1,klon
      ! forecast of frost temperature ztcondsol
-     !ztcondsol(ig) = tcond_n2(zplev(ig),zqn2(ig,1))
-     ztcondsol(ig) = tcond_n2(zplev(ig),vmrn2(ig))
+     ztcondsol(ig) = tcond_n2(zplev(ig),zqn2(ig,1))
+     !ztcondsol(ig) = tcond_n2(zplev(ig),vmrn2(ig))
 
 !     Loop over where we have condensation / sublimation
@@ -633 +634 @@
              if((pq(ig,l,iq) +(pdqc(ig,l,iq)+ pdq(ig,l,iq))*ptimestep) &
                      .lt.0.01) then ! if n2 < 1 %  !
-               pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq)
+                     write(*,*) 'Warning: n2 < 1%'
+                     pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq)
              end if
             end if

trunk/LMDZ.PLUTO/libf/phypluto/inifis_mod.F90

@@ -176 +176 @@
      if (is_master) write(*,*) trim(rname)//": season = ",season
 
+     if (is_master) then
+       write(*,*) trim(rname)//": Fast mode (nogcm) ?"
+     endif
+     fast=.false. ! default value
+     call getin_p("fast",fast)
+     if (is_master) write(*,*) trim(rname)//": fast = ",fast
+
      if (is_master) write(*,*) trim(rname)//&
        ": No seasonal cycle: initial day to lock the run during restart"

trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90

@@ -45 +45 @@
       use comcstfi_mod, only: pi, g, rcp, r, rad, mugaz, cpp
       use time_phylmdz_mod, only: daysec
-      use callkeys_mod, only: albedo_spectral_mode, calladj, calldifv, &
+      use callkeys_mod, only: fast,albedo_spectral_mode, calladj, calldifv, &
                               callrad, callsoil, nosurf, &
                               aerohaze, corrk, diagdtau,&
@@ -293 +293 @@
       ! For Surface Tracers : (kg/m2/s)
       real dqsurf(ngrid,nq)                 ! Cumulated tendencies.
-      real zdqsurfc(ngrid)                  ! Condense_n2 routine.
+      !real zdqsurfc(ngrid)                  ! Condense_n2 routine.
       REAL zdqsc(ngrid,nq)                  ! Condense_n2 routine.
       real zdqsdif(ngrid,nq)                ! Turbdiff/vdifc routines.
@@ -768 +768 @@
      ! w = F / (rho*area) and rho = P/(r*T)
      ! But first linearly interpolate mass flux to mid-layers
-      do l=1,nlayer-1
+      if (.not.fast) then 
+       do l=1,nlayer-1
          pw(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
-      enddo
-      pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
-      do l=1,nlayer
+       enddo
+       pw(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
+       do l=1,nlayer
          pw(1:ngrid,l)=(pw(1:ngrid,l)*r*pt(1:ngrid,l)) /  &
                        (pplay(1:ngrid,l)*cell_area(1:ngrid))
-      enddo
-      ! omega in Pa/s
-      do l=1,nlayer-1
+       enddo
+       ! omega in Pa/s
+       do l=1,nlayer-1
          omega(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
-      enddo
-      omega(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
-      do l=1,nlayer
+       enddo
+       omega(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
+       do l=1,nlayer
          omega(1:ngrid,l)=g*omega(1:ngrid,l)/cell_area(1:ngrid)
-      enddo
-
+       enddo
+      endif
 !---------------------------------
 ! II. Compute radiative tendencies
@@ -922 +923 @@
                albedo_equivalent(1:ngrid)=albedo(1:ngrid,1)
                fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,1))
-               ! TB24:
                fluxabs_sw(1:ngrid)=fluxsurfabs_sw(1:ngrid)
                fluxrad_sky(1:ngrid)    = fluxsurfabs_sw(1:ngrid)
@@ -1138 +1138 @@
 
          pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer)+zdtc(1:ngrid,1:nlayer)
+         pdv(1:ngrid,1:nlayer) = pdv(1:ngrid,1:nlayer)+zdvc(1:ngrid,1:nlayer)
+         pdu(1:ngrid,1:nlayer) = pdu(1:ngrid,1:nlayer)+zduc(1:ngrid,1:nlayer)
          zdtsurf(1:ngrid)      = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid)
 
          pdq(1:ngrid,1:nlayer,1:nq)   = pdq(1:ngrid,1:nlayer,1:nq)+ zdqc(1:ngrid,1:nlayer,1:nq)
-         ! dqsurf(1:ngrid,igcm_n2_ice) = dqsurf(1:ngrid,igcm_n2_ice) + zdqsurfc(1:ngrid)
+         dqsurf(1:ngrid,igcm_n2) = dqsurf(1:ngrid,igcm_n2) + zdqsc(1:ngrid,igcm_n2)
 
 !!         call writediagfi(ngrid,"condense_n2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
@@ -1324 +1326 @@
 
       tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
-
       ! Compute soil temperatures and subsurface heat flux.
       if (callsoil) then
@@ -1378 +1379 @@
       ! Surface pressure.
       ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
-
 
 
@@ -1645 +1645 @@
       call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
       call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
-      call writediagfi(ngrid,"temp","temperature","K",3,zt)
-      call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
-      call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
-      call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
-      call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
-      call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
+
+      if (.not.fast) then
+       call writediagfi(ngrid,"temp","temperature","K",3,zt)
+       call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
+       call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
+       call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
+       call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
+       call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
+      endif
 
 !     Subsurface temperatures

trunk/LMDZ.PLUTO/libf/phypluto/soil.F

@@ -69 +69 @@
         do ik=0,nsoil-1
           mthermdiff(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa
-!         write(*,*),'soil: ik: ',ik,' mthermdiff:',mthermdiff(ig,ik)
+          !write(*,*),'soil: ik: ',ik,' mthermdiff:',mthermdiff(ig,ik)
         enddo
       enddo

trunk/LMDZ.PLUTO/libf/phypluto/surfini.F

@@ -44 +44 @@
       ! Step 2 : We get the bare ground albedo from the start files.
       DO ig=1,ngrid
-         albedo_bareground(ig)=0.1 ! TB24 albedodat(ig)
+         albedo_bareground(ig)=0.7 ! albedodat(ig)
          DO nw=1,L_NSPECTV
-            albedo(ig,nw)=0.1 !albedo_bareground(ig)
+            albedo(ig,nw)=0.7 !albedo_bareground(ig)
          ENDDO
       ENDDO

trunk/LMDZ.PLUTO/libf/phypluto/turbdiff_mod.F90

@@ -138 +138 @@
 
       IF (firstcall) THEN
-         ivap=1 !TB24
+         ivap=1 
          iliq=0
          iliq_surf=0