module sfluxi_PLUTO_mod implicit none contains SUBROUTINE SFLUXI_PLUTO(PLEV,TLEV,DTAUI,TAUCUMI,UBARI,RSFI,WNOI, * DWNI,COSBI,WBARI,NFLUXTOPI,NFLUXTOPI_nu, * FMNETI,fmneti_nu,fluxupi,fluxdni,fluxupi_nu, * FZEROI,TAUGSURF) use radinc_h, only: NTfac, NTstart, L_LEVELS, L_NSPECTI, L_NGAUSS use radinc_h, only: L_NLAYRAD, L_NLEVRAD use radcommon_h, only: planckir, tlimit,sigma, gweight use comcstfi_mod, only: pi use gfluxi_old_mod, only: gfluxi_old implicit none integer NLEVRAD, L, NW, NG, NTS, NTT real*8 TLEV(L_LEVELS), PLEV(L_LEVELS) real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS) real*8 FMNETI(L_NLAYRAD) real*8 FMNETI_NU(L_NLAYRAD,L_NSPECTI) real*8 WNOI(L_NSPECTI), DWNI(L_NSPECTI) real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS) real*8 FMUPI(L_NLEVRAD), FMDI(L_NLEVRAD) real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS) real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS) real*8 NFLUXTOPI real*8 NFLUXTOPI_nu(L_NSPECTI) real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) real*8 FTOPUP real*8 UBARI, RSFI, TSURF, BSURF, TTOP, BTOP, TAUTOP real*8 PLANCK, PLTOP real*8 fluxupi(L_NLAYRAD), fluxdni(L_NLAYRAD) real*8 FZEROI(L_NSPECTI) real*8 taugsurf(L_NSPECTI,L_NGAUSS-1), fzero real*8 BSURFtest ! by RW for test real*8 fup_tmp(L_NSPECTI),fdn_tmp(L_NSPECTI) real*8 XX,YY ! calcul intermediate C======================================================================C NLEVRAD = L_NLEVRAD C ZERO THE NET FLUXES NFLUXTOPI = 0.0 DO NW=1,L_NSPECTI NFLUXTOPI_nu(NW) = 0.0 DO L=1,L_NLAYRAD FLUXUPI_nu(L,NW) = 0.0 FMNETI_nu(L,NW) = 0.0 fup_tmp(nw)=0.0 fdn_tmp(nw)=0.0 END DO END DO DO L=1,L_NLAYRAD FMNETI(L) = 0.0 FLUXUPI(L) = 0.0 FLUXDNI(L) = 0.0 END DO C WE NOW ENTER A MAJOR LOOP OVER SPECTRAL INTERVALS IN THE INFRARED C TO CALCULATE THE NET FLUX IN EACH SPECTRAL INTERVAL TTOP = TLEV(2) TSURF = TLEV(L_LEVELS) !PRINT*, 'TB18 NTstart: ',NTstart,TSURF NTS = int(TSURF*10.0D0)-NTstart + 1 !PRINT*, 'TB18 NTS: ',NTS !PRINT*, 'TB18 TLEV sfluxi: ',TLEV NTT = int(TTOP *10.0D0)-NTstart +1 !PRINT*, 'TB18 NTT: ',NTT ! if temperatures superior to 50K ! NTS = TSURF*10.0D0-499 ! NTT = TTOP *10.0D0-499 ! analyse c write(*,*)'TTOP',TTOP c write(*,*)'TSURF',TSURF c write(*,*)'NTS',NTS c write(*,*)'NTT',NTT c DO L=1,L_LEVELS c write(*,*)'TLEV',TLEV c ENDDO BSURFtest=0.0 DO 501 NW=1,L_NSPECTI C SURFACE EMISSIONS - INDEPENDENT OF GAUSS POINTS BSURF = (1.-RSFI)*PLANCKIR(NW,NTS) ! interpolate for accuracy?? PLTOP = PLANCKIR(NW,NTT) !BSURFtest=BSURFtest+BSURF*DWNI(NW) !if(NW.eq.L_NSPECTI)then ! print*,'eps*sigma*T^4',5.67e-8*(1.-RSFI)*TSURF**4 ! print*,'BSURF in sfluxi=',pi*BSURFtest !endif C If FZEROI(NW) = 1, then the k-coefficients are zero - skip to the C special Gauss point at the end. FZERO = FZEROI(NW) IF(FZERO.ge.0.99) goto 40 DO NG=1,L_NGAUSS-1 if(TAUGSURF(NW,NG).lt. TLIMIT) then fzero = fzero + (1.0-FZEROI(NW))*GWEIGHT(NG) goto 30 end if C SET UP THE UPPER AND LOWER BOUNDARY CONDITIONS ON THE IR C CALCULATE THE DOWNWELLING RADIATION AT THE TOP OF THE MODEL C OR THE TOP LAYER WILL COOL TO SPACE UNPHYSICALLY TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2)) BTOP = (1.0-EXP(-TAUTOP/UBARI))*PLTOP C WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM C CALL A SUBROUTINE THAT SOLVES FOR THE FLUX TERMS C WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER CALL GFLUXI_OLD(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG), * TAUCUMI(1,NW,NG), * WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP, * BSURF,FTOPUP,FMUPI,FMDI) C NOW CALCULATE THE CUMULATIVE IR NET FLUX XX= FTOPUP*DWNI(NW)*GWEIGHT(NG)*(1.0-FZEROI(NW)) NFLUXTOPI = NFLUXTOPI+ XX NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW) + XX !same by spectral band. DO L=1,L_NLEVRAD-1 C CORRECT FOR THE WAVENUMBER INTERVALS YY = DWNI(NW)*GWEIGHT(NG)*(1.0-FZEROI(NW)) FMNETI(L) = FMNETI(L)+(FMUPI(L)-FMDI(L))* YY FMNETI_nu(L,NW) = FMNETI_nu(L,NW)+(FMUPI(L)-FMDI(L))* YY !same by spectral band FLUXDNI(L) = FLUXDNI(L) + FMDI(L)* YY FLUXUPI(L) = FLUXUPI(L) + FMUPI(L)* YY FLUXUPI_nu(L,NW) = FLUXUPI_nu(L,NW) + FMUPI(L)* YY !same by spectral band. END DO !fup_tmp(nw)=fup_tmp(nw)+FMUPI(L_NLEVRAD-1)*DWNI(NW)*GWEIGHT(NG) !fdn_tmp(nw)=fdn_tmp(nw)+FMDI(L_NLEVRAD-1)*DWNI(NW)*GWEIGHT(NG) !fup_tmp(nw)=fup_tmp(nw)+FMUPI(1)*DWNI(NW)*GWEIGHT(NG) !fdn_tmp(nw)=fdn_tmp(nw)+FMDI(1)*DWNI(NW)*GWEIGHT(NG) 30 CONTINUE END DO !End NGAUSS LOOP ! print*,'-----------------------------------' !print*,'FMDI(',nw,')=',FMDI(L_NLEVRAD-1) !print*,'FMUPI(',nw,')=',FMUPI(L_NLEVRAD-1) !print*,'DWNI(',nw,')=',DWNI(nw) ! print*,'-----------------------------------' 40 CONTINUE C SPECIAL 17th Gauss point ! print*,'fzero for ng=17',fzero NG = L_NGAUSS TAUTOP = DTAUI(1,NW,NG)*PLEV(2)/(PLEV(4)-PLEV(2)) BTOP = (1.0-EXP(-TAUTOP/UBARI))*PLTOP C WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM C CALL A SUBROUTINE THAT SOLVES FOR THE FLUX TERMS C WITHIN EACH INTERVAL AT THE MIDPOINT WAVENUMBER CALL GFLUXI_OLD(NLEVRAD,TLEV,NW,DWNI(NW),DTAUI(1,NW,NG), * TAUCUMI(1,NW,NG), * WBARI(1,NW,NG),COSBI(1,NW,NG),UBARI,RSFI,BTOP, * BSURF,FTOPUP,FMUPI,FMDI) C NOW CALCULATE THE CUMULATIVE IR NET FLUX NFLUXTOPI = NFLUXTOPI+FTOPUP*DWNI(NW)*FZERO NFLUXTOPI_nu(NW) = NFLUXTOPI_nu(NW)+FTOPUP*DWNI(NW)*FZERO !same by spectral band. DO L=1,L_NLEVRAD-1 C CORRECT FOR THE WAVENUMBER INTERVALS FLUXDNI(L) = FLUXDNI(L) + FMDI(L)*DWNI(NW)*FZERO XX = (FMUPI(L)-FMDI(L))*DWNI(NW)*FZERO FMNETI(L) = FMNETI(L)+ XX FMNETI_nu(L,NW) = FMNETI_nu(L,NW)+XX !same by spectral band. XX= FMUPI(L)*DWNI(NW)*FZERO FLUXUPI(L) = FLUXUPI(L) + XX FLUXUPI_nu(L,NW) = FLUXUPI_nu(L,NW) + XX !same by spectral band. END DO ! print*,'-----------------------------------' ! print*,'nw=',nw ! print*,'ng=',ng ! print*,'FMDI=',FMDI(L_NLEVRAD-1) ! print*,'FMUPI=',FMUPI(L_NLEVRAD-1) ! print*,'-----------------------------------' 501 CONTINUE !End Spectral Interval LOOP C *** END OF MAJOR SPECTRAL INTERVAL LOOP IN THE INFRARED**** !print*,'-----------------------------------' !print*,'gweight=',gweight !print*,'FLUXDNI=',FLUXDNI(L_NLEVRAD-1) !print*,'FLUXUPI=',FLUXUPI(L_NLEVRAD-1) !print*,'-----------------------------------' !do nw=1,L_NSPECTI ! print*,'fup_tmp(',nw,')=',fup_tmp(nw) ! print*,'fdn_tmp(',nw,')=',fdn_tmp(nw) !enddo END SUBROUTINE SFLUXI_PLUTO end module sfluxi_pluto_mod