Changeset 3412 for trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90

Timestamp:

Aug 21, 2024, 3:46:54 PM (3 months ago)

Author:

afalco

Message:

Pluto PCM: Imported glaciers & conserv mass routines from pluto.old.
AF

File:

• trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90 (modified) (14 diffs)

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

@@ -54 +54 @@
                               fast,fasthaze,haze,metcloud,monoxcloud,&
                               n2cond,nearn2cond,noseason_day,conservn2, &
-                              convergeps,kbo,triton,paleo,paleoyears, &
+                              convergeps,kbo,triton,paleo,paleoyears,glaflow, &
                               carbox, methane,&
                               oldplutovdifc,oldplutocorrk,oldplutosedim, &
-                              aerohaze,haze_proffix,source_haze,&
+                              aerohaze,haze_proffix,source_haze, tsurfmax, &
+                              albmin_ch4, &
                               season, sedimentation,generic_condensation, &
                               specOLR, &
@@ -291 +292 @@
       REAL zplev(ngrid,nlayer+1),zplay(ngrid,nlayer)
 
+      REAL,SAVE :: ptime0    ! store the first time
+      REAL dstep
+      REAL,SAVE :: glastep=20   ! step in pluto day to spread glacier
 
 
@@ -617 +621 @@
          albedo_snow_SPECTV(:)=0.0
          albedo_n2_ice_SPECTV(:)=0.0
+
+         ptime0=ptime
+         write (*,*) 'In physiq ptime0 =', ptime
+
          call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_n2_ice_SPECTV)
 
@@ -759 +767 @@
       taux(1:ngrid) = 0.0
       tauy(1:ngrid) = 0.0
+
+      if (conservn2) then
+         write(*,*) 'conservn2 iniloop'
+         call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
+      endif
 
       zday=pday+ptime ! Compute time, in sols (and fraction thereof).
@@ -848 +861 @@
        enddo
       endif
+
+      if (conservn2) then
+         write(*,*) 'conservn2 thermo'
+         call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
+      endif
+
 !---------------------------------
 ! II. Compute radiative tendencies
@@ -1094 +1113 @@
 !!      call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_gas))
 
+   if (conservn2) then
+      write(*,*) 'conservn2 radiat'
+      call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1))
+   endif
 
 !  --------------------------------------------
@@ -1229 +1252 @@
       endif ! end of 'calldifv'
 
+      if (conservn2) then
+        write(*,*) 'conservn2 diff'
+        call testconservmass(ngrid,nlayer,pplev(:,1),qsurf(:,1)+ &
+                                       dqsurf(:,1)*ptimestep)
+      endif
 
 !-------------------
@@ -1322 +1350 @@
       endif  ! end of 'n2cond'
 
+      if (conservn2) then
+       write(*,*) 'conservn2 n2cond'
+       call testconservmass(ngrid,nlayer,pplev(:,1)+ &
+           pdpsrf(:)*ptimestep,qsurf(:,1)+dqsurf(:,1)*ptimestep)
+      endif
 
 !---------------------------------------------
@@ -1592 +1625 @@
       if (conservn2) then
          write(*,*) 'conservn2 tracer'
-         ! call testconservmass(ngrid,nlayer,pplev(:,1)+   &
-         !    pdpsrf(:)*ptimestep,qsurf(:,1))
+         call testconservmass(ngrid,nlayer,pplev(:,1)+   &
+            pdpsrf(:)*ptimestep,qsurf(:,1))
       endif
 
@@ -1646 +1679 @@
 
 
-      ! ! Increment surface temperature
-      ! if(ok_slab_ocean)then  !AF24: removed
-
+! VII.1 Increment surface temperature
+!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
-      ! Compute soil temperatures and subsurface heat flux.
+
+      ! Prevent surface (.e.g. non volatile ch4) to exceed max temperature
+      ! Lellouch et al., 2000,2011
+      IF (tsurfmax) THEN
+        DO ig=1,ngrid
+         if (albedo_equivalent(ig).gt.albmin_ch4.and. &
+                           qsurf(ig,igcm_n2).eq.0.) then
+              tsurf(ig)=min(tsurf(ig),54.)
+         endif
+        ENDDO
+      ENDIF
+
+! VII.2 Compute soil temperatures and subsurface heat flux.
+!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
       if (callsoil) then
          call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat,   &
@@ -1656 +1701 @@
       endif
 
-
-!       if (ok_slab_ocean) then !AF24: removed
+      ! ! For output :
+      ! tidat_out(:,:)=0.
+      ! DO l=1,min(nlayermx,nsoilmx)
+      !    tidat_out(:,l)=tidat(:,l)
+      ! ENDDO
 
       ! Test energy conservation
@@ -1665 +1713 @@
       endif
 
+
+
+!   VII.3 multiply tendencies of cond/subli for paleo loop only in the
+!       last Pluto year of the simulation
+!       Year day must be adapted in the startfi for each object
+!       Paleo uses year_day to calculate the annual mean tendancies
+!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      IF (paleo) then
+         if (zday.gt.day_ini+ptime0+nday-year_day) then
+            DO iq=1,nq
+             DO ig=1,ngrid
+               qsurfyear(ig,iq)=qsurfyear(ig,iq)+ &
+                             (qsurf(ig,iq)-qsurf1(ig,iq))  !kg m-2 !ptimestep
+             ENDDO
+            ENDDO
+         endif
+      endif
+
+!   VII.4 Glacial flow at each timestep glastep or at lastcall
+!   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      IF (fast.and.glaflow) THEN
+         if ((mod(zday-day_ini-ptime0,glastep)).lt.1. &
+                                                   .or.lastcall) then
+           IF (lastcall) then
+            dstep=mod(zday-day_ini-ptime0,glastep)*daysec
+           else
+            dstep=glastep*daysec
+           endif
+           zdqflow(:,:)=qsurf(:,:)
+           IF (paleo) then
+             call spreadglacier_paleo(ngrid,nq,qsurf, &
+                                    phisfinew,dstep,tsurf)
+           else
+             call spreadglacier_simple(ngrid,nq,qsurf,dstep)
+           endif
+           zdqflow(:,:)=(zdqflow(:,:)-qsurf(:,:))/dstep
+
+           if (conservn2) then
+            write(*,*) 'conservn2 glaflow'
+            call testconservmass(ngrid,nlayer,pplev(:,1)+ &
+            pdpsrf(:)*ptimestep,qsurf(:,1))
+           endif
+
+         endif
+      ENDIF
 
 !---------------------------------------------------
@@ -1816 +1909 @@
 
 
-      ! ! Test for water conservation. !AF24: removed
 
       ! Calculate RH_generic (Generic Relative Humidity) for diagnostic.
@@ -2271 +2363 @@
 
          do iq=1,nq
-          ! call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
+          call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
           !  call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
           !                   'kg m^-2',2,qsurf_hist(1,iq) )
-          !  call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
-          !                  'kg m^-2',2,qcol(1,iq) )
+           call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
+                           'kg m^-2',2,qcol(1,iq) )
 !          call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
 !                          'kg m^-2',2,qsurf(1,iq) )