Changeset 4153 for trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90

Timestamp:

Mar 25, 2026, 2:01:37 PM (6 days ago)

Author:

aarfaux

Message:

Titan PCM:

Add a correction to the shortwave heating rates based on htrdr calculations, and
additional related modifications.

File:

• trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90 (modified) (9 diffs)

trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90

@@ -1 +1 @@
       subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,zday,      &
-          albedo,albedo_equivalent,emis,mu0,pplev,pplay,zzlev, &
+          albedo,albedo_equivalent,emis,mu0,                   &
+          pplev,pplay,zzlev,zzlay,                             &
           pt,tsurf,fract,dist_star,                            &
           dtlw,dtsw,fluxsurf_lw,                               &
@@ -7 +8 @@
           OLR_nu,OSR_nu,                                       &
           int_dtaui,int_dtauv,popthi,popthv,poptti,popttv,     &
-          lastcall)
+          lastcall, zlss, zls)
 
       use mod_phys_lmdz_para, only : is_master
@@ -15 +16 @@
       USE tracer_h
       use comcstfi_mod, only: pi, mugaz, cpp
-      use callkeys_mod, only: global1d, szangle,                      &
+      use callkeys_mod, only: global1d, szangle, updatecorrhtrdr,           &
                               diurnal,tracer,seashaze,corrk_recombin, &
                               strictboundcorrk,specOLR,diagdtau,      &
                               tplanckmin,tplanckmax,callclouds,Fcloudy
       use geometry_mod, only: latitude
+      use phys_state_var_mod, only: htrdr_dtauv, htrdr_ssav, htrdr_asymv
+      use hrcorr_mod, only: hr_write_oprop_nc, hr_optcv
 
       implicit none
@@ -67 +70 @@
       REAL,INTENT(IN) :: pplay(ngrid,nlayer)       ! Mid-layer pressure (Pa).
       REAL,INTENT(IN) :: zzlev(ngrid,nlayer+1)     ! Altitude at the atmospheric layer boundaries (ref : local surf).
+      REAL,INTENT(IN) :: zzlay(ngrid,nlayer)       ! Mid-layer altitude
       REAL,INTENT(IN) :: pt(ngrid,nlayer)          ! Air temperature (K).
       REAL,INTENT(IN) :: tsurf(ngrid)              ! Surface temperature (K).
@@ -72 +76 @@
       REAL,INTENT(IN) :: dist_star                 ! Distance star-planet (AU).
       logical,intent(in) :: lastcall               ! Signals last call to physics.
+      REAL,INTENT(IN) :: zlss                      ! sub-solar longitude
+      REAL,INTENT(IN) :: zls                       ! solar longitude
       
       ! OUTPUT
@@ -140 +146 @@
       ! Optical diagnostics
       ! Haze
-      REAL*8 diag_opthi(L_LEVELS,L_NSPECTI,6)
-      REAL*8 diag_opthv(L_LEVELS,L_NSPECTV,6)
+      REAL*8 diag_opthi(L_LEVELS,L_NSPECTI,3)
+      REAL*8 diag_opthv(L_LEVELS,L_NSPECTV,3)
       ! Clouds
       REAL*8 diag_optti(L_LEVELS,L_NSPECTI,3)
@@ -199 +205 @@
 !=======================================================================
 
+      if (lastcall .and. updatecorrhtrdr) then
+          ALLOCATE(htrdr_dtauv(ngrid,L_NLAYRAD,L_NSPECTV,L_NGAUSS))
+          ALLOCATE(htrdr_ssav(ngrid,L_NLAYRAD,L_NSPECTV,L_NGAUSS))
+          ALLOCATE(htrdr_asymv(ngrid,L_NLAYRAD,L_NSPECTV,L_NGAUSS))
+      endif
+         
 
       ! How much light do we get ?
@@ -503 +515 @@
 !-----------------------------------------------------------------------
 
+         if (lastcall .and. updatecorrhtrdr) then
+             ! Clear column :
+             cdcolumn = 0
+             call hr_optcv(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tmid,pmid,    &
+                  dtauv_cc,tauv_cc,taucumv_cc,wbarv_cc,cosbv_cc,tauray,taugsurf,seashazefact,&
+                  diag_opthv,diag_opttv_cc,cdcolumn)
+             ! Dark column :
+             cdcolumn = 1
+             call hr_optcv(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tmid,pmid,    &
+                  dtauv_dc,tauv_dc,taucumv_dc,wbarv_dc,cosbv_dc,tauray,taugsurf,seashazefact,&
+                  diag_opthv,diag_opttv_dc,cdcolumn)
+
+             ! Mean opacity, ssa and asf :
+             where (dtauv_cc(:,:,:) .le. 100 .and. dtauv_dc(:,:,:) .le. 100.)
+                dtauv(:,:,:) = (1-Fcloudy)*dtauv_cc(:,:,:) + Fcloudy*dtauv_dc(:,:,:)
+             elsewhere
+                dtauv(:,:,:) = dtauv_dc(:,:,:) ! No need to average...
+             endwhere
+             do ng = 1, L_NGAUSS
+                do nw = 1, L_NSPECTV
+                   taucumv(1,nw,ng) = 0.0d0
+                   do k = 2, L_LEVELS
+                      if ((taucumv_cc(k,nw,ng).le.100.) .and. (taucumv_dc(k,nw,ng).le.100.)) then
+                         taucumv(k,nw,ng) = taucumv(k-1,nw,ng) + (1-Fcloudy)*taucumv_cc(k,nw,ng) + Fcloudy*taucumv_dc(k,nw,ng)
+                      else
+                         taucumv(k,nw,ng) = taucumv(k-1,nw,ng) + taucumv_dc(k,nw,ng) ! No need to average...
+                      end if
+                   end do
+                   do l = 1, L_NLAYRAD
+                      tauv(l,nw,ng) = taucumv(2*l,nw,ng)
+                   end do
+                   tauv(l,nw,ng) = taucumv(2*L_NLAYRAD+1,nw,ng)
+                end do
+             end do
+
+             wbarv = ((1-Fcloudy) * wbarv_cc*dtauv_cc            + Fcloudy * wbarv_dc *dtauv_dc)
+             wbarv = wbarv /((1-Fcloudy) * dtauv_cc              + Fcloudy * dtauv_dc  + 1.e-30)
+
+             cosbv = ((1-Fcloudy) * cosbv_cc * wbarv_cc*dtauv_cc + Fcloudy * cosbv_dc  * wbarv_dc *dtauv_dc)
+             cosbv = cosbv /((1-Fcloudy) * wbarv_cc*dtauv_cc     + Fcloudy * wbarv_dc *dtauv_dc + 1.e-30)
+             !------------------------------------------------------------------------------
+             htrdr_dtauv(ig,:,:,:) = dtauv(L_NLAYRAD:1:-1,:,:)
+             htrdr_ssav(ig,:,:,:)  = wbarv(L_NLAYRAD:1:-1,:,:)
+             htrdr_asymv(ig,:,:,:) = cosbv(L_NLAYRAD:1:-1,:,:)
+         endif
+
          ! Clear column :
          cdcolumn = 0
@@ -863 +921 @@
 
       if (lastcall) then
+
+         if(updatecorrhtrdr) then
+             call hr_write_oprop_nc(ngrid, nlayer, zzlev, zzlay, pplev, tsurf, modulo(zlss,2*pi), zls)
+         endif
+
         IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
         IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )