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

Timestamp:

May 15, 2024, 8:29:23 PM (6 months ago)

Author:

afalco

Message:

Pluto PCM:
Include cooling, hazes in radiative module
AF

File:

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

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

@@ -54 +54 @@
                               fast,fasthaze,haze,metcloud,monoxcloud,&
                               n2cond,nearn2cond,noseason_day,conservn2, &
-                              kbo,triton,paleo,paleoyears, &
+                              convergeps,kbo,triton,paleo,paleoyears, &
                               carbox, methane, oldplutovdifc, oldplutocorrk, &
                               aerohaze,haze_proffix,source_haze, &
@@ -67 +67 @@
       use phys_state_var_mod
       use callcorrk_mod, only: callcorrk
-    !   use callcorrk_pluto_mod, only: callcorrk_pluto
+      use callcorrk_pluto_mod, only: callcorrk_pluto
       use vdifc_mod, only: vdifc
       use vdifc_pluto_mod, only: vdifc_pluto
@@ -260 +260 @@
       REAL,save  :: tcond1p4Pa
       DATA tcond1p4Pa/38/
+      real :: tidat(ngrid,nsoilmx)    ! thermal inertia soil
 
 
@@ -306 +307 @@
       real zlss                      ! Sub solar point longitude (radians).
       real zday                      ! Date (time since Ls=0, calculated in sols).
+      REAL,save :: saveday           ! saved date
+      REAL,save :: savedeclin        ! saved declin
       real zzlay(ngrid,nlayer)       ! Altitude at the middle of the atmospheric layers.
       real zzlev(ngrid,nlayer+1)     ! Altitude at the atmospheric layer boundaries.
@@ -619 +622 @@
          call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
 
+         savedeclin=0.
+         saveday=pday
+         !adjust=0. ! albedo adjustment for convergeps
 
 !        Initialize soil.
@@ -871 +877 @@
          if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
             mu0 = cos(pi*szangle/180.0)
+            fract= 1/(4*mu0) ! AF24: from pluto.old
          endif
 
       endif
 
-      ! AF24: TODO insert surfprop for pluto & triton around here
+
+!     Pluto albedo /IT changes depending on surface ices (only in 3D)
+!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      if (ngrid.ne.1) then
+
+         !! Specific to change albedo of N2 so that Psurf
+         !! converges toward 1.4 Pa in "1989" seasons for Triton
+         !! converges toward 1.1 Pa in "2015" seasons for Pluto
+         if (convergeps) then
+            if (triton) then
+               ! 1989 declination
+               if (declin*180./pi.gt.-46..and.declin*180./pi.lt.-45.   &
+            .and.zday.gt.saveday+1000.   &
+            .and.declin.lt.savedeclin) then
+               call globalaverage2d(ngrid,pplev(:,1),globave)
+               if (globave.gt.1.5) then
+                     adjust=adjust+0.005
+               else if (globave.lt.1.3) then
+                     adjust=adjust-0.005
+               endif
+               saveday=zday
+               endif
+            else
+               ! Pluto : 2015 declination current epoch
+               if (declin*180./pi.gt.50.and.declin*180./pi.lt.51.  &
+            .and.zday.gt.saveday+10000.  &
+            .and.declin.gt.savedeclin) then
+               call globalaverage2d(ngrid,pplev(:,1),globave)
+               if (globave.gt.1.2) then
+                     adjust=adjust+0.005
+               else if (globave.lt.1.) then
+                     adjust=adjust-0.005
+               endif
+               saveday=zday
+               endif
+            endif
+         end if
+      end if ! if ngrid ne 1
+
+      call surfprop(ngrid,nq,fract,qsurf,tsurf,tidat,   &
+      capcal,adjust,dist_star,albedo,emis,flusurfold,ptimestep,   &
+      zls)
+      ! do i=2,ngrid
+      !    albedo(i,:) = albedo(1,:)
+      ! enddo
 
       if (callrad) then
@@ -920 +971 @@
 
                ! standard callcorrk
-               ! clearsky=.false.
-            !    if (oldplutocorrk) then
-            !       call callcorrk_pluto(icount,ngrid,nlayer,pq,nq,qsurf,          &
-            !                    albedo,emis,mu0,pplev,pplay,pt,                   &
-            !                    tsurf,fract,dist_star,aerosol, &
-            !                    zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,   &
-            !                    fluxabs_sw,fluxtop_dn,reffrad,tau_col,ptime,pday, &
-            !                    firstcall,lastcall,zzlay)
-            !    else
-                call callcorrk(ngrid,nlayer,pq,nq,qsurf,                          &
+               if (oldplutocorrk) then
+                  call callcorrk_pluto(icount,ngrid,nlayer,pq,nq,qsurf,          &
+                               albedo(:,1),emis,mu0,pplev,pplay,pt,                   &
+                               zzlay,tsurf,fract,dist_star,aerosol,              &
+                               zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,   &
+                               fluxabs_sw,fluxtop_dn,reffrad,tau_col,ptime,pday, &
+                               cloudfrac,totcloudfrac,.false.,                   &
+                               firstcall,lastcall)
+               else
+                call callcorrk(ngrid,nlayer,pq,nq,qsurf,  &
                               albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,   &
                               zzlay,tsurf,fract,dist_star,aerosol,muvar,                &
@@ -938 +989 @@
                               tau_col,cloudfrac,totcloudfrac,                     &
                               .false.,firstcall,lastcall)
-            !    endif ! oldplutocorrk
+               endif ! oldplutocorrk
                 !GG (feb2021): Option to "artificially" decrease the raditive time scale in
                 !the deep atmosphere  press > 0.1 bar. Suggested by MT.
@@ -956 +1007 @@
 
                 ! AF24: removed CLFvarying Option
-
-               ! if(ok_slab_ocean) then
-               !    tsurf(:)=tsurf2(:)
-               ! endif
-
 
                ! Radiative flux from the sky absorbed by the surface (W.m-2).
@@ -1052 +1098 @@
 
          if (oldplutovdifc) then
-            zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
             zdum1(:,:) = 0
             zdum2(:,:) = 0
@@ -1066 +1111 @@
                     zdqdif,zdqsdif,qsat_ch4,qsat_ch4_l1) !,zq1temp_ch4,qsat_ch4)
 
-            pdv(1:ngrid,1:nlayer)=pdv(1:ngrid,1:nlayer)+zdvdif(1:ngrid,1:nlayer)
-            pdu(1:ngrid,1:nlayer)=pdu(1:ngrid,1:nlayer)+zdudif(1:ngrid,1:nlayer)
-            pdt(1:ngrid,1:nlayer)=pdt(1:ngrid,1:nlayer)+zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer)
-
             zdtdif(1:ngrid,1:nlayer)=zdhdif(1:ngrid,1:nlayer)*zpopsk(1:ngrid,1:nlayer) ! for diagnostic only
-            zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
 
             bcond=1./tcond1p4Pa
             acond=r/lw_n2
-
-            if (tracer) then
-               pdq(1:ngrid,1:nlayer,1:nq)=pdq(1:ngrid,1:nlayer,1:nq)+ zdqdif(1:ngrid,1:nlayer,1:nq)
-               dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
-            end if ! of if (tracer)
 
             !------------------------------------------------------------------
@@ -1379 +1414 @@
          zdqphot_prec(:,:)=0.
          zdqphot_ch4(:,:)=0.
+         zdqhaze(:,:,:)=0
          ! Forcing to a fixed haze profile if haze_proffix
          if (haze_proffix.and.i_haze.gt.0.) then
-         call haze_prof(ngrid,nlayer,zzlay,pplay,pt,  &
-                                  reffrad,profmmr)
-         zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_haze))  &
-                                                  /ptimestep
+            call haze_prof(ngrid,nlayer,zzlay,pplay,pt,  &
+                           reffrad,profmmr)
+            zdqhaze(:,:,i_haze)=(profmmr(:,:)-pq(:,:,igcm_haze))  &
+                                 /ptimestep
          else
-         call hazecloud(ngrid,nlayer,nq,ptimestep, &
-            pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,zdqhaze,   &
-            zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
+            call hazecloud(ngrid,nlayer,nq,ptimestep, &
+               pplay,pplev,pq,pdq,dist_star,mu0,zfluxuv,zdqhaze,   &
+               zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
          endif
 
@@ -1589 +1625 @@
       ! todo: should be placed in haze and use tendency of n2 instead of flusurf
       IF (source_haze) THEN
+         write(*,*) "Source haze not supported yet."
+         stop
             !  call hazesource(ngrid,nlayer,nq,ptimestep,  &
             !                 pplev,flusurf,mu0,zdq_source)