MODULE phyparam_mod
  USE callkeys
  USE comgeomfi
  IMPLICIT NONE
  PRIVATE
  SAVE

  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
    
  REAL, ALLOCATABLE :: tsurf(:),tsoil(:,:),rnatur(:), &
       capcal(:),fluxgrd(:), &
       dtrad(:,:),fluxrad(:), &
       q2(:,:),q2l(:,:), &
       albedo(:),emissiv(:),z0(:),inertie(:)
  !$OMP THREADPRIVATE( tsurf,tsoil,rnatur)
  !$OMP THREADPRIVATE( capcal,fluxgrd,dtrad,fluxrad)
  !$OMP THREADPRIVATE( q2,q2l)
  !$OMP THREADPRIVATE( 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)
    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:
    !   raditive transfer (long and shortwave) for CO2 and dust.
    !   vertical turbulent mixing
    !   convective adjsutment
    !
    !   author: Frederic Hourdin 15 / 10 /93
    !=======================================================================

    INTEGER, INTENT(IN) ::      &
         ngrid,                 & ! Size of the horizontal grid.
         nlayer                   ! Number of vertical layers.
    LOGICAL, INTENT(IN) ::      &
         firstcall,             & ! True at the first call 
         lastcall                 ! True at the last call
    REAL, INTENT(IN)    ::      &
         rjourvrai,             & ! Number of days counted from the North. Spring equinox
         gmtime,                & ! time of the day in seconds
         ptimestep,             & ! timestep (s)
         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,fract
    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)
    
    print*,'OK DANS PHYPARAM'
    print*,'latitude0',ngrid,lati(1:2),lati(ngrid-1:ngrid)
    print*,'nlayer',nlayer
    
    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_phyparam(ngrid, nlayer, igout)

       !     print*,'OK PHYPARAM 1 '
       IF(callsoil) THEN
          CALL soil(ngrid,nsoilmx,firstcall,inertie, &
               &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
          ! NB : this call to soil also performs some calculations, see surface.F90
       ELSE
          PRINT*,'WARNING!!! Thermal conduction in the soil turned off'
          DO ig=1,ngrid
             capcal(ig)  = capcal_nosoil
             fluxgrd(ig) = 0.
          ENDDO
       ENDIF
    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. Calcul of the radiative tendencies :
    !    ---------------------------------------
    
    IF(callrad) CALL radiative_tendencies(ngrid, igout, nlayer, &
         gmtime, ptimestep*float(iradia), zday, pplev, pplay, pt, &
         pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, 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,q2,q2l,   &
       &        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
          PRINT*,'Surface Heat capacity,conduction Flux, Ts, dTs, dt'
          PRINT*,capcal(igout),fluxgrd(igout),tsurf(igout), &
               &          zdtsrf(igout),ptimestep
       ENDIF
    ENDIF
    
    !-----------------------------------------------------------------------
    !   sorties:
    !   --------
    
    print*,'zday, zday_last ',zday,zday_last,icount
    if(abs(zday-zday_last-period_sort)<=ptimestep/unjours/10.) then
       print*,'zday, zday_last SORTIE ',zday,zday_last
       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('fract','fract',' ',fract)
       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, 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), zfluxlw(ngrid), &
         zdtsw(ngrid,nlayer), zfluxsw(ngrid), mu0(ngrid)
    REAL, DIMENSION(ngrid) :: fract
    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)
          print*,'ZENANG '
       ENDIF
       
       IF(lverbose) THEN
          PRINT*,'day, declin, sinlon,coslon,sinlat,coslat'
          PRINT*,zday, declin,                       &
               &            sinlon(igout),coslon(igout),  &
               &            sinlat(igout),coslat(igout)
       ENDIF
    ELSE
       print*,'declin,ngrid,planet_rad',declin,ngrid,planet_rad
       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
       PRINT*,'Diagnostics for radiation'
       PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck'
       PRINT*,albedo(igout),emissiv(igout),mu0(igout), &
            &           fract(igout),                       &
            &           fluxrad(igout),zplanck(igout)
       PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)'
       PRINT*,'unjours',unjours
       DO l=1,nlayer
          WRITE(*,'(3f15.5,2e15.2)') pt(igout,l),   &
               &         pplay(igout,l),pplev(igout,l),  &
               &         zdtsw(igout,l),zdtlw(igout,l)
       ENDDO
    ENDIF
    
  END SUBROUTINE radiative_tendencies

  SUBROUTINE alloc_phyparam(ngrid, nlayer, igout)
    USE surface
    USE astronomy, ONLY : iniorbit
    INTEGER, INTENT(IN) :: ngrid, nlayer, igout
    LOGICAL, PARAMETER :: firstcall=.TRUE.

    print*,'AKk',ngrid,nsoilmx
    allocate(tsurf(ngrid))
    print*,'AKa'
    allocate (tsoil(ngrid,nsoilmx))
    print*,'AKb'
    allocate (rnatur(ngrid))
    print*,'AK2'
    allocate(capcal(ngrid),fluxgrd(ngrid))
    print*,'AK3'
    allocate(dtrad(ngrid,nlayer),fluxrad(ngrid))
    print*,'AK4'
    allocate(q2(ngrid,nlayer+1),q2l(ngrid,nlayer+1))
    print*,'AK5'
    allocate(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid))
    print*,'AK6'
    
    
    PRINT*,'FIRSTCALL  '
    
    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
    q2=1.e-10
    q2l=1.e-10
    tsurf=300.
    tsoil=300.
    print*,tsoil(igout,nsoilmx/2+2)
    print*,'TS ',tsurf(igout),tsoil(igout,5)
    CALL iniorbit
    
    if (.not.callrad) fluxrad(:)=0.
    
    icount=0
    
  END SUBROUTINE alloc_phyparam

END MODULE phyparam_mod