Changeset 3585 for trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90

Timestamp:

Jan 20, 2025, 10:17:54 AM (12 days ago)

Author:

debatzbr

Message:

Connecting microphysics to radiative transfer + miscellaneous cleans

File:

• trunk/LMDZ.PLUTO/libf/phypluto/callcorrk.F90 (modified) (11 diffs)

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

@@ -7 +7 @@
       subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
           albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,    &
-          zzlay,tsurf,fract,dist_star,aerosol,muvar,           &
+          zzlay,zzlev,tsurf,fract,dist_star,dtau_aer,muvar,    &
           dtlw,dtsw,fluxsurf_lw,                               &
           fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,               &
@@ -36 +36 @@
                               optichaze,haze_radproffix,&
                               methane,carbox,cooling,nlte,strobel,&
-                              ch4fix,vmrch4_proffix,vmrch4fix
+                              ch4fix,vmrch4_proffix,vmrch4fix,&
+                              callmufi
       use optcv_mod, only: optcv
       use optci_mod, only: optci
@@ -44 +45 @@
       use planetwide_mod, only: planetwide_maxval, planetwide_minval
       use radcommon_h, only: wavev,wavei
+      use mp2m_diagnostics
       implicit none
 
@@ -81 +83 @@
       REAL,INTENT(IN) :: pt(ngrid,nlayer)          ! Air temperature (K).
       REAL,INTENT(IN) :: zzlay(ngrid,nlayer)       ! Mid-layer altitude
+      REAL,INTENT(IN) :: zzlev(ngrid,nlayer)       ! Altitude at the layer boundaries.
       REAL,INTENT(IN) :: tsurf(ngrid)              ! Surface temperature (K).
       REAL,INTENT(IN) :: fract(ngrid)              ! Fraction of day.
@@ -89 +92 @@
 
       ! OUTPUT
