source: dynamico_lmdz/simple_physics/phyparam/physics/phyparam_mod.F90 @ 4240

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