Changeset 4153


Ignore:
Timestamp:
Mar 25, 2026, 2:01:37 PM (7 days ago)
Author:
aarfaux
Message:

Titan PCM:

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

Location:
trunk/LMDZ.TITAN/libf/phytitan
Files:
1 added
8 edited

Legend:

Unmodified
Added
Removed

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

    r4148 r4153  
    11      subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,zday,      &
    2           albedo,albedo_equivalent,emis,mu0,pplev,pplay,zzlev, &
     2          albedo,albedo_equivalent,emis,mu0,                   &
     3          pplev,pplay,zzlev,zzlay,                             &
    34          pt,tsurf,fract,dist_star,                            &
    45          dtlw,dtsw,fluxsurf_lw,                               &
    … …  
    78          OLR_nu,OSR_nu,                                       &
    89          int_dtaui,int_dtauv,popthi,popthv,poptti,popttv,     &
    9           lastcall)
     10          lastcall, zlss, zls)
    1011
    1112      use mod_phys_lmdz_para, only : is_master
    … …  
    1516      USE tracer_h
    1617      use comcstfi_mod, only: pi, mugaz, cpp
    17       use callkeys_mod, only: global1d, szangle,                      &
     18      use callkeys_mod, only: global1d, szangle, updatecorrhtrdr,           &
    1819                              diurnal,tracer,seashaze,corrk_recombin, &
    1920                              strictboundcorrk,specOLR,diagdtau,      &
    2021                              tplanckmin,tplanckmax,callclouds,Fcloudy
    2122      use geometry_mod, only: latitude
     23      use phys_state_var_mod, only: htrdr_dtauv, htrdr_ssav, htrdr_asymv
     24      use hrcorr_mod, only: hr_write_oprop_nc, hr_optcv
    2225
    2326      implicit none
    … …  
    6770      REAL,INTENT(IN) :: pplay(ngrid,nlayer)       ! Mid-layer pressure (Pa).
    6871      REAL,INTENT(IN) :: zzlev(ngrid,nlayer+1)     ! Altitude at the atmospheric layer boundaries (ref : local surf).
     72      REAL,INTENT(IN) :: zzlay(ngrid,nlayer)       ! Mid-layer altitude
    6973      REAL,INTENT(IN) :: pt(ngrid,nlayer)          ! Air temperature (K).
    7074      REAL,INTENT(IN) :: tsurf(ngrid)              ! Surface temperature (K).
    … …  
    7276      REAL,INTENT(IN) :: dist_star                 ! Distance star-planet (AU).
    7377      logical,intent(in) :: lastcall               ! Signals last call to physics.
     78      REAL,INTENT(IN) :: zlss                      ! sub-solar longitude
     79      REAL,INTENT(IN) :: zls                       ! solar longitude
    7480     
    7581      ! OUTPUT
    … …  
    140146      ! Optical diagnostics
    141147      ! Haze
    142       REAL*8 diag_opthi(L_LEVELS,L_NSPECTI,6)
    143       REAL*8 diag_opthv(L_LEVELS,L_NSPECTV,6)
     148      REAL*8 diag_opthi(L_LEVELS,L_NSPECTI,3)
     149      REAL*8 diag_opthv(L_LEVELS,L_NSPECTV,3)
    144150      ! Clouds
    145151      REAL*8 diag_optti(L_LEVELS,L_NSPECTI,3)
    … …  
    199205!=======================================================================
    200206
     207      if (lastcall .and. updatecorrhtrdr) then
     208          ALLOCATE(htrdr_dtauv(ngrid,L_NLAYRAD,L_NSPECTV,L_NGAUSS))
     209          ALLOCATE(htrdr_ssav(ngrid,L_NLAYRAD,L_NSPECTV,L_NGAUSS))
     210          ALLOCATE(htrdr_asymv(ngrid,L_NLAYRAD,L_NSPECTV,L_NGAUSS))
     211      endif
     212         
    201213
    202214      ! How much light do we get ?
    … …  
    503515!-----------------------------------------------------------------------
    504516
     517         if (lastcall .and. updatecorrhtrdr) then
     518             ! Clear column :
     519             cdcolumn = 0
     520             call hr_optcv(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tmid,pmid,    &
     521                  dtauv_cc,tauv_cc,taucumv_cc,wbarv_cc,cosbv_cc,tauray,taugsurf,seashazefact,&
     522                  diag_opthv,diag_opttv_cc,cdcolumn)
     523             ! Dark column :
     524             cdcolumn = 1
     525             call hr_optcv(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tmid,pmid,    &
     526                  dtauv_dc,tauv_dc,taucumv_dc,wbarv_dc,cosbv_dc,tauray,taugsurf,seashazefact,&
     527                  diag_opthv,diag_opttv_dc,cdcolumn)
     528
     529             ! Mean opacity, ssa and asf :
     530             where (dtauv_cc(:,:,:) .le. 100 .and. dtauv_dc(:,:,:) .le. 100.)
     531                dtauv(:,:,:) = (1-Fcloudy)*dtauv_cc(:,:,:) + Fcloudy*dtauv_dc(:,:,:)
     532             elsewhere
     533                dtauv(:,:,:) = dtauv_dc(:,:,:) ! No need to average...
     534             endwhere
     535             do ng = 1, L_NGAUSS
     536                do nw = 1, L_NSPECTV
     537                   taucumv(1,nw,ng) = 0.0d0
     538                   do k = 2, L_LEVELS
     539                      if ((taucumv_cc(k,nw,ng).le.100.) .and. (taucumv_dc(k,nw,ng).le.100.)) then
     540                         taucumv(k,nw,ng) = taucumv(k-1,nw,ng) + (1-Fcloudy)*taucumv_cc(k,nw,ng) + Fcloudy*taucumv_dc(k,nw,ng)
     541                      else
     542                         taucumv(k,nw,ng) = taucumv(k-1,nw,ng) + taucumv_dc(k,nw,ng) ! No need to average...
     543                      end if
     544                   end do
     545                   do l = 1, L_NLAYRAD
     546                      tauv(l,nw,ng) = taucumv(2*l,nw,ng)
     547                   end do
     548                   tauv(l,nw,ng) = taucumv(2*L_NLAYRAD+1,nw,ng)
     549                end do
     550             end do
     551
     552             wbarv = ((1-Fcloudy) * wbarv_cc*dtauv_cc            + Fcloudy * wbarv_dc *dtauv_dc)
     553             wbarv = wbarv /((1-Fcloudy) * dtauv_cc              + Fcloudy * dtauv_dc  + 1.e-30)
     554
     555             cosbv = ((1-Fcloudy) * cosbv_cc * wbarv_cc*dtauv_cc + Fcloudy * cosbv_dc  * wbarv_dc *dtauv_dc)
     556             cosbv = cosbv /((1-Fcloudy) * wbarv_cc*dtauv_cc     + Fcloudy * wbarv_dc *dtauv_dc + 1.e-30)
     557             !------------------------------------------------------------------------------
     558             htrdr_dtauv(ig,:,:,:) = dtauv(L_NLAYRAD:1:-1,:,:)
     559             htrdr_ssav(ig,:,:,:)  = wbarv(L_NLAYRAD:1:-1,:,:)
     560             htrdr_asymv(ig,:,:,:) = cosbv(L_NLAYRAD:1:-1,:,:)
     561         endif
     562
    505563         ! Clear column :
    506564         cdcolumn = 0
    … …  
    863921
    864922      if (lastcall) then
     923
     924         if(updatecorrhtrdr) then
     925             call hr_write_oprop_nc(ngrid, nlayer, zzlev, zzlay, pplev, tsurf, modulo(zlss,2*pi), zls)
     926         endif
     927
    865928        IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
    866929        IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )

  • trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90

    r4047 r4153  
    44      logical,save :: callrad,corrk,calldifv,UseTurbDiff
    55!$OMP THREADPRIVATE(callrad,corrk,calldifv,UseTurbDiff)
     6      logical,save :: corrhtrdr, updatecorrhtrdr
     7!$OMP THREADPRIVATE(corrhtrdr,updatecorrhtrdr)
    68      logical,save :: calladj,callsoil
    79!$OMP THREADPRIVATE(calladj,callsoil)

  • trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90

    r4047 r4153  
    297297     write(*,*) " corrk = ",corrk
    298298     
     299     write(*,*) "call the 3D correction of the sw heating rate ?"
     300     corrhtrdr=.false. ! default value
     301     call getin_p("corrhtrdr",corrhtrdr)
     302     write(*,*) " corrhtrdr = ",corrhtrdr
     303
     304     write(*,*) "update the correction factor every week ?"
     305     corrhtrdr=.false. ! default value
     306     call getin_p("updatecorrhtrdr",updatecorrhtrdr)
     307     write(*,*) " updatecorrhtrdr = ",updatecorrhtrdr
     308
    299309     if (corrk) then
    300310       ! default path is set in datadir

  • trunk/LMDZ.TITAN/libf/phytitan/optci.F90

    r4013 r4153  
    400400                  enddo
    401401
    402                   elseif(igas.eq.igas_CH4)then
     402               elseif(igas.eq.igas_CH4)then
    403403
    404404                  ! first do self-induced absorption

  • trunk/LMDZ.TITAN/libf/phytitan/optcv.F90

    r4013 r4153  
    133133  logical,save :: firstcall=.true.
    134134!$OMP THREADPRIVATE(firstcall)
    135 
    136135
    137136  !! AS: to save time in computing continuum (see bilinearbig)
    … …  
    537536   DO L=1,L_NLAYRAD-1
    538537      K = 2*L+1
    539            atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW) + ASF_T(K+1,NW)*DHAZES_T(K+1,NW)
     538      atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW) + ASF_T(K+1,NW)*DHAZES_T(K+1,NW)
    540539      btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
    541            ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ?
    542            btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
     540      ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ?
     541      btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
    543542      COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
    544543   END DO
    … …  
    568567      L = L_NLAYRAD
    569568      K = 2*L+1
    570            DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
     569      DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
    571570      WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG)
    572571   END DO

  • trunk/LMDZ.TITAN/libf/phytitan/phyredem.F90

    r1903 r4153  
    6868  tab_cntrl(10) = daysec  ! length of a sol (s)  ~88775
    6969
    70   tab_cntrl(11) = phystep  ! time step in the physics
     70  tab_cntrl(11) = phystep ! time step in the physics
    7171  tab_cntrl(12) = 0.
    7272  tab_cntrl(13) = 0.
    … …  
    9393  tab_cntrl(34) = dtemisice(2) ! time scale for snow metamorphism (south)
    9494
    95   tab_cntrl(28) = 0. 
     95  tab_cntrl(28) = 0.
    9696  tab_cntrl(29) = 0.
    9797  tab_cntrl(30) = 0.
    … …  
    186186  ! Methane tank depth
    187187  call put_field("tankCH4","Depth of methane tank",tankCH4)
    188  
     188
    189189  ! Tracers
    190190  if (nq>0) then

  • trunk/LMDZ.TITAN/libf/phytitan/phys_state_var_mod.F90

    r4013 r4153  
    7777      real,dimension(:,:,:,:),allocatable,save :: zpopttv ! VI optical properties [haze+clouds] within narrowbands for diags (dtau,tau,k,wbar,gbar,drayaer,taugaz,dcont).
    7878!$OMP THREADPRIVATE(zpopthi,zpopthv,zpoptti,zpopttv)
     79      real,dimension(:,:,:,:),allocatable,save :: htrdr_dtauv ! opacity of layers for htrdr coupling
     80      real,dimension(:,:,:,:),allocatable,save :: htrdr_ssav  ! single scattering albedo for htrdr coupling
     81      real,dimension(:,:,:,:),allocatable,save :: htrdr_asymv ! essymetrie parameter for htrdr coupling
     82!$OMP THREADPRIVATE(htrdr_dtauv,htrdr_ssav,htrdr_asymv)
    7983
    8084      real,dimension(:,:),allocatable,save :: u_ref ! Reference profile for the zonal wind nudging (m.s-1).
    … …  
    209213        DEALLOCATE(zpoptti)
    210214        DEALLOCATE(zpopttv)
     215        DEALLOCATE(htrdr_dtauv)
     216        DEALLOCATE(htrdr_ssav)
     217        DEALLOCATE(htrdr_asymv)
    211218        DEALLOCATE(u_ref)
    212219        DEALLOCATE(dycchi)

  • trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90

    r4051 r4153  
    5959#endif
    6060      use muphy_diag
     61      use hrcorr_mod
    6162      implicit none
    6263
    … …  
    248249      real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine.
    249250      real zdtlc(ngrid,nlayer)                        ! Condensation heating rate.
     251      real refCorr, time                              ! for hrcorr_mod routines.
    250252                             
    251253      ! For Surface Tracers : (kg/m2/s)
    … …  
    464466           call setspv            ! Basic visible properties.
    465467           call sugas_corrk       ! Set up gaseous absorption properties.
    466          
     468
    467469           OLR_nu(:,:) = 0.D0
    468470           OSR_nu(:,:) = 0.D0
    … …  
    483485              endif
    484486           ENDIF
     487
     488           if (corrhtrdr) then
     489               call ini_hrcorr
     490           endif
    485491!#endif           
    486492
    … …  
    820826      ! omega in Pa/s
    821827      do l=1,nlayer-1
    822          omega(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
    823        enddo
    824        omega(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
    825        do l=1,nlayer
    826          omega(1:ngrid,l)=g*omega(1:ngrid,l)/cell_area(1:ngrid)
    827        enddo
     828          omega(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1))
     829      enddo
     830          omega(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0
     831      do l=1,nlayer
     832          omega(1:ngrid,l)=g*omega(1:ngrid,l)/cell_area(1:ngrid)
     833      enddo
    828834
    829835!---------------------------------
    … …  
    876882               ! standard callcorrk
    877883               call callcorrk(ngrid,nlayer,pq,nq,qsurf,zday,                      &
    878                               albedo,albedo_equivalent,emis,mu0,pplev,pplay,zzlev,&
     884                              albedo,albedo_equivalent,emis,mu0,                  &
     885                              pplev,pplay,zzlev,zzlay,                            &
    879886                              pt,tsurf,fract,dist_star,                           &
    880887                              zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,                &
    … …  
    882889                              fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,                &
    883890                              int_dtaui,int_dtauv,zpopthi,zpopthv,zpoptti,zpopttv,&
    884                               lastcall)
     891                              lastcall,zlss,zls)
    885892
    886893               ! Radiative flux from the sky absorbed by the surface (W.m-2).
    … …  
    895902
    896903               ! Net atmospheric radiative heating rate (K.s-1)
     904               if(corrhtrdr) then
     905                   call hr_average(ngrid, nlayer, pplay, zzlev, zdtsw, refCorr)
     906                   do ig=1,ngrid
     907                       time = modulo(ptime+longitude(ig)/(2.*pi), 1.)
     908                       if(is_master) write(*,*) longitude(ig)*180./pi, time
     909                       call apply_hrcorr(ig, nlayer, zls*180./pi, time , zdtsw(ig,:), refCorr)
     910                   enddo
     911                   if (lastcall) then
     912                       call deallocate_hrcorr
     913                   endif
     914               endif
    897915               dtrad(:,:)=zdtsw(:,:)+zdtlw(:,:)
    898                
     916
    899917            elseif(newtonian)then
    900918           
    … …  
    973991         ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff.
    974992         if (UseTurbDiff) then
    975          
     993
    976994            call turbdiff(ngrid,nlayer,nq,                       &
    977995                          ptimestep,capcal,lwrite,               &
Note: See TracChangeset for help on using the changeset viewer.