-      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau at reference wavelenght.
-      REAL,INTENT(OUT) :: dtlw(ngrid,nlayer)             ! Heating rate (K/s) due to LW radiation.
-      REAL,INTENT(OUT) :: dtsw(ngrid,nlayer)             ! Heating rate (K/s) due to SW radiation.
-      REAL,INTENT(OUT) :: fluxsurf_lw(ngrid)             ! Incident LW flux to surf (W/m2).
-      REAL,INTENT(OUT) :: fluxsurf_sw(ngrid)             ! Incident SW flux to surf (W/m2)
-      REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid)          ! Absorbed SW flux by the surface (W/m2). By MT2015.
-      REAL,INTENT(OUT) :: fluxtop_lw(ngrid)              ! Outgoing LW flux to space (W/m2).
-      REAL,INTENT(OUT) :: fluxabs_sw(ngrid)              ! SW flux absorbed by the planet (W/m2).
-      REAL,INTENT(OUT) :: fluxtop_dn(ngrid)              ! Incident top of atmosphere SW flux (W/m2).
-      REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI)        ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1).
-      REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV)        ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1).
-      REAL,INTENT(OUT) :: GSR_nu(ngrid,L_NSPECTV)        ! Surface SW radiation in each band (Normalized to the band width (W/m2/cm-1).
-      REAL,INTENT(OUT) :: tau_col(ngrid)                 ! Diagnostic from aeropacity.
-      REAL,INTENT(OUT) :: albedo_equivalent(ngrid)       ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
+      REAL,INTENT(OUT) :: dtau_aer(ngrid,nlayer,naerkind) ! Aerosol tau at reference wavelenght.
+      REAL,INTENT(OUT) :: dtlw(ngrid,nlayer)              ! Heating rate (K/s) due to LW radiation.
+      REAL,INTENT(OUT) :: dtsw(ngrid,nlayer)              ! Heating rate (K/s) due to SW radiation.
+      REAL,INTENT(OUT) :: fluxsurf_lw(ngrid)              ! Incident LW flux to surf (W/m2).
+      REAL,INTENT(OUT) :: fluxsurf_sw(ngrid)              ! Incident SW flux to surf (W/m2)
+      REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid)           ! Absorbed SW flux by the surface (W/m2). By MT2015.
+      REAL,INTENT(OUT) :: fluxtop_lw(ngrid)               ! Outgoing LW flux to space (W/m2).
+      REAL,INTENT(OUT) :: fluxabs_sw(ngrid)               ! SW flux absorbed by the planet (W/m2).
+      REAL,INTENT(OUT) :: fluxtop_dn(ngrid)               ! Incident top of atmosphere SW flux (W/m2).
+      REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI)         ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1).
+      REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV)         ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1).
+      REAL,INTENT(OUT) :: GSR_nu(ngrid,L_NSPECTV)         ! Surface SW radiation in each band (Normalized to the band width (W/m2/cm-1).
+      REAL,INTENT(OUT) :: tau_col(ngrid)                  ! Diagnostic from aeropacity.
+      REAL,INTENT(OUT) :: albedo_equivalent(ngrid)        ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
       REAL,INTENT(OUT) :: int_dtaui(ngrid,nlayer,L_NSPECTI) ! VI optical thickness of layers within narrowbands for diags ().
       REAL,INTENT(OUT) :: int_dtauv(ngrid,nlayer,L_NSPECTV) ! IR optical thickness of layers within narrowbands for diags ().
@@ -186 +189 @@
 !$OMP THREADPRIVATE(QREFvis3d,QREFir3d)
 
-
       ! Miscellaneous :
-      real*8  temp,temp1,temp2,pweight
+      real*8  temp,temp1,temp2,pweight,sig
       character(len=10) :: tmp1
       character(len=10) :: tmp2
@@ -331 +333 @@
          if (is_master) call system('rm -f surf_vals_long.out')
 
-         call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
-
-
 !--------------------------------------------------
 !             Set up correlated k
@@ -454 +453 @@
 !     Effective radius and variance of the aerosols
 !--------------------------------------------------
+      ! Radiative Hazes
       if (optichaze) then
-         do iaer=1,naerkind
-            if ((iaer.eq.iaero_haze)) then
-               call su_aer_radii(ngrid,nlayer,reffrad(1,1,iaer), &
-                            nueffrad(1,1,iaer))
+         if (callmufi) then
+            ! Spherical aerosols
+            sig = 0.2
+            where (mp2m_rc_sph(:,:) > 1e-10)
+               reffrad(:,:,1) = mp2m_rc_sph(:,:) * exp(5.*sig**2 / 2.)
+            elsewhere
+               reffrad(:,:,1) = 0d0     
+            endwhere
+            if (exp(sig**2) - 1 > 0.1) then
+               nueffrad(:,:,1) = exp(sig**2) - 1
+            else
+               nueffrad(:,:,1) = 0.1    
             endif
-         end do !iaer=1,naerkind.
-         if (haze_radproffix) then
+            ! Fractal aerosols
+            sig = 0.35
+            where (mp2m_rc_fra(:,:) > 1e-10)
+               reffrad(:,:,2) = mp2m_rc_fra(:,:) * exp(5.*sig**2 / 2.)
+               elsewhere
+               reffrad(:,:,2) = 0d0
+               endwhere
+            if (exp(sig**2) - 1 > 0.1) then
+               nueffrad(:,:,2) = exp(sig**2) - 1
+            else
+               nueffrad(:,:,2) = 0.1    
+            endif
+           
+         else
+            do iaer=1,naerkind
+               if ((iaer.eq.iaero_haze)) then
+               call su_aer_radii(ngrid,nlayer,reffrad(1,1,iaer),nueffrad(1,1,iaer))
+               endif
+            end do ! iaer = 1, naerkind.
             if (haze_radproffix) then
-                call haze_reffrad_fix(ngrid,nlayer,zzlay,&
-                    reffrad,nueffrad)
-            endif
-
-            print*, 'haze_radproffix=T : fixed profile for haze rad'
-         else
-            print*,'reffrad haze:',reffrad(1,1,iaero_haze)
-            print*,'nueff haze',nueffrad(1,1,iaero_haze)
-         endif
-      endif
-
+               call haze_reffrad_fix(ngrid,nlayer,zzlay,reffrad,nueffrad)
+               if (is_master) print*, 'haze_radproffix=T : fixed profile for haze rad'
+            else
+               if (is_master) print*,'reffrad haze:',reffrad(1,1,iaero_haze)
+               if (is_master) print*,'nueff haze',nueffrad(1,1,iaero_haze)
+            endif ! end haze_radproffix
+         endif ! end callmufi
+      endif ! end radiative haze
 
       ! How much light do we get ?
@@ -489 +511 @@
 
          ! Get aerosol optical depths.
-         call aeropacity(ngrid,nlayer,nq,pplay,pplev, pt,pq,aerosol,      &
+         call aeropacity(ngrid,nlayer,nq,pplay,pplev,zzlev,pt,pq,dtau_aer,      &
             reffrad,nueffrad,QREFvis3d,QREFir3d,                             &
             tau_col)
@@ -665 +687 @@
                pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/   &
                        (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))
-               ! As 'aerosol' is at reference (visible) wavelenght we scale it as
+               ! As 'dtau_aer' is at reference (visible) wavelenght we scale it as
                ! it will be multplied by qxi/v in optci/v
-               temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
+               temp=dtau_aer(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer)
                tauaero(2*k+2,iaer)=max(temp*pweight,0.d0)
                tauaero(2*k+3,iaer)=max(temp-tauaero(2*k+2,iaer),0.d0)
@@ -986 +1008 @@
             print*,'fluxtop_dn=',fluxtop_dn(ig)
             print*,'acosz=',acosz
-            print*,'aerosol=',aerosol(ig,:,:)
+            print*,'dtau_aer=',dtau_aer(ig,:,:)
             print*,'temp=   ',pt(ig,:)
             print*,'pplay=  ',pplay(ig,:)