MODULE phyparam_mod #include "use_logging.h" USE callkeys USE comgeomfi IMPLICIT NONE PRIVATE SAVE REAL :: nan ! initialized to NaN with adequate compiler options REAL, PARAMETER :: pi=2*ASIN(1.), solarcst=1370., stephan=5.67e-08, & ps_rad=1.e5, height_scale=10000., ref_temp=285., & capcal_nosoil=1e5, tsoil_init=300. ! internal state, written to / read from disk at checkpoint / restart REAL, ALLOCATABLE :: tsurf(:), tsoil(:,:), capcal(:), fluxgrd(:) !$OMP THREADPRIVATE(tsurf, tsoil, capcal, fluxgrd) ! precomputed arrays REAL, ALLOCATABLE :: rnatur(:), albedo(:),emissiv(:), z0(:), inertie(:) !$OMP THREADPRIVATE( rnatur, albedo, emissiv, z0, inertie) INTEGER :: icount REAL :: zday_last !$OMP THREADPRIVATE( icount,zday_last) PUBLIC :: phyparam CONTAINS SUBROUTINE phyparam(ngrid,nlayer, & & firstcall,lastcall, & & rjourvrai,gmtime,ptimestep, & & pplev,pplay,pphi, & & pu,pv,pt, & & pdu,pdv,pdt,pdpsrf) BIND(C, name='phyparam_phyparam') USE phys_const, ONLY : g, rcp, r, unjours USE surface, ONLY : soil USE turbulence, ONLY : vdif USE convection, ONLY : convadj USE writefield_mod, ONLY : writefield !======================================================================= ! Top routine of the physical parametrisations of the LMD ! 20 parameters GCM for planetary atmospheres. ! It includes: ! radiative transfer (long and shortwave) for CO2 and dust. ! vertical turbulent mixing ! convective adjsutment ! heat diffusion in the soil ! ! author: Frederic Hourdin 15 / 10 /93 !======================================================================= INTEGER, INTENT(IN), VALUE :: & ngrid, & ! Size of the horizontal grid. nlayer ! Number of vertical layers. LOGICAL, INTENT(IN), VALUE :: & firstcall, & ! True at the first call lastcall ! True at the last call REAL, INTENT(IN), VALUE :: & rjourvrai, & ! Number of days counted from the North. Spring equinox gmtime, & ! time of the day in seconds ptimestep ! timestep (s) REAL, INTENT(IN) :: & pplev(ngrid,nlayer+1), & ! Pressure at interfaces between layers (pa) pplay(ngrid,nlayer), & ! Pressure at the middle of the layers (Pa) pphi(ngrid,nlayer), & ! Geopotential at the middle of the layers (m2s-2) pu(ngrid,nlayer), & ! u component of the wind (ms-1) pv(ngrid,nlayer), & ! v component of the wind (ms-1) pt(ngrid,nlayer) ! Temperature (K) REAL, INTENT(OUT) :: & ! output : physical tendencies pdu(ngrid,nlayer), & pdv(ngrid,nlayer), & pdt(ngrid,nlayer), & pdpsrf(ngrid) ! Local variables : REAL, DIMENSION(ngrid) :: mu0 INTEGER :: j,l,ig,nlevel,igout ! REAL :: zday, zdtime REAL zh(ngrid,nlayer),z1,z2 REAL zzlev(ngrid,nlayer+1),zzlay(ngrid,nlayer) REAL zdvfr(ngrid,nlayer),zdufr(ngrid,nlayer) REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid) REAL zflubid(ngrid),zpmer(ngrid) REAL zpopsk(ngrid,nlayer) REAL zdum1(ngrid,nlayer) REAL zdum2(ngrid,nlayer) REAL zdum3(ngrid,nlayer) REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer) REAL zfluxsw(ngrid),zfluxlw(ngrid), fluxrad(ngrid) WRITELOG(*,*) 'latitude0', ngrid, lati(1:2), lati(ngrid-1:ngrid) WRITELOG(*,*) 'nlayer',nlayer LOG_DBG('phyparam') IF (ngrid.NE.ngridmax) THEN PRINT*,'STOP in inifis' PRINT*,'Probleme de dimenesions :' PRINT*,'ngrid = ',ngrid PRINT*,'ngridmax = ',ngridmax STOP ENDIF nlevel=nlayer+1 igout=ngrid/2+1 zday=rjourvrai+gmtime !----------------------------------------------------------------------- ! 0. Allocate and initialize at first call ! -------------------- IF(firstcall) THEN ! zday_last=rjourvrai zday_last=zday-ptimestep/unjours CALL alloc(ngrid, nlayer) CALL precompute CALL coldstart(ngrid, ptimestep) ENDIF !----------------------------------------------------------------------- ! 1. Initialisations : ! -------------------- icount=icount+1 pdv(:,:) = 0. pdu(:,:) = 0. pdt(:,:) = 0. pdpsrf(:) = 0. zflubid(:)= 0. zdtsrf(:) = 0. !----------------------------------------------------------------------- ! calcul du geopotentiel aux niveaux intercouches ! ponderation des altitudes au niveau des couches en dp/p DO l=1,nlayer DO ig=1,ngrid zzlay(ig,l)=pphi(ig,l)/g ENDDO ENDDO DO ig=1,ngrid zzlev(ig,1)=0. ENDDO DO l=2,nlayer DO ig=1,ngrid z1=(pplay(ig,l-1)+pplev(ig,l))/(pplay(ig,l-1)-pplev(ig,l)) z2=(pplev(ig,l)+pplay(ig,l))/(pplev(ig,l)-pplay(ig,l)) zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2) ENDDO ENDDO !----------------------------------------------------------------------- ! Transformation de la temperature en temperature potentielle DO l=1,nlayer DO ig=1,ngrid zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp ! surface pressure is used as reference pressure zh(ig,l)=pt(ig,l)/zpopsk(ig,l) ENDDO ENDDO !----------------------------------------------------------------------- ! 2. Compute radiative tendencies : ! --------------------------------------- IF(callrad) CALL radiative_tendencies(ngrid, igout, nlayer, & gmtime, ptimestep*float(iradia), zday, pplev, pplay, pt, & pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, fluxrad, mu0) !----------------------------------------------------------------------- ! 3. Vertical diffusion (turbulent mixing): ! ----------------------------------------- ! IF(calldifv) THEN DO ig=1,ngrid zflubid(ig)=fluxrad(ig)+fluxgrd(ig) ENDDO zdum1(:,:)=0. zdum2(:,:)=0. do l=1,nlayer do ig=1,ngrid zdum3(ig,l)=pdt(ig,l)/zpopsk(ig,l) enddo enddo CALL vdif(ngrid,nlayer,zday, & & ptimestep,capcal,z0, & & pplay,pplev,zzlay,zzlev, & & pu,pv,zh,tsurf,emissiv, & & zdum1,zdum2,zdum3,zflubid, & & zdufr,zdvfr,zdhfr,zdtsrfr, & & lverbose) DO l=1,nlayer DO ig=1,ngrid pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l) pdu(ig,l)=pdu(ig,l)+zdufr(ig,l) pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l) ENDDO ENDDO DO ig=1,ngrid zdtsrf(ig)=zdtsrf(ig)+zdtsrfr(ig) ENDDO ELSE DO ig=1,ngrid zdtsrf(ig)=zdtsrf(ig)+ & & (fluxrad(ig)+fluxgrd(ig))/capcal(ig) ENDDO ENDIF ! !----------------------------------------------------------------------- ! 4. Dry convective adjustment: ! ----------------------------- IF(calladj) THEN DO l=1,nlayer DO ig=1,ngrid zdum1(ig,l)=pdt(ig,l)/zpopsk(ig,l) ENDDO ENDDO zdufr(:,:)=0. zdvfr(:,:)=0. zdhfr(:,:)=0. CALL convadj(ngrid,nlayer,ptimestep, & & pplay,pplev,zpopsk, & & pu,pv,zh, & & pdu,pdv,zdum1, & & zdufr,zdvfr,zdhfr) DO l=1,nlayer DO ig=1,ngrid pdu(ig,l)=pdu(ig,l)+zdufr(ig,l) pdv(ig,l)=pdv(ig,l)+zdvfr(ig,l) pdt(ig,l)=pdt(ig,l)+zdhfr(ig,l)*zpopsk(ig,l) ENDDO ENDDO ENDIF !----------------------------------------------------------------------- ! On ajoute les tendances physiques a la temperature du sol: ! --------------------------------------------------------------- DO ig=1,ngrid tsurf(ig)=tsurf(ig)+ptimestep*zdtsrf(ig) ENDDO WRITE(55,'(2e15.5)') zday,tsurf(ngrid/2+1) !----------------------------------------------------------------------- ! soil temperatures: ! -------------------- IF (callsoil) THEN CALL soil(ngrid,nsoilmx,.false.,inertie, & & ptimestep,tsurf,tsoil,capcal,fluxgrd) IF(lverbose) THEN WRITELOG(*,*) 'Surface Heat capacity,conduction Flux, Ts, dTs, dt' WRITELOG(*,*) capcal(igout), fluxgrd(igout), tsurf(igout), & & zdtsrf(igout), ptimestep LOG_DBG('phyparam') ENDIF ENDIF !----------------------------------------------------------------------- ! sorties: ! -------- WRITELOG(*,*) 'zday, zday_last ',zday,zday_last,icount LOG_DBG('phyparam') if(abs(zday-zday_last-period_sort)<=ptimestep/unjours/10.) then WRITELOG(*,*) 'zday, zday_last SORTIE ', zday, zday_last LOG_INFO('phyparam') zday_last=zday ! Ecriture/extension de la coordonnee temps do ig=1,ngridmax zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(r*ref_temp)) enddo call writefield('u','Vent zonal moy','m/s',pu) call writefield('v','Vent meridien moy','m/s',pv) call writefield('temp','Temperature','K',pt) call writefield('geop','Geopotential','m2/s2',pphi) call writefield('plev','plev','Pa',pplev(:,1:nlayer)) call writefield('du','du',' ',pdu) call writefield('dv','du',' ',pdv) call writefield('dt','du',' ',pdt) call writefield('dtsw','dtsw',' ',zdtsw) call writefield('dtlw','dtlw',' ',zdtlw) call writefield('ts','Surface temper','K',tsurf) call writefield('coslon','coslon',' ',coslon) call writefield('sinlon','sinlon',' ',sinlon) call writefield('coslat','coslat',' ',coslat) call writefield('sinlat','sinlat',' ',sinlat) call writefield('mu0','mu0',' ',mu0) call writefield('alb','alb',' ',albedo) call writefield('ps','Surface pressure','Pa',pplev(:,1)) call writefield('slp','Sea level pressure','Pa',zpmer) call writefield('swsurf','SW surf','Pa',zfluxsw) call writefield('lwsurf','LW surf','Pa',zfluxlw) endif END SUBROUTINE phyparam SUBROUTINE radiative_tendencies(ngrid, igout, nlayer, & gmtime, zdtime, zday, pplev, pplay, pt, & pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, fluxrad, mu0) USE planet USE phys_const, ONLY : planet_rad, unjours USE astronomy, ONLY : orbite, solarlong USE solar, ONLY : solang, zenang, mucorr USE radiative_sw, ONLY : sw USE radiative_lw, ONLY : lw INTEGER, INTENT(IN) :: ngrid, igout, nlayer REAL, INTENT(IN) :: gmtime, zdtime, zday, & & pplev(ngrid,nlayer+1), pplay(ngrid, nlayer), & & pt(ngrid, nlayer+1) REAL, INTENT(INOUT) :: pdt(ngrid,nlayer) REAL, INTENT(OUT) :: zdtlw(ngrid,nlayer), & ! long-wave temperature tendency & zfluxlw(ngrid), & ! long-wave flux at surface & zdtsw(ngrid,nlayer), & ! short-wave temperature tendency & zfluxsw(ngrid), & ! short-wave flux at surface & fluxrad(ngrid), & ! surface flux mu0(ngrid) ! ?? ! local variables REAL :: fract(ngrid), dtrad(ngrid, nlayer) REAL :: zls, zinsol, tsurf2, ztim1,ztim2,ztim3, dist_sol, declin REAL :: zplanck(ngrid) INTEGER :: ig, l ! 2.1 Insolation ! -------------------------------------------------- CALL solarlong(zday,zls) CALL orbite(zls,dist_sol,declin) IF(diurnal) THEN IF ( .TRUE. ) then ztim1=SIN(declin) ztim2=COS(declin)*COS(2.*pi*(zday-.5)) ztim3=-COS(declin)*SIN(2.*pi*(zday-.5)) CALL solang(ngrid,sinlon,coslon,sinlat,coslat, & & ztim1,ztim2,ztim3, & & mu0,fract) ELSE CALL zenang(ngrid,zls,gmtime,zdtime,lati,long,mu0,fract) ENDIF IF(lverbose) THEN WRITELOG(*,*) 'day, declin, sinlon,coslon,sinlat,coslat' WRITELOG(*,*) zday, declin, & & sinlon(igout),coslon(igout), & & sinlat(igout),coslat(igout) LOG_DBG('radiative_tendencies') ENDIF ELSE WRITELOG(*,*) 'declin,ngrid,planet_rad', declin, ngrid, planet_rad LOG_DBG('radiative_tendencies') CALL mucorr(ngrid,declin,lati,mu0,fract,height_scale,planet_rad) ENDIF zinsol=solarcst/(dist_sol*dist_sol) ! 2.2 Radiative tendencies and fluxes: ! -------------------------------------------------- CALL sw(ngrid,nlayer,diurnal,coefvis,albedo, & & pplev,ps_rad, & & mu0,fract,zinsol, & & zfluxsw,zdtsw, & & lverbose) CALL lw(ngrid,nlayer,coefir,emissiv, & & pplev,ps_rad,tsurf,pt, & & zfluxlw,zdtlw, & & lverbose) ! 2.4 surface fluxes ! ------------------------------ DO ig=1,ngrid fluxrad(ig)=emissiv(ig)*zfluxlw(ig) & & +zfluxsw(ig)*(1.-albedo(ig)) tsurf2 = tsurf(ig)*tsurf(ig) zplanck(ig)=emissiv(ig)*stephan*tsurf2*tsurf2 fluxrad(ig)=fluxrad(ig)-zplanck(ig) ENDDO ! 2.5 Temperature tendencies ! -------------------------- DO l=1,nlayer DO ig=1,ngrid dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l) pdt(ig,l)=pdt(ig,l)+dtrad(ig,l) ENDDO ENDDO IF(lverbose) THEN WRITELOG(*,*) 'Diagnostics for radiation' WRITELOG(*,*) 'albedo, emissiv, mu0,fract,Frad,Planck' WRITELOG(*,*) albedo(igout),emissiv(igout),mu0(igout), & & fract(igout), & & fluxrad(igout),zplanck(igout) WRITELOG(*,*) 'Tlay Play Plev dT/dt SW dT/dt LW (K/day)' WRITELOG(*,*) 'unjours',unjours DO l=1,nlayer WRITELOG(*,'(3f15.5,2e15.2)') pt(igout,l), & & pplay(igout,l),pplev(igout,l), & & zdtsw(igout,l),zdtlw(igout,l) ENDDO LOG_DBG('radiative_tendencies') ENDIF END SUBROUTINE radiative_tendencies SUBROUTINE alloc(ngrid, nlayer) BIND(C, name='phyparam_alloc') !$cython header void phyparam_alloc(int, int); !$cython wrapper def alloc(ngrid, nlayer) : phy.phyparam_alloc(ngrid, nlayer) USE astronomy, ONLY : iniorbit USE surface, ONLY : zc,zd INTEGER, INTENT(IN), VALUE :: ngrid, nlayer LOGICAL, PARAMETER :: firstcall=.TRUE. ! allocate arrays for internal state ALLOCATE(tsurf(ngrid)) ALLOCATE(tsoil(ngrid,nsoilmx)) ! we could avoid the arrays below with a different implementation of surface / radiation / turbulence coupling ALLOCATE(capcal(ngrid),fluxgrd(ngrid)) ALLOCATE(zc(ngrid,nsoilmx),zd(ngrid,nsoilmx)) ! allocate precomputed arrays ALLOCATE(rnatur(ngrid)) ALLOCATE(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid)) CALL iniorbit END SUBROUTINE alloc SUBROUTINE precompute() BIND(C, name='phyparam_precompute') !$cython header void phyparam_precompute(); !$cython wrapper def precompute() : phy.phyparam_precompute() USE surface ! precompute time-independent arrays rnatur(:) = 1. emissiv(:) = (1.-rnatur(:))*emi_mer+rnatur(:)*emi_ter inertie(:) = (1.-rnatur(:))*I_mer+rnatur(:)*I_ter albedo(:) = (1.-rnatur(:))*alb_mer+rnatur(:)*alb_ter z0(:) = (1.-rnatur(:))*Cd_mer+rnatur(:)*Cd_ter END SUBROUTINE precompute SUBROUTINE coldstart(ngrid, ptimestep) BIND(C, name='phyparam_coldstart') !$cython header void phyparam_coldstart(int, double); !$cython wrapper def coldstart (ngrid, timestep): phy.phyparam_coldstart(ngrid, timestep) ! create internal state to start a run without a restart file USE surface, ONLY : soil INTEGER, INTENT(IN), VALUE :: ngrid REAL, INTENT(IN), VALUE :: ptimestep tsurf(:) = tsoil_init tsoil(:,:) = tsoil_init icount=0 IF(callsoil) THEN ! initializes zc, zd, capcal, fluxgrd CALL soil(ngrid,nsoilmx,.TRUE.,inertie, & & ptimestep,tsurf,tsoil,capcal,fluxgrd) ELSE WRITELOG(*,*) 'WARNING!!! Thermal conduction in the soil turned off' LOG_WARN('coldstart') capcal(:) = capcal_nosoil fluxgrd(:) = 0. ENDIF END SUBROUTINE coldstart END MODULE phyparam_mod