Index: trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90 (revision 4152) +++ trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90 (revision 4153) @@ -1,4 +1,5 @@ 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,5 +8,5 @@ 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,9 +16,11 @@ 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,4 +70,5 @@ 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,4 +76,6 @@ 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,6 +146,6 @@ ! 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,4 +205,10 @@ !======================================================================= + 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,4 +515,50 @@ !----------------------------------------------------------------------- + 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,4 +921,9 @@ 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 ) Index: trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90 (revision 4152) +++ trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90 (revision 4153) @@ -4,4 +4,6 @@ logical,save :: callrad,corrk,calldifv,UseTurbDiff !$OMP THREADPRIVATE(callrad,corrk,calldifv,UseTurbDiff) + logical,save :: corrhtrdr, updatecorrhtrdr +!$OMP THREADPRIVATE(corrhtrdr,updatecorrhtrdr) logical,save :: calladj,callsoil !$OMP THREADPRIVATE(calladj,callsoil) Index: trunk/LMDZ.TITAN/libf/phytitan/hrcorr_mod.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/hrcorr_mod.F90 (revision 4153) +++ trunk/LMDZ.TITAN/libf/phytitan/hrcorr_mod.F90 (revision 4153) @@ -0,0 +1,1159 @@ +MODULE hrcorr_mod + IMPLICIT NONE +PRIVATE + INTEGER, SAVE :: NALT ! number of altitudes + INTEGER, SAVE :: NLAT ! number of latitudes + INTEGER, SAVE :: NLS ! number of solar longitudes + INTEGER, SAVE :: NT ! number of times +!$OMP THREADPRIVATE(NALT,NLAT,NLS) + REAL, SAVE, ALLOCATABLE :: tau_MC_2S(:,:,:,:) ! table containing the correction, shape (NALT,NLAT,NLS,NT) + REAL, SAVE, ALLOCATABLE :: solLon(:) ! solar longitudes, shape (NLS) + REAL, SAVE, ALLOCATABLE :: latCorr(:) ! latitudes of the correction table, shape (NLAT) + REAL, SAVE, ALLOCATABLE :: times(:) ! latitudes of the correction table, shape (NT) +!$OMP THREADPRIVATE(tau_MC_2S,solLon,latCorr) + REAL, SAVE :: intLatInf ! lower boundary for latitudinal band integral + REAL, SAVE :: intLatSup ! upper boundary for latitudinal band integral +!$OMP THREADPRIVATE(intLatInf,intLatSup) + REAL, SAVE :: intPrsInf ! lower boundary for pressure integral + REAL, SAVE :: intPrsSup ! upper boundary for pressure integral +!$OMP THREADPRIVATE(intPrsInf,intPrsSup) + CHARACTER(len=30) :: corr_method ! correction method to use +!$OMP THREADPRIVATE(corr_method) +PUBLIC :: ini_hrcorr, apply_hrcorr, hr_average, deallocate_hrcorr, hr_write_oprop_nc, hr_optcv + +!=============================================================================== +! +! Purpose +! ------- +! This module provides the necessary tools to apply a correction on the heating rate. +! The correction should be contained in a file that must be read iby 'ini_hrcorr' at the initilization step. +! The correction can be applied by 'apply_hrcorr'. +! In the eventuality of using the 'refEq' method, the 'hr_average' subroutine calculates the reference value by averaging the +! 2-stream heating rate. +! The routine 'hr_optcv' is a copy of the 'optcv' routine, but with a half-layer shift, as it seems that the original behavior of +! The routine is wrong (but compensated for in the radiative transfer routine). +! The routine 'hr_write_oprop_nc' is called by callcorrk at the end of every week (if the updatecorrhtrdr key is set to .true.) to +! write an 'optical_properties.nc' file in the working directory. This can be used to update the correction lookup table. +! Finally, 'deallocate_hrcorr' deallocates the allocated tables +! +! Author +! ------ +! Anthony Arfaux (03/2026) +! +!=============================================================================== + + CONTAINS + + subroutine ini_hrcorr + use datafile_mod, only: datadir + + implicit none + +!=============================================================================== +! +! Purpose +! ------- +! Initialize the module by reading the 'hr_corr.txt' input file, +! which contains the necessary information for the correction. +! Allocates the tables and assign the values to the private parameters. +! The file is "hr_corr.txt" and should be placed in the datadir repository. +! +! Author +! ------ +! Anthony Arfaux (03/2026) +! +!=============================================================================== + + integer :: nLines + integer :: iAlt + integer :: iLat + integer :: iLS + integer :: iT + integer :: i + real :: pres + real :: lat + real :: ls + real :: t + real :: corr + + open(25, file=trim(datadir)//"/hr_corr.txt") + read(25,*) corr_method + read(25,*) nAlt, nLat, nLS, nT + read(25,*) intLatInf, intLatSup + read(25,*) intPrsInf, intPrsSup + read(25,*) + allocate(tau_MC_2S(nAlt,nLat,nLS,nT),solLon(nLS),latCorr(nLat),times(nT)) + + nLines = nAlt * nLat * nLS * nT + + do i=1,nLines + read(25,*) iAlt, iLat, iLS, iT, pres, lat, ls, t, corr + tau_MC_2S(iAlt, iLat, iLS, iT) = corr + solLon(iLS) = ls + latCorr(iLat) = lat + times(iT) = t + enddo + close(25) + end subroutine ini_hrcorr + + subroutine apply_hrcorr(ig, nlayer, ls, t, tau, tau_eq) + use geometry_mod, only: latitude_deg + use mod_phys_lmdz_para, only : is_master + + implicit none + +!=============================================================================== +! +! Purpose +! ------- +! Apply the correction on a column. +! +! INPUT +! ----- +! ig = Index of the column [-] +! nlayer = Number of altitude layer [-] +! ls = Solar longitude [°] +! t = Local time [-] +! tau = Heating rate in column 'iLat' [K.s-1] +! tau_eq = Averaged heating rate for the reference (see 'hr_average') [K.s-1] +! +! OUTPUT +! ------ +! tau = Corrected heating rate in column 'ig' [K.s-1] +! +! Author +! ------ +! Anthony Arfaux (03/2026) +! +!=============================================================================== + + integer, intent(in ) :: ig + integer, intent(in ) :: nlayer + real, intent(in ) :: ls + real, intent(in ) :: t + real, intent(inout) :: tau(nlayer) + real, intent(in ) :: tau_eq + + real, allocatable :: corrLS(:,:,:) + real, allocatable :: corrLat(:,:) + real, allocatable :: corr(:) + real :: lat + integer :: iLS + integer :: iLat + integer :: iT + real :: lsTest + real :: latTest + real :: tTest + real :: alphaLS + real :: alphaLAT + real :: alphaT + + allocate(corrLS(nAlt,nLat,nT),corrLat(nAlt,nT),corr(nAlt)) + + lat = latitude_deg(ig) + + iLS = 1 + lsTest = solLon(iLS + 1) + do while (lsTest < ls) + iLS = iLS + 1 + lsTest = solLon(iLS + 1) + end do + + alphaLS = (ls - solLon(iLS)) / (solLon(iLS + 1) - solLon(iLS)) + corrLS = alphaLS * tau_MC_2S(:,:,iLS+1,:) + (1 - alphaLS) * tau_MC_2S(:,:,iLS,:) + + iLat = 1 + latTest = latCorr(iLat + 1) + do while (latTest < lat) + iLat = iLat + 1 + latTest = latCorr(iLat + 1) + end do + + alphaLAT = (lat - latCorr(iLat)) / (latCorr(iLat + 1) - latCorr(iLat)) + corrLat = alphaLAT * corrLS(:,iLat+1,:) + (1 - alphaLAT) * corrLS(:,iLat,:) + + iT = 1 + tTest = times(iT + 1) + do while (tTest < t) + iT = iT + 1 + tTest = times(iT + 1) + end do + + alphaT = (t - times(iT)) / (times(iT + 1) - times(iT)) + corr = alphaT * corrLAT(:,iT+1) + (1 - alphaT) * corrLAT(:,iT) + + + if (trim(corr_method) == "refEq") then + tau = tau_eq * corr + elseif (trim(corr_method) == "diff") then + tau = tau + corr + elseif (trim(corr_method) == "prod") then + tau = tau * corr + else + stop "Wrong correction method" + endif + + deallocate(corrLS,corr) + end subroutine apply_hrcorr + + subroutine hr_average(ngrid, nlayer, pressure, altitude, dt, dt_average) + use geometry_mod, only: latitude_deg + use comcstfi_mod, only: rad + use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat + use mod_phys_lmdz_transfert_para, only: bcast, reduce_sum + use mod_phys_lmdz_para, only : is_master + + implicit none + +!=============================================================================== +! +! Purpose +! ------- +! Calculate the reference heating rate by integarting the shortwave heating +! rate in a range of latitude and pressure as provided in the 'hr_corr.txt' +! file. +! +! INPUT +! ----- +! ngrid = Number of columns [-] +! nlayer = Number of altitude layer [-] +! pressure = Layer pressure [Pa] +! altitude = Level altitude [m] +! dt = heating rates [K.s-1] +! +! OUTPUT +! ------ +! dt_average = Heating rates averaged according to the reference [K.s-1] +! +! Author +! ------ +! Anthony Arfaux (03/2026) +! +!=============================================================================== + + integer, intent(in) :: ngrid + integer, intent(in) :: nlayer + real, intent(in) :: pressure(ngrid,nlayer) + real, intent(in) :: altitude(ngrid,nlayer+1) + real, intent(in) :: dt(ngrid,nlayer) + real, intent(out) :: dt_average + + real :: pi + real :: deltaLon + real :: deltaLat + integer :: i + integer :: j + real :: volume + real :: radInf + real :: radSup + real :: dt_sum + real :: volume_sum + real :: dt_sum_tot + real :: volume_sum_tot + + dt_average = 0. + dt_sum_tot = 0. + volume_sum_tot = 0. + if(corr_method /= "refEq") return + + ! write(*,*) "@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" + ! write(*,*) "AnthonyAfx" + + pi = 2. * asin(1.) + deltaLon = 2. * pi / nbp_lon + deltaLat = pi / nbp_lat + + dt_sum = 0. + volume_sum = 0. + do i = 1, ngrid + do j = 1, nlayer + if( latitude_deg(i) >= intLatInf .and. latitude_deg(i) <= intLatSup & + .and. pressure(i,j) >= intPrsInf .and. pressure(i,j) <= intPrsSup ) then + radInf = altitude(i,j) + rad + radSup = altitude(i,j+1) + rad + volume = (2. / 3.) * deltaLon * (radSup**3. - radInf**3.) * cos(latitude_deg(i) * pi/180.) * sin(deltaLat/2.) + dt_sum = dt_sum + volume * dt(i,j) + volume_sum = volume_sum + volume + endif + enddo + enddo + + call reduce_sum(dt_sum, dt_sum_tot) + call reduce_sum(volume_sum, volume_sum_tot) + if (is_master) then + dt_average = dt_sum_tot / volume_sum_tot + endif + call bcast(dt_average) + ! write(*,*) "dt_average: ", dt_average + ! dt_average = dt_sum/ volume_sum + end subroutine hr_average + + subroutine hr_write_oprop_nc(ngrid, nlayer, zlev, zlay, plev, tsurf, zlss, ls) + USE netcdf, only: NF90_CREATE, NF90_CLOBBER, NF90_64BIT_OFFSET, & + NF90_NOERR, nf90_strerror, & + NF90_PUT_ATT, NF90_GLOBAL, NF90_DEF_DIM, & + NF90_UNLIMITED, NF90_ENDDEF, & + NF90_WRITE, NF90_OPEN, NF90_CLOSE, & + NF90_DOUBLE, NF90_FLOAT, NF90_PUT_VAR, NF90_DEF_VAR + use radinc_h, only : L_NSPECTV, L_NLAYRAD, L_NGAUSS + use comsaison_h, only: dist_star, declin + use comcstfi_mod, only: rad, g, cpp + use radcommon_h, only: WNOV, DWNV, GWEIGHT + use geometry_mod, only: latitude_deg, longitude_deg + use mod_grid_phy_lmdz, only : klon_glo, nbp_lon, nbp_lat, grid1dTo2d_glo + use phys_state_var_mod, only: htrdr_dtauv, htrdr_ssav, htrdr_asymv, zdtsw, albedo + use mod_phys_lmdz_para, only : is_master + use mod_phys_lmdz_transfert_para, only: gather + + implicit none + +!=============================================================================== +! +! Purpose +! ------- +! Create an output file containing all necessary information to update the +! the correction lookup table. +! +! INPUT +! ----- +! ngrid = Number of columns [-] +! nlayer = Number of altitude layer [-] +! zlev = Level altitude [m] +! zlay = Layer altitude [m] +! plev = Level pressure [Pa] +! tsurf = Surface temperature [K] +! zlss = Sub-solar longitude [rad] +! ls = Solar longitude [rad] +! +! Author +! ------ +! Anthony Arfaux (03/2026) +! +!=============================================================================== + + integer, intent(in) :: ngrid, nlayer + real, intent(in) :: zlay(ngrid,nlayer) + real, intent(in) :: zlev(ngrid,nlayer+1) + real, intent(in) :: plev(ngrid,nlayer+1) + real, intent(in) :: tsurf(ngrid) + real, intent(in) :: zlss + real, intent(in) :: ls + + real :: htrdr_dtauv_glo(klon_glo, L_NLAYRAD, L_NSPECTV, L_NGAUSS) + real :: htrdr_ssav_glo(klon_glo, L_NLAYRAD, L_NSPECTV, L_NGAUSS) + real :: htrdr_asymv_glo(klon_glo, L_NLAYRAD, L_NSPECTV, L_NGAUSS) + real :: surf_albedo_glo(klon_glo, L_NSPECTV) + real :: zdtsw_glo(klon_glo, nlayer) + real :: zlay_glo(klon_glo, nlayer) + real :: zlev_glo(klon_glo, nlayer+1) + real :: plev_glo(klon_glo, nlayer+1) + real :: tsurf_glo(klon_glo) + real :: lat_glo(klon_glo) + real :: lon_glo(klon_glo) + + real :: htrdr_dtauv_glo_lonlat(nbp_lon, nbp_lat, L_NLAYRAD, L_NSPECTV, L_NGAUSS) + real :: htrdr_ssav_glo_lonlat(nbp_lon, nbp_lat, L_NLAYRAD, L_NSPECTV, L_NGAUSS) + real :: htrdr_asymv_glo_lonlat(nbp_lon, nbp_lat, L_NLAYRAD, L_NSPECTV, L_NGAUSS) + real :: surf_albedo_glo_lonlat(nbp_lon, nbp_lat, L_NSPECTV) + real :: zdtsw_glo_lonlat(nbp_lon, nbp_lat, nlayer) + real :: zlay_glo_lonlat(nbp_lon, nbp_lat, nlayer) + real :: zlev_glo_lonlat(nbp_lon, nbp_lat, nlayer+1) + real :: plev_glo_lonlat(nbp_lon, nbp_lat, nlayer+1) + real :: tsurf_glo_lonlat(nbp_lon, nbp_lat) + real :: lat_glo_lonlat(nbp_lon, nbp_lat) + real :: lon_glo_lonlat(nbp_lon, nbp_lat) + + integer :: iLon, iLat, iLay, iLev, iWave, iWeight + integer :: id_dtauv, id_ssav, id_asymv, id_salb + integer :: id_g, id_wave, id_wavebin, id_dtsw + integer :: id_zlay, id_zlev, id_plev, id_lat, id_lon, id_tsurf + integer :: id_dist, id_sslon, id_sslat, id_rad, id_grav, id_cp, id_ls + + INTEGER :: ierr, nid + + character(len=50) :: filename + + filename = "optical_properties.nc" + + call gather(htrdr_dtauv, htrdr_dtauv_glo) + call gather(htrdr_ssav, htrdr_ssav_glo) + call gather(htrdr_asymv, htrdr_asymv_glo) + call gather(albedo, surf_albedo_glo) + call gather(zdtsw, zdtsw_glo) + call gather(zlay, zlay_glo) + call gather(zlev, zlev_glo) + call gather(plev, plev_glo) + call gather(tsurf, tsurf_glo) + call gather(latitude_deg, lat_glo) + call gather(longitude_deg, lon_glo) + + if (is_master) then + call grid1dTo2d_glo(htrdr_dtauv_glo, htrdr_dtauv_glo_lonlat) + call grid1dTo2d_glo(htrdr_ssav_glo, htrdr_ssav_glo_lonlat) + call grid1dTo2d_glo(htrdr_asymv_glo, htrdr_asymv_glo_lonlat) + call grid1dTo2d_glo(surf_albedo_glo, surf_albedo_glo_lonlat) + call grid1dTo2d_glo(zdtsw_glo, zdtsw_glo_lonlat) + call grid1dTo2d_glo(zlay_glo, zlay_glo_lonlat) + call grid1dTo2d_glo(zlev_glo, zlev_glo_lonlat) + call grid1dTo2d_glo(plev_glo, plev_glo_lonlat) + call grid1dTo2d_glo(tsurf_glo, tsurf_glo_lonlat) + call grid1dTo2d_glo(lat_glo, lat_glo_lonlat) + call grid1dTo2d_glo(lon_glo, lon_glo_lonlat) + + ierr=NF90_CREATE(filename,IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), nid) + IF (ierr/=NF90_NOERR) THEN + write(*,*)'hr_write_oprop_nc: problem creating file '//trim(filename) + write(*,*)trim(nf90_strerror(ierr)) + CALL ABORT + ENDIF + + ierr=NF90_DEF_DIM(nid,"longitude",nbp_lon,iLon) + IF (ierr/=NF90_NOERR) THEN + write(*,*)'hr_write_oprop_nc: problem defining longitude dimension ' + write(*,*)trim(nf90_strerror(ierr)) + CALL ABORT + ENDIF + + ierr=NF90_DEF_DIM(nid,"latitude",nbp_lat,iLat) + IF (ierr/=NF90_NOERR) THEN + write(*,*)'hr_write_oprop_nc: problem defining latitude dimension ' + write(*,*)trim(nf90_strerror(ierr)) + CALL ABORT + ENDIF + + ierr=NF90_DEF_DIM(nid,"layer",nlayer,iLay) + IF (ierr/=NF90_NOERR) THEN + write(*,*)'hr_write_oprop_nc: problem defining layer dimension ' + write(*,*)trim(nf90_strerror(ierr)) + CALL ABORT + ENDIF + + ierr=NF90_DEF_DIM(nid,"level",nlayer+1,iLev) + IF (ierr/=NF90_NOERR) THEN + write(*,*)'hr_write_oprop_nc: problem defining level dimension ' + write(*,*)trim(nf90_strerror(ierr)) + CALL ABORT + ENDIF + + ierr=NF90_DEF_DIM(nid,"wavenumber",L_NSPECTV,iWave) + IF (ierr/=NF90_NOERR) THEN + write(*,*)'hr_write_oprop_nc: problem defining wavenumber dimension ' + write(*,*)trim(nf90_strerror(ierr)) + CALL ABORT + ENDIF + + ierr=NF90_DEF_DIM(nid,"gauss_point",L_NGAUSS,iWeight) + IF (ierr/=NF90_NOERR) THEN + write(*,*)'hr_write_oprop_nc: problem defining gauss_point dimension ' + write(*,*)trim(nf90_strerror(ierr)) + CALL ABORT + ENDIF + +#ifdef NC_DOUBLE + ierr=NF90_DEF_VAR(nid, 'ls', NF90_DOUBLE, id_ls ) + ierr=NF90_DEF_VAR(nid, 'dist', NF90_DOUBLE, id_dist ) + ierr=NF90_DEF_VAR(nid, 'star_lon', NF90_DOUBLE, id_sslon ) + ierr=NF90_DEF_VAR(nid, 'star_lat', NF90_DOUBLE, id_sslat ) + ierr=NF90_DEF_VAR(nid, 'radius', NF90_DOUBLE, id_rad ) + ierr=NF90_DEF_VAR(nid, 'gravity', NF90_DOUBLE, id_grav ) + ierr=NF90_DEF_VAR(nid, 'cp', NF90_DOUBLE, id_cp ) + ierr=NF90_DEF_VAR(nid, 'wavenumber', NF90_DOUBLE, (/iWave/), id_wave ) + ierr=NF90_DEF_VAR(nid, 'wvn_width', NF90_DOUBLE, (/iWave/), id_wavebin) + ierr=NF90_DEF_VAR(nid, 'gauss_point', NF90_DOUBLE, (/iWeight/), id_g ) + ierr=NF90_DEF_VAR(nid, 'latitude', NF90_DOUBLE, (/iLat/), id_lat ) + ierr=NF90_DEF_VAR(nid, 'longitude', NF90_DOUBLE, (/iLon/), id_lon ) + ierr=NF90_DEF_VAR(nid, 'tsurf', NF90_DOUBLE, (/iLon,iLat/), id_tsurf ) + ierr=NF90_DEF_VAR(nid, 'pressure', NF90_DOUBLE, (/iLon,iLat,iLev/), id_plev ) + ierr=NF90_DEF_VAR(nid, 'level', NF90_DOUBLE, (/iLon,iLat,iLev/), id_zlev ) + ierr=NF90_DEF_VAR(nid, 'layer', NF90_DOUBLE, (/iLon,iLat,iLay/), id_zlay ) + ierr=NF90_DEF_VAR(nid, 'dtsw', NF90_DOUBLE, (/iLon,iLat,iLay/), id_dtsw ) + ierr=NF90_DEF_VAR(nid, 'surf_albedo', NF90_DOUBLE, (/iLon,iLat,iWave/), id_salb ) + ierr=NF90_DEF_VAR(nid, 'dtauv', NF90_DOUBLE, (/iLon,iLat,iLay,iWave,iWeight/),id_dtauv ) + ierr=NF90_DEF_VAR(nid, 'ssav', NF90_DOUBLE, (/iLon,iLat,iLay,iWave,iWeight/),id_ssav ) + ierr=NF90_DEF_VAR(nid, 'asymv', NF90_DOUBLE, (/iLon,iLat,iLay,iWave,iWeight/),id_asymv ) +#else + ierr=NF90_DEF_VAR(nid, 'ls', NF90_FLOAT, id_ls ) + ierr=NF90_DEF_VAR(nid, 'dist', NF90_FLOAT, id_dist ) + ierr=NF90_DEF_VAR(nid, 'star_lon', NF90_FLOAT, id_sslon ) + ierr=NF90_DEF_VAR(nid, 'star_lat', NF90_FLOAT, id_sslat ) + ierr=NF90_DEF_VAR(nid, 'radius', NF90_FLOAT, id_rad ) + ierr=NF90_DEF_VAR(nid, 'gravity', NF90_FLOAT, id_grav ) + ierr=NF90_DEF_VAR(nid, 'cp', NF90_FLOAT, id_cp ) + ierr=NF90_DEF_VAR(nid, 'wavenumber', NF90_FLOAT, (/iWave/), id_wave ) + ierr=NF90_DEF_VAR(nid, 'wvn_width', NF90_FLOAT, (/iWave/), id_wavebin) + ierr=NF90_DEF_VAR(nid, 'gauss_point', NF90_FLOAT, (/iWeight/), id_g ) + ierr=NF90_DEF_VAR(nid, 'latitude', NF90_FLOAT, (/iLat/), id_lat ) + ierr=NF90_DEF_VAR(nid, 'longitude', NF90_FLOAT, (/iLon/), id_lon ) + ierr=NF90_DEF_VAR(nid, 'tsurf', NF90_FLOAT, (/iLon,iLat/), id_tsurf ) + ierr=NF90_DEF_VAR(nid, 'pressure', NF90_FLOAT, (/iLon,iLat,iLev/), id_plev ) + ierr=NF90_DEF_VAR(nid, 'level', NF90_FLOAT, (/iLon,iLat,iLev/), id_zlev ) + ierr=NF90_DEF_VAR(nid, 'layer', NF90_FLOAT, (/iLon,iLat,iLay/), id_zlay ) + ierr=NF90_DEF_VAR(nid, 'dtsw', NF90_FLOAT, (/iLon,iLat,iLay/), id_dtsw ) + ierr=NF90_DEF_VAR(nid, 'surf_albedo', NF90_FLOAT, (/iLon,iLat,iWave/), id_salb ) + ierr=NF90_DEF_VAR(nid, 'dtauv', NF90_FLOAT, (/iLon,iLat,iLay,iWave,iWeight/),id_dtauv ) + ierr=NF90_DEF_VAR(nid, 'ssav', NF90_FLOAT, (/iLon,iLat,iLay,iWave,iWeight/),id_ssav ) + ierr=NF90_DEF_VAR(nid, 'asymv', NF90_FLOAT, (/iLon,iLat,iLay,iWave,iWeight/),id_asymv ) +#endif + + ierr = NF90_PUT_ATT(nid, id_ls , 'units', 'rad') + ierr = NF90_PUT_ATT(nid, id_dist , 'units', 'm') + ierr = NF90_PUT_ATT(nid, id_sslon , 'units', 'rad') + ierr = NF90_PUT_ATT(nid, id_sslat , 'units', 'rad') + ierr = NF90_PUT_ATT(nid, id_rad , 'units', 'm') + ierr = NF90_PUT_ATT(nid, id_grav , 'units', 'm/s2') + ierr = NF90_PUT_ATT(nid, id_cp , 'units', 'J/kg/K') + ierr = NF90_PUT_ATT(nid, id_wave , 'units', 'cm-1') + ierr = NF90_PUT_ATT(nid, id_wavebin, 'units', 'cm-1') + ierr = NF90_PUT_ATT(nid, id_g , 'units', '') + ierr = NF90_PUT_ATT(nid, id_lat , 'units', '°') + ierr = NF90_PUT_ATT(nid, id_lon , 'units', '°') + ierr = NF90_PUT_ATT(nid, id_tsurf , 'units', 'K') + ierr = NF90_PUT_ATT(nid, id_plev , 'units', 'Pa') + ierr = NF90_PUT_ATT(nid, id_zlev , 'units', 'm') + ierr = NF90_PUT_ATT(nid, id_zlay , 'units', 'm') + ierr = NF90_PUT_ATT(nid, id_salb , 'units', '') + ierr = NF90_PUT_ATT(nid, id_dtauv , 'units', '') + ierr = NF90_PUT_ATT(nid, id_ssav , 'units', '') + ierr = NF90_PUT_ATT(nid, id_asymv , 'units', '') + ierr = NF90_PUT_ATT(nid, id_dtsw , 'units', 'K/s') + + ierr = NF90_PUT_ATT(nid, id_ls , 'long_name', 'Ls') + ierr = NF90_PUT_ATT(nid, id_dist , 'long_name', 'Stellar distance') + ierr = NF90_PUT_ATT(nid, id_sslon , 'long_name', 'Sub-solar longitude') + ierr = NF90_PUT_ATT(nid, id_sslat , 'long_name', 'Sub-solar latitude') + ierr = NF90_PUT_ATT(nid, id_rad , 'long_name', 'Planet radius') + ierr = NF90_PUT_ATT(nid, id_grav , 'long_name', 'Planet surface gravity') + ierr = NF90_PUT_ATT(nid, id_cp , 'long_name', 'Planet surface gravity') + ierr = NF90_PUT_ATT(nid, id_wave , 'long_name', 'Wavenumber') + ierr = NF90_PUT_ATT(nid, id_wavebin, 'long_name', 'Wavenumber bin width') + ierr = NF90_PUT_ATT(nid, id_g , 'long_name', 'Gauss point weigth') + ierr = NF90_PUT_ATT(nid, id_lat , 'long_name', 'Latitude') + ierr = NF90_PUT_ATT(nid, id_lon , 'long_name', 'Longitude') + ierr = NF90_PUT_ATT(nid, id_tsurf , 'long_name', 'Surface Temperature') + ierr = NF90_PUT_ATT(nid, id_plev , 'long_name', 'Pressure') + ierr = NF90_PUT_ATT(nid, id_zlev , 'long_name', 'Altitude') + ierr = NF90_PUT_ATT(nid, id_zlay , 'long_name', 'Altitude') + ierr = NF90_PUT_ATT(nid, id_salb , 'long_name', 'Surface Albedo') + ierr = NF90_PUT_ATT(nid, id_dtauv , 'long_name', 'Layer opacity') + ierr = NF90_PUT_ATT(nid, id_ssav , 'long_name', 'Single scattering albedo') + ierr = NF90_PUT_ATT(nid, id_asymv , 'long_name', 'Asymmetry parameter') + ierr = NF90_PUT_ATT(nid, id_dtsw , 'long_name', 'Heating Rate') + + ierr=NF90_ENDDEF(nid) + + ierr = NF90_PUT_VAR(nid,id_ls ,ls) + ierr = NF90_PUT_VAR(nid,id_dist ,dist_star) + ierr = NF90_PUT_VAR(nid,id_sslon ,zlss) + ierr = NF90_PUT_VAR(nid,id_sslat ,declin) + ierr = NF90_PUT_VAR(nid,id_rad ,rad) + ierr = NF90_PUT_VAR(nid,id_grav ,g) + ierr = NF90_PUT_VAR(nid,id_cp ,cpp) + ierr = NF90_PUT_VAR(nid,id_wave ,wnov) + ierr = NF90_PUT_VAR(nid,id_wavebin,dwnv) + ierr = NF90_PUT_VAR(nid,id_g ,gweight) + ierr = NF90_PUT_VAR(nid,id_lat ,lat_glo_lonlat(0,:)) + ierr = NF90_PUT_VAR(nid,id_lon ,lon_glo_lonlat(:,1)) + ierr = NF90_PUT_VAR(nid,id_tsurf ,tsurf_glo_lonlat) + ierr = NF90_PUT_VAR(nid,id_plev ,plev_glo_lonlat) + ierr = NF90_PUT_VAR(nid,id_zlev ,zlev_glo_lonlat) + ierr = NF90_PUT_VAR(nid,id_zlay ,zlay_glo_lonlat) + ierr = NF90_PUT_VAR(nid,id_salb ,surf_albedo_glo_lonlat) + ierr = NF90_PUT_VAR(nid,id_dtauv ,htrdr_dtauv_glo_lonlat) + ierr = NF90_PUT_VAR(nid,id_ssav ,htrdr_ssav_glo_lonlat) + ierr = NF90_PUT_VAR(nid,id_asymv ,htrdr_asymv_glo_lonlat) + ierr = NF90_PUT_VAR(nid,id_dtsw ,zdtsw_glo_lonlat) + + ierr = NF90_CLOSE (nid) + endif + + end subroutine hr_write_oprop_nc + + subroutine deallocate_hrcorr + deallocate(tau_MC_2S,solLon,latCorr) + end subroutine deallocate_hrcorr + + SUBROUTINE hr_optcv(PQMO,NLAY,ZLEV,PLEV,TMID,PMID, & + DTAUV,TAUV,TAUCUMV,WBARV,COSBV,TAURAY,TAUGSURF,SEASHAZEFACT,& + DIAG_OPTH,DIAG_OPTT,CDCOLUMN) + + use radinc_h + use radcommon_h, only: gasv,gasv_recomb,tlimit,Cmk,gzlat_ig, & + tgasref,pfgasref,wnov,scalep,indv + use gases_h + use datafile_mod, only: haze_opt_file + use comcstfi_mod, only: pi,r + use callkeys_mod, only: continuum,graybody,callgasvis,corrk_recombin, & + callclouds,callmufi,seashaze,uncoupl_optic_haze,& + opt4clouds,FHVIS,FCVIS + use tracer_h, only: nmicro,nice,ices_indx + + implicit none + + !================================================================== + ! + ! Purpose + ! ------- + ! Calculates shortwave optical constants at each level. + ! + ! Authors + ! ------- + ! Adapted from the NASA Ames code by R. Wordsworth (2009) + ! + ! Modified + ! -------- + ! J. Vatant d'Ollone (2016-17) + ! --> Clean and adaptation to Titan + ! B. de Batz de Trenquelléon (2022-2023) + ! --> Clean and correction + ! --> New optics added for Titan's clouds + ! A. Arfaux (2026) + ! --> add the half-layer shift + ! + !================================================================== + ! + ! THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL + ! IT CALCULATES FOR EACH LAYER, FOR EACH SPECTRAL INTERVAL IN THE VISUAL + ! LAYER: WBAR, DTAU, COSBAR + ! LEVEL: TAU + ! + ! TAUV(L,NW,NG) is the cumulative optical depth at the top of radiation code + ! layer L. NW is spectral wavelength interval, ng the Gauss point index. + ! + ! TLEV(L) - Temperature at the layer boundary + ! PLEV(L) - Pressure at the layer boundary (i.e. level) + ! GASV(NT,NPS,NW,NG) - Visible k-coefficients + ! + !------------------------------------------------------------------- + + + !========================================================== + ! Input/Output + !========================================================== + REAL*8, INTENT(IN) :: PQMO(nlay,nmicro) ! Tracers for microphysics optics (X/m2). + INTEGER, INTENT(IN) :: NLAY ! Number of pressure layers (for pqmo) + REAL*8, INTENT(IN) :: ZLEV(NLAY+1) + REAL*8, INTENT(IN) :: PLEV(L_LEVELS) + REAL*8, INTENT(IN) :: TMID(L_LEVELS), PMID(L_LEVELS) + REAL*8, INTENT(IN) :: TAURAY(L_NSPECTV) + REAL*8, INTENT(IN) :: SEASHAZEFACT(L_LEVELS) + INTEGER, INTENT(IN) :: CDCOLUMN + + REAL*8, INTENT(OUT) :: DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS) + REAL*8, INTENT(OUT) :: TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS) + REAL*8, INTENT(OUT) :: TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS) + REAL*8, INTENT(OUT) :: COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS) + REAL*8, INTENT(OUT) :: WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS) + REAL*8, INTENT(OUT) :: TAUGSURF(L_NSPECTV,L_NGAUSS-1) + REAL*8, INTENT(OUT) :: DIAG_OPTH(L_LEVELS,L_NSPECTV,6) + REAL*8, INTENT(OUT) :: DIAG_OPTT(L_LEVELS,L_NSPECTV,6) + ! ========================================================== + + real*8 DTAUKV(L_LEVELS,L_NSPECTV,L_NGAUSS) + + ! Titan customisation + ! J. Vatant d'Ollone (2016) + real*8 DHAZE_T(L_LEVELS,L_NSPECTV) + real*8 DHAZES_T(L_LEVELS,L_NSPECTV) + real*8 SSA_T(L_LEVELS,L_NSPECTV) + real*8 ASF_T(L_LEVELS,L_NSPECTV) + ! ========================== + + integer L, NW, NG, K, LK, IAER + integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS) + real*8 ANS, TAUGAS + real*8 TRAY(L_LEVELS,L_NSPECTV) + real*8 DPR(L_LEVELS), U(L_LEVELS) + real*8 LCOEF(4), LKCOEF(L_LEVELS,4) + + real*8 DCONT + real*8 DRAYAER + double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc + double precision p_cross + + real*8 KCOEF(4) + + ! temporary variable to reduce memory access time to gasv + real*8 tmpk(2,2) + + ! temporary variables for multiple aerosol calculation + real*8 atemp(L_NLAYRAD,L_NSPECTV) + real*8 btemp(L_NLAYRAD,L_NSPECTV) + real*8 ctemp(L_NLAYRAD,L_NSPECTV) + + ! variables for k in units m^-1 + real*8 dz(L_LEVELS) + + integer igas, jgas, ilay + + integer interm + + ! Variables for haze optics + character(len=200) file_path + logical file_ok + integer dumch + real*8 dumwvl + integer ilev_cutoff ! pressure level index of cutoff + real*8 corr_haze + real*8 corr_alb + + ! Variables for new optics + integer iq, iw, FTYPE, CTYPE + real*8 m0as,m0af,m0ccn,m3as,m3af,m3ccn,m3cld + real*8 dtauaer_s,dtauaer_f,dtau_ccn,dtau_cld + real*8,save :: rhoaer_s(L_NSPECTV),ssa_s(L_NSPECTV),asf_s(L_NSPECTV) + real*8,save :: rhoaer_f(L_NSPECTV),ssa_f(L_NSPECTV),asf_f(L_NSPECTV) + real*8,save :: ssa_ccn(L_NSPECTV),asf_ccn(L_NSPECTV) + real*8,save :: ssa_cld(L_NSPECTV),asf_cld(L_NSPECTV) + !$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,ssa_cld,asf_s,asf_f,asf_cld) + + logical,save :: firstcall=.true. + !$OMP THREADPRIVATE(firstcall) + + !! AS: to save time in computing continuum (see bilinearbig) + IF (.not.ALLOCATED(indv)) THEN + ALLOCATE(indv(L_NSPECTV,ngasmx,ngasmx)) + indv = -9999 ! this initial value means "to be calculated" + ENDIF + + ! Some initialisation beacause there's a pb with disr_haze at the limits (nw=1) + ! I should check this - For now we set vars to zero : better than nans - JVO 2017 + DHAZE_T(:,:) = 0.0 + SSA_T(:,:) = 0.0 + ASF_T(:,:) = 0.0 + + ! Load tabulated haze optical properties if needed. + ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + IF (firstcall .AND. callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN + OPEN(12,file=TRIM(haze_opt_file),form='formatted') ! The file has been inquired in physiq_mod firstcall + READ(12,*) ! dummy header + DO NW=1,L_NSPECTI + READ(12,*) ! there's IR 1st + ENDDO + DO NW=1,L_NSPECTV + READ(12,*) dumch, dumwvl, rhoaer_f(nw), ssa_f(nw), asf_f(nw), rhoaer_s(nw), ssa_s(nw), asf_s(nw) + ENDDO + CLOSE(12) + ENDIF + + !======================================================================= + ! Determine the total gas opacity throughout the column, for each + ! spectral interval, NW, and each Gauss point, NG. + ! Calculate the continuum opacities, i.e., those that do not depend on + ! NG, the Gauss index. + + taugsurf(:,:) = 0.0 + dpr(:) = 0.0 + lkcoef(:,:) = 0.0 + + ! Level of cutoff + !~~~~~~~~~~~~~~~~ + ilev_cutoff = 26 + + do K=2,L_LEVELS + ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed) + DPR(k) = PLEV(K)-PLEV(K-1) + + ! if we have continuum opacities, we need dz + dz(k) = dpr(k)*R*TMID(K)/(gzlat_ig(ilay)*PMID(K)) + U(k) = Cmk(ilay)*DPR(k) ! only Cmk line in optcv.F + + call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K)) + + do LK=1,4 + LKCOEF(K,LK) = LCOEF(LK) + end do + end do ! L_LEVELS + + ! Rayleigh scattering + do NW=1,L_NSPECTV + TRAY(1:4,NW) = 1.d-30 + do K=5,L_LEVELS + TRAY(K,NW) = TAURAY(NW) * DPR(K) + end do ! L_LEVELS + end do + + DIAG_OPTH(:,:,:) = 0.D0 + DIAG_OPTT(:,:,:) = 0.D0 + + do NW=1,L_NSPECTV + ! We ignore K=1... + do K=2,L_LEVELS + ! int. arithmetic => gives the gcm layer index (reversed) + ilay = L_NLAYRAD+1 - k/2 + + ! Optics coupled with the microphysics : + IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN + + !========================================================================================== + ! Old optics (must have callclouds = .false.): + !========================================================================================== + IF (.NOT. opt4clouds) THEN + m3as = pqmo(ilay,2) / 2.0 + m3af = pqmo(ilay,4) / 2.0 + ! Cut-off (here for p = 2.7e3Pa / alt = 70km) + IF (ilay .lt. ilev_cutoff) THEN + m3as = pqmo(ilev_cutoff,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff)) + m3af = pqmo(ilev_cutoff,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff)) + ENDIF + + dtauaer_s = m3as*rhoaer_s(nw) + dtauaer_f = m3af*rhoaer_f(nw) + + !========================================================================================== + ! New optics : + !========================================================================================== + ELSE + iw = (L_NSPECTV + 1) - NW ! Visible first and return + !----------------------------- + ! HAZE (Spherical + Fractal) : + !----------------------------- + FTYPE = 1 + + ! Spherical aerosols : + !--------------------- + CTYPE = 5 + m0as = pqmo(ilay,1) / 2.0 + m3as = pqmo(ilay,2) / 2.0 + ! If not callclouds : must have a cut-off (here for p = 2.7e3Pa / alt = 70km) + IF (.NOT. callclouds) THEN + IF (ilay .lt. ilev_cutoff) THEN + m0as = pqmo(ilev_cutoff,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff)) + m3as = pqmo(ilev_cutoff,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff)) + ENDIF + ENDIF + call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0as,m3as,iw,dtauaer_s,ssa_s(nw),asf_s(nw)) + + ! Fractal aerosols : + !------------------- + CTYPE = FTYPE + m0af = pqmo(ilay,3) / 2.0 + m3af = pqmo(ilay,4) / 2.0 + ! If not callclouds : must have a cut-off (here for p = 2.7e3Pa / alt = 70km) + IF (.NOT. callclouds) THEN + IF (ilay .lt. ilev_cutoff) THEN + m0af = pqmo(ilev_cutoff,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff)) + m3af = pqmo(ilev_cutoff,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(ilev_cutoff+1)-zlev(ilev_cutoff)) + ENDIF + ENDIF + call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0af,m3af,iw,dtauaer_f,ssa_f(nw),asf_f(nw)) + ENDIF + + ! Tuning of optical properties for haze : + !dtauaer_s = dtauaer_s * (FHVIS * (1-ssa_s(nw)) + ssa_s(nw)) + !ssa_s(nw) = ssa_s(nw) / (FHVIS * (1-ssa_s(nw)) + ssa_s(nw)) + dtauaer_f = dtauaer_f * (FHVIS * (1-ssa_f(nw)) + ssa_f(nw)) + ssa_f(nw) = ssa_f(nw) / (FHVIS * (1-ssa_f(nw)) + ssa_f(nw)) + + ! Total of Haze opacity (dtau), SSA (w) and ASF (COS) : + DHAZE_T(k,nw) = dtauaer_s + dtauaer_f + IF (dtauaer_s + dtauaer_f .GT. 1.D-30) THEN + SSA_T(k,nw) = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) ) / ( dtauaer_s+dtauaer_f ) + ASF_T(k,nw) = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) ) & + / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f ) + ELSE + DHAZE_T(k,nw) = 0.D0 + SSA_T(k,nw) = 1.0 + ASF_T(k,nw) = 1.0 + ENDIF + + ! Opacity and albedo adjustement below cutoff + !---------------------------------------------- + corr_haze=0.85-0.15*TANH((PMID(k)*100.-2500.)/250.) + corr_alb =0.85-0.15*TANH((PMID(k)*100.-4500.)/1000.) + IF (.NOT. callclouds) THEN + IF (ilay .lt. ilev_cutoff) THEN + DHAZE_T(k,nw) = DHAZE_T(k,nw) * corr_haze + SSA_T(k,nw) = SSA_T(k,nw) * corr_alb + ENDIF + ENDIF + + ! Diagnostics for the haze : + DIAG_OPTH(k,nw,1) = DHAZE_T(k,nw) ! dtau + DIAG_OPTH(k,nw,2) = SSA_T(k,nw) ! wbar + DIAG_OPTH(k,nw,3) = ASF_T(k,nw) ! gbar + + !--------------------- + ! CLOUDS (Spherical) : + !--------------------- + IF (callclouds) THEN + CTYPE = 0 + m0ccn = pqmo(ilay,5) / 2.0 + m3ccn = pqmo(ilay,6) / 2.0 + m3cld = pqmo(ilay,6) / 2.0 + + ! Clear / Dark column method : + !----------------------------- + + ! CCN's SSA : + call get_haze_and_cloud_opacity(FTYPE,FTYPE,m0ccn,m3ccn,iw,dtau_ccn,ssa_ccn(nw),asf_ccn(nw)) + + ! Clear column (CCN + minor ices): + IF (CDCOLUMN == 0) THEN + DO iq = 2, nice + m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0) + ENDDO + call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw)) + + ! Dark column (CCN + CH4 ice + minor ices): + ELSEIF (CDCOLUMN == 1) THEN + DO iq = 1, nice + m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0) + ENDDO + call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw)) + + ELSE + WRITE(*,*) 'WARNING OPTCV.F90 : CDCOLUMN MUST BE 0 OR 1' + WRITE(*,*) 'WE USE DARK COLUMN ...' + DO iq = 1, nice + m3cld = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0) + ENDDO + call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw)) + ENDIF + + ! For small dropplets, opacity of nucleus dominates + IF ((m3ccn + m3cld) .le. tiny(m3ccn)) THEN ! No cloud + dtau_cld = 0. + ssa_cld(nw) = 0. + ELSE + dtau_cld = (dtau_cld*m3ccn + dtau_cld*m3cld) / (m3ccn + m3cld) + ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3cld) / (m3ccn + m3cld) + ENDIF + + ! Tuning of optical properties for clouds : + dtau_cld = dtau_cld * (FCVIS * (1-ssa_cld(nw)) + ssa_cld(nw)) + ssa_cld(nw) = ssa_cld(nw) / (FCVIS * (1-ssa_cld(nw)) + ssa_cld(nw)) + + ! Total of Haze + Clouds opacity (dtau), SSA (w) and ASF (COS) : + DHAZE_T(k,nw) = dtauaer_s + dtauaer_f + dtau_cld + IF (DHAZE_T(k,nw) .GT. 1.D-30) THEN + SSA_T(k,nw) = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) + dtau_cld*ssa_cld(nw) ) / ( dtauaer_s+dtauaer_f+dtau_cld ) + ASF_T(k,nw) = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) + dtau_cld*ssa_cld(nw)*asf_cld(nw) ) & + / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f + ssa_cld(nw)*dtau_cld ) + ELSE + DHAZE_T(k,nw) = 0.D0 + SSA_T(k,nw) = 1.0 + ASF_T(k,nw) = 1.0 + ENDIF + + ! Diagnostics for clouds : + DIAG_OPTT(k,nw,1) = DHAZE_T(k,nw) ! dtau + DIAG_OPTT(k,nw,2) = SSA_T(k,nw) ! wbar + DIAG_OPTT(k,nw,3) = ASF_T(k,nw) ! gbar + + ELSE + ! Diagnostics for clouds : + DIAG_OPTT(k,nw,1) = 0.D0 ! dtau + DIAG_OPTT(k,nw,2) = 1.0 ! wbar + DIAG_OPTT(k,nw,3) = 1.0 ! gbar + ENDIF + + ! Optics and microphysics no coupled : + ELSE + ! Call fixed vertical haze profile of extinction - same for all columns + call disr_haze(dz(k),plev(k),wnov(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw)) + if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k) + ! Diagnostics for the haze : + DIAG_OPTH(k,nw,1) = DHAZE_T(k,nw) ! dtau + DIAG_OPTH(k,nw,2) = SSA_T(k,nw) ! wbar + DIAG_OPTH(k,nw,3) = ASF_T(k,nw) ! gbar + ! Diagnostics for clouds : + DIAG_OPTT(k,nw,1) = 0.D0 ! dtau + DIAG_OPTT(k,nw,2) = 1.0 ! wbar + DIAG_OPTT(k,nw,3) = 1.0 ! gbar + ENDIF ! ENDIF callmufi + + !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR + ! but visible does not handle very well diffusion in first layer. + ! The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise) + ! in the 4 first semilayers in optcv, but not optci. + ! This solves random variations of the sw heating at the model top. + if (k<5) DHAZE_T(K,:) = 0.0 + + DRAYAER = TRAY(K,NW) + ! DRAYAER is Tau RAYleigh scattering, plus AERosol opacity + DRAYAER = DRAYAER + DHAZE_T(K,NW) ! Titan's aerosol + + DCONT = 0.0 ! continuum absorption + + if(continuum.and.(.not.graybody).and.callgasvis)then + ! include continua if necessary + wn_cont = dble(wnov(nw)) + T_cont = dble(TMID(k)) + do igas=1,ngasmx + + p_cont = dble(PMID(k)*scalep*gfrac(igas,ilay)) + + dtemp=0.0 + if(igas.eq.igas_N2)then + + interm = indv(nw,igas,igas) + ! call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm) + indv(nw,igas,igas) = interm + ! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible + + elseif(igas.eq.igas_H2)then + + ! first do self-induced absorption + interm = indv(nw,igas,igas) + call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm) + indv(nw,igas,igas) = interm + + ! then cross-interactions with other gases + do jgas=1,ngasmx + p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay)) + dtempc = 0.0 + if(jgas.eq.igas_N2)then + interm = indv(nw,igas,jgas) + call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) + indv(nw,igas,jgas) = interm + ! should be irrelevant in the visible + endif + dtemp = dtemp + dtempc + enddo + + elseif(igas.eq.igas_CH4)then + + ! first do self-induced absorption + interm = indv(nw,igas,igas) + call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm) + indv(nw,igas,igas) = interm + + ! then cross-interactions with other gases + do jgas=1,ngasmx + p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay)) + dtempc = 0.0 + if(jgas.eq.igas_N2)then + interm = indv(nw,igas,jgas) + call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm) + indv(nw,igas,jgas) = interm + endif + dtemp = dtemp + dtempc + enddo + + endif + + DCONT = DCONT + dtemp + + enddo + + DCONT = DCONT*dz(k) + + endif + + do ng=1,L_NGAUSS-1 + + ! Now compute TAUGAS + + ! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically + ! the execution time of optci/v -> ~ factor 2 on the whole radiative + ! transfer on the tested simulations ! + + if (corrk_recombin) then + tmpk = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NW,NG) + else + tmpk = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG) + endif + + KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG) + KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG) + KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG) + KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG) + + ! Interpolate the gaseous k-coefficients to the requested T,P values + + ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) + & + LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4) + + + TAUGAS = U(k)*ANS + + TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT + DTAUKV(K,nw,ng) = TAUGAS & + + DRAYAER & ! DRAYAER includes all scattering contributions + + DCONT ! For parameterized continuum aborption + end do + + ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS), + ! which holds continuum opacity only + + NG = L_NGAUSS + DTAUKV(K,nw,ng) = DRAYAER + DCONT ! Scattering + parameterized continuum absorption, including Titan's haze + + DIAG_OPTH(K,nw,4) = DRAYAER + DIAG_OPTH(K,nw,5) = TAUGAS + DIAG_OPTH(K,nw,6) = DCONT + DIAG_OPTT(K,nw,4) = DRAYAER + DIAG_OPTT(K,nw,5) = TAUGAS + DIAG_OPTT(K,nw,6) = DCONT + + end do ! K = L_LEVELS + end do ! nw = L_NSPECTV + + !======================================================================= + ! Now the full treatment for the layers, where besides the opacity + ! we need to calculate the scattering albedo and asymmetry factors + ! ====================================================================== + + ! Haze scattering + !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR + ! but not in the visible + ! The dhaze_s is thus set to 0 in the 4 first semilayers in optcv, but not optci. + ! This solves random variations of the sw heating at the model top. + DO NW=1,L_NSPECTV + DHAZES_T(1:4,NW) = 0.d0 + DO K=5,L_LEVELS + DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW) ! effect of scattering albedo on haze + ENDDO + ENDDO + + ! NW spectral loop + DO NW=1,L_NSPECTV + ! L vertical loop + DO L=1,L_NLAYRAD-1 + K = 2*L + atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW) + ASF_T(K+1,NW)*DHAZES_T(K+1,NW) + btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW) + ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ? + btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW) + COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW) + END DO + + ! Last level + L = L_NLAYRAD + K = 2*L + atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW) + btemp(L,NW) = DHAZES_T(K,NW) + ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW) ! JVO 2017 : does this 0.999 is really meaningful ? + btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW) + END DO + + ! NG Gauss loop + DO NG=1,L_NGAUSS + ! NW spectral loop + DO NW=1,L_NSPECTV + ! L vertical loop + DO L=1,L_NLAYRAD-1 + K = 2*L + DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG) + WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng) + END DO + + ! Last level + L = L_NLAYRAD + K = 2*L + DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG) + END DO + END DO + + TAUCUMV(:,:,:) = DTAUKV(:,:,:) + DO L=1,L_NLAYRAD + TAUV(L,:,:)=TAUCUMV(2*L,:,:) + END DO + TAUV(L,:,:)=TAUCUMV(2*L_NLAYRAD+1,:,:) + + if(firstcall) firstcall = .false. + end subroutine hr_optcv + +END MODULE hrcorr_mod Index: trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90 (revision 4152) +++ trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90 (revision 4153) @@ -297,4 +297,14 @@ write(*,*) " corrk = ",corrk + write(*,*) "call the 3D correction of the sw heating rate ?" + corrhtrdr=.false. ! default value + call getin_p("corrhtrdr",corrhtrdr) + write(*,*) " corrhtrdr = ",corrhtrdr + + write(*,*) "update the correction factor every week ?" + corrhtrdr=.false. ! default value + call getin_p("updatecorrhtrdr",updatecorrhtrdr) + write(*,*) " updatecorrhtrdr = ",updatecorrhtrdr + if (corrk) then ! default path is set in datadir Index: trunk/LMDZ.TITAN/libf/phytitan/optci.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/optci.F90 (revision 4152) +++ trunk/LMDZ.TITAN/libf/phytitan/optci.F90 (revision 4153) @@ -400,5 +400,5 @@ enddo - elseif(igas.eq.igas_CH4)then + elseif(igas.eq.igas_CH4)then ! first do self-induced absorption Index: trunk/LMDZ.TITAN/libf/phytitan/optcv.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/optcv.F90 (revision 4152) +++ trunk/LMDZ.TITAN/libf/phytitan/optcv.F90 (revision 4153) @@ -133,5 +133,4 @@ logical,save :: firstcall=.true. !$OMP THREADPRIVATE(firstcall) - !! AS: to save time in computing continuum (see bilinearbig) @@ -537,8 +536,8 @@ DO L=1,L_NLAYRAD-1 K = 2*L+1 - atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW) + ASF_T(K+1,NW)*DHAZES_T(K+1,NW) + atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW) + ASF_T(K+1,NW)*DHAZES_T(K+1,NW) btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW) - ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ? - btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW) + ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ? + btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW) COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW) END DO @@ -568,5 +567,5 @@ L = L_NLAYRAD K = 2*L+1 - DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG) END DO Index: trunk/LMDZ.TITAN/libf/phytitan/phyredem.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/phyredem.F90 (revision 4152) +++ trunk/LMDZ.TITAN/libf/phytitan/phyredem.F90 (revision 4153) @@ -68,5 +68,5 @@ tab_cntrl(10) = daysec ! length of a sol (s) ~88775 - tab_cntrl(11) = phystep ! time step in the physics + tab_cntrl(11) = phystep ! time step in the physics tab_cntrl(12) = 0. tab_cntrl(13) = 0. @@ -93,5 +93,5 @@ tab_cntrl(34) = dtemisice(2) ! time scale for snow metamorphism (south) - tab_cntrl(28) = 0. + tab_cntrl(28) = 0. tab_cntrl(29) = 0. tab_cntrl(30) = 0. @@ -186,5 +186,5 @@ ! Methane tank depth call put_field("tankCH4","Depth of methane tank",tankCH4) - + ! Tracers if (nq>0) then Index: trunk/LMDZ.TITAN/libf/phytitan/phys_state_var_mod.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/phys_state_var_mod.F90 (revision 4152) +++ trunk/LMDZ.TITAN/libf/phytitan/phys_state_var_mod.F90 (revision 4153) @@ -77,4 +77,8 @@ real,dimension(:,:,:,:),allocatable,save :: zpopttv ! VI optical properties [haze+clouds] within narrowbands for diags (dtau,tau,k,wbar,gbar,drayaer,taugaz,dcont). !$OMP THREADPRIVATE(zpopthi,zpopthv,zpoptti,zpopttv) + real,dimension(:,:,:,:),allocatable,save :: htrdr_dtauv ! opacity of layers for htrdr coupling + real,dimension(:,:,:,:),allocatable,save :: htrdr_ssav ! single scattering albedo for htrdr coupling + real,dimension(:,:,:,:),allocatable,save :: htrdr_asymv ! essymetrie parameter for htrdr coupling +!$OMP THREADPRIVATE(htrdr_dtauv,htrdr_ssav,htrdr_asymv) real,dimension(:,:),allocatable,save :: u_ref ! Reference profile for the zonal wind nudging (m.s-1). @@ -209,4 +213,7 @@ DEALLOCATE(zpoptti) DEALLOCATE(zpopttv) + DEALLOCATE(htrdr_dtauv) + DEALLOCATE(htrdr_ssav) + DEALLOCATE(htrdr_asymv) DEALLOCATE(u_ref) DEALLOCATE(dycchi) Index: trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90 =================================================================== --- trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90 (revision 4152) +++ trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90 (revision 4153) @@ -59,4 +59,5 @@ #endif use muphy_diag + use hrcorr_mod implicit none @@ -248,4 +249,5 @@ real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine. real zdtlc(ngrid,nlayer) ! Condensation heating rate. + real refCorr, time ! for hrcorr_mod routines. ! For Surface Tracers : (kg/m2/s) @@ -464,5 +466,5 @@ call setspv ! Basic visible properties. call sugas_corrk ! Set up gaseous absorption properties. - + OLR_nu(:,:) = 0.D0 OSR_nu(:,:) = 0.D0 @@ -483,4 +485,8 @@ endif ENDIF + + if (corrhtrdr) then + call ini_hrcorr + endif !#endif @@ -820,10 +826,10 @@ ! omega in Pa/s do l=1,nlayer-1 - omega(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1)) - enddo - omega(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0 - do l=1,nlayer - omega(1:ngrid,l)=g*omega(1:ngrid,l)/cell_area(1:ngrid) - enddo + omega(1:ngrid,l)=0.5*(flxw(1:ngrid,l)+flxw(1:ngrid,l+1)) + enddo + omega(1:ngrid,nlayer)=0.5*flxw(1:ngrid,nlayer) ! since flxw(nlayer+1)=0 + do l=1,nlayer + omega(1:ngrid,l)=g*omega(1:ngrid,l)/cell_area(1:ngrid) + enddo !--------------------------------- @@ -876,5 +882,6 @@ ! standard callcorrk call 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, & zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw, & @@ -882,5 +889,5 @@ fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu, & int_dtaui,int_dtauv,zpopthi,zpopthv,zpoptti,zpopttv,& - lastcall) + lastcall,zlss,zls) ! Radiative flux from the sky absorbed by the surface (W.m-2). @@ -895,6 +902,17 @@ ! Net atmospheric radiative heating rate (K.s-1) + if(corrhtrdr) then + call hr_average(ngrid, nlayer, pplay, zzlev, zdtsw, refCorr) + do ig=1,ngrid + time = modulo(ptime+longitude(ig)/(2.*pi), 1.) + if(is_master) write(*,*) longitude(ig)*180./pi, time + call apply_hrcorr(ig, nlayer, zls*180./pi, time , zdtsw(ig,:), refCorr) + enddo + if (lastcall) then + call deallocate_hrcorr + endif + endif dtrad(:,:)=zdtsw(:,:)+zdtlw(:,:) - + elseif(newtonian)then @@ -973,5 +991,5 @@ ! JL12 the following if test is temporarily there to allow us to compare the old vdifc with turbdiff. if (UseTurbDiff) then - + call turbdiff(ngrid,nlayer,nq, & ptimestep,capcal,lwrite, &