! ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/radlwsw.F,v 1.2 2004/10/27 10:14:46 lmdzadmin Exp $ ! SUBROUTINE radlwsw(dist, rmu0, fract, zzlev, . paprs, pplay,tsol, t) c c====================================================================== c Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719 c Objet: interface entre le modele et les rayonnements c Arguments: c dist-----input-R- distance astronomique terre-soleil c rmu0-----input-R- cosinus de l'angle zenithal c fract----input-R- duree d'ensoleillement normalisee c zzlev----input-R- altitude a inter-couche (m) c paprs----input-R- pression a inter-couche (Pa) c pplay----input-R- pression au milieu de couche (Pa) c tsol-----input-R- temperature du sol (en K) c t--------input-R- temperature (K) c MODIFS pour multimatrices ksi SPECIFIQUE VENUS c S. Lebonnois 20/12/2006 c corrections 13/07/2007 c New ksi matrix: possibility of different cloud model fct of lat 05/2014 c With extension NLTE (G. Gilli, 2014) c====================================================================== use dimphy USE geometry_mod, ONLY: latitude USE phys_state_var_mod, only: heat,cool,radsol, . topsw,toplw,solsw,sollw,sollwdown,lwnet,swnet use write_field_phy IMPLICIT none #include "YOMCST.h" #include "clesphys.h" #include "comcstVE.h" #include "nlteparams.h" !=========== ! Arguments !=========== real rmu0(klon), fract(klon), dist REAL zzlev(klon,klev+1) real paprs(klon,klev+1), pplay(klon,klev) real tsol(klon) real t(klon,klev) !=========== ! Local !=========== INTEGER k, kk, i, j, band REAL PPB(klev+1) REAL zfract, zrmu0 REAL zheat(klev), zcool(klev) real temp(klev),znivs(klev+1) REAL ZFSNET(klev+1),ZFLNET(klev+1) REAL ztopsw, ztoplw REAL zsolsw, zsollw cIM BEG REAL zsollwdown cIM END real,save,allocatable :: ksive(:,:,:,:) ! ksi matrixes in Vincent's file real psi(0:klev+1,0:klev+1) real deltapsi(0:klev+1,0:klev+1) real pt0(0:klev+1) real bplck(0:klev+1,nnuve) ! Planck luminances in table layers real y(0:klev,nnuve) ! temporary variable for Planck real zdblay(0:klev+1,nnuve) ! temperature gradient of planck function integer mat0,lat,ips,isza,ips0,isza0 real factp,factz,ksi logical firstcall data firstcall/.true./ save firstcall c------------------------------------------- c Initialisations c----------------- if (firstcall) then c ---------- ksive -------------- allocate(ksive(0:klev+1,0:klev+1,nnuve,nbmat)) call load_ksi(ksive) endif ! firstcall c------------------------------------------- DO k = 1, klev DO i = 1, klon heat(i,k)=0. cool(i,k)=0. ENDDO ENDDO c+++++++ BOUCLE SUR LA GRILLE +++++++++++++++++++++++++ DO j = 1, klon c====================================================================== c Initialisations c --------------- DO k = 1, klev zheat(k) = 0.0 zcool(k) = 0.0 ENDDO DO k = 1, klev+1 ZFLNET(k) = 0.0 ZFSNET(k) = 0.0 ENDDO ztopsw = 0.0 ztoplw = 0.0 zsolsw = 0.0 zsollw = 0.0 zsollwdown = 0.0 zfract = fract(j) zrmu0 = rmu0(j) DO k = 1, klev+1 PPB(k) = paprs(j,k)/1.e5 ENDDO pt0(0) = tsol(j) DO k = 1, klev pt0(k) = t(j,k) ENDDO pt0(klev+1) = 0. DO k = 0,klev+1 DO i = 0,klev+1 psi(i,k) = 0. ! positif quand nrj de i->k deltapsi(i,k) = 0. ENDDO ENDDO c====================================================================== c Getting psi and deltapsi c ------------------------ c Planck function c --------------- do band=1,nnuve do k=0,klev c B(T,l) = al/(exp(bl/T)-1) y(k,band) = exp(bl(band)/pt0(k))-1. bplck(k,band) = al(band)/(y(k,band)) zdblay(k,band)= al(band)*bl(band)*exp(bl(band)/pt0(k))/ . ((pt0(k)*pt0(k))*(y(k,band)*y(k,band))) enddo bplck(klev+1,band) = 0.0 zdblay(klev+1,band)= 0.0 enddo c finding the right matrixes c -------------------------- mat0 = 0 if (nlatve.eq.1) then ! clouds are taken as uniform lat=1 else if (abs(latitude(j)).le.50.) then lat=1 elseif (abs(latitude(j)).le.60.) then lat=2 elseif (abs(latitude(j)).le.70.) then lat=3 elseif (abs(latitude(j)).le.80.) then lat=4 else lat=5 endif endif ips0=0 if (nbpsve(lat).gt.1) then do ips=1,nbpsve(lat)-1 if ( (psurfve(ips,lat).ge.paprs(j,1)) . .and.(psurfve(ips+1,lat).lt.paprs(j,1)) ) then ips0 = ips c print*,'ig=',j,' ips0=',ips factp = (paprs(j,1) -psurfve(ips0,lat)) . /(psurfve(ips0+1,lat)-psurfve(ips0,lat)) exit endif enddo else ! Only one ps, no interpolation ips0=1 endif isza0=0 if (nbszave(lat).gt.1) then do isza=1,nbszave(lat)-1 if ( (szave(isza,lat).ge.zrmu0) . .and.(szave(isza+1,lat).lt.zrmu0) ) then isza0 = isza c print*,'ig=',j,' isza0=',isza factz = (zrmu0 -szave(isza0,lat)) . /(szave(isza0+1,lat)-szave(isza0,lat)) exit endif enddo else ! Only one sza, no interpolation isza0=-99 endif c -------- Probleme aux bords if ((ips0.eq.0).and.(psurfve(1,lat).gt.paprs(j,1))) then ips0 = 1 print*,'Extrapolation! ig=',j,' ips0=',ips0 factp = (paprs(j,1) -psurfve(ips0,lat)) . /(psurfve(ips0+1,lat)-psurfve(ips0,lat)) endif if ((ips0.eq.0).and.(psurfve(nbpsve(lat),lat).le.paprs(j,1))) . then ips0 = nbpsve(lat)-1 print*,'Extrapolation! ig=',j,' ips0=',ips0 factp = (paprs(j,1) -psurfve(ips0,lat)) . /(psurfve(ips0+1,lat)-psurfve(ips0,lat)) endif c --------- if ((ips0.eq.0).or.(isza0.eq.0)) then write(*,*) 'Finding the right matrix in radlwsw' print*,'Interpolation problem, grid point ig=',j print*,'psurf = ',paprs(j,1),' mu0 = ',zrmu0 stop endif if (isza0.eq.-99) then mat0 = indexve(lat)+ips0 else mat0 = indexve(lat)+(isza0-1)*nbpsve(lat)+ips0 endif c interpolation of ksi and computation of psi,deltapsi c ---------------------------------------------------- do band=1,nnuve do k=0,klev+1 do i=0,klev+1 if (isza0.eq.-99) then ksi = ksive(i,k,band,mat0)*(1-factp) . +ksive(i,k,band,mat0+1)*factp else ksi = ksive(i,k,band,mat0)*(1-factp)*(1-factz) . +ksive(i,k,band,mat0+1)*factp *(1-factz) . +ksive(i,k,band,mat0+nbpsve(lat))*(1-factp)*factz . +ksive(i,k,band,mat0+nbpsve(lat)+1)*factp *factz endif psi(i,k) = psi(i,k) + . RPI*ksi*(bplck(i,band)-bplck(k,band)) deltapsi(i,k) = deltapsi(i,k) + RPI*ksi*zdblay(i,band) enddo enddo enddo c====================================================================== c LW call c--------- temp(1:klev)=t(j,1:klev) CALL LW_venus_ve( . PPB,temp,psi,deltapsi, . zcool, . ztoplw,zsollw, . zsollwdown,ZFLNET) c--------- c SW call c--------- znivs=zzlev(j,:) c CALL SW_venus_ve_1Dglobave(zrmu0, zfract, ! pour moy globale c CALL SW_venus_ve(zrmu0, zfract, c S PPB,temp,znivs, c S zheat, c S ztopsw,zsolsw,ZFSNET) c CALL SW_venus_cl_1Dglobave(zrmu0, zfract, ! pour moy globale c CALL SW_venus_cl(zrmu0, zfract, c CALL SW_venus_dc_1Dglobave(zrmu0, zfract, ! pour moy globale CALL SW_venus_dc(zrmu0, zfract, S PPB,temp, S zheat, S ztopsw,zsolsw,ZFSNET) c====================================================================== radsol(j) = zsolsw - zsollw ! + vers bas topsw(j) = ztopsw ! + vers bas toplw(j) = ztoplw ! + vers haut solsw(j) = zsolsw ! + vers bas sollw(j) = -zsollw ! + vers bas sollwdown(j) = zsollwdown ! + vers bas DO k = 1, klev+1 lwnet (j,k) = ZFLNET(k) swnet (j,k) = ZFSNET(k) ENDDO c C heat/cool with upper atmosphere C IF(callnlte) THEN DO k = 1,nlaylte heat(j,k) = zheat(k) cool(j,k) = zcool(k) ENDDO c Zero tendencies for any remaining layers between nlaylte and klev if (klev.gt.nlaylte) then do k = nlaylte+1, klev heat(j,k) = 0. cool(j,k) = 0. enddo endif ELSE DO k = 1, klev heat(j,k) = zheat(k) cool(j,k) = zcool(k) ENDDO ENDIF ! callnlte ENDDO ! of DO j = 1, klon c+++++++ FIN BOUCLE SUR LA GRILLE +++++++++++++++++++++++++ ! for tests: write output fields... ! call writefield_phy('radlwsw_heat',heat,klev) ! call writefield_phy('radlwsw_cool',cool,klev) ! call writefield_phy('radlwsw_radsol',radsol,1) ! call writefield_phy('radlwsw_topsw',topsw,1) ! call writefield_phy('radlwsw_toplw',toplw,1) ! call writefield_phy('radlwsw_solsw',solsw,1) ! call writefield_phy('radlwsw_sollw',sollw,1) ! call writefield_phy('radlwsw_sollwdown',sollwdown,1) ! call writefield_phy('radlwsw_swnet',swnet,klev+1) ! call writefield_phy('radlwsw_lwnet',lwnet,klev+1) c tests c j = klon/2 c j = 1 c print*,'mu0=',rmu0(j) c print*,' net flux vis HEAT(K/Eday)' c do k=1,klev c print*,k,ZFSNET(k),heat(j,k)*86400. c enddo c print*,' net flux IR COOL(K/Eday)' c do k=1,klev c print*,k,ZFLNET(k),cool(j,k)*86400. c enddo firstcall = .false. RETURN END