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,nq, &
       &            firstcall,lastcall, &
       &            rjourvrai,gmtime,ptimestep, &
       &            pplev,pplay,pphi,pphis,presnivs, &
       &            pu,pv,pt,pq, &
       &            pw, &
       &            pdu,pdv,pdt,pdq,pdpsrf)
    USE phys_const, ONLY : g, rcp, r, unjours
    USE surface,    ONLY : soil
    USE turbulence, ONLY : vdif
    USE convection, ONLY : convadj

    !=======================================================================
    !   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.
         nq                       ! Number of advected fields (tracers)
    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)
         pphis(ngrid),          & ! surface geopotential (unused)
         presnivs(nlayer),      & 
         pu(ngrid,nlayer),      & ! u component of the wind (ms-1)
         pv(ngrid,nlayer),      & ! v component of the wind (ms-1)
         pw(ngrid,nlayer),      & ! vertical velocity (unused)
         pt(ngrid,nlayer),      & ! Temperature (K)
         pq(ngrid,nlayer,nq)      ! Advected fields (unused)
    REAL, INTENT(OUT)   ::      & ! output : physical tendencies
         pdu(ngrid,nlayer),     &
         pdv(ngrid,nlayer),     &
         pdt(ngrid,nlayer),     &
         pdq(ngrid,nlayer,nq),  &
         pdpsrf(ngrid)
    
    !    Local variables :
    REAL, DIMENSION(ngrid) :: mu0,fract
    INTEGER :: j,l,ig,ierr,aslun,nlevel,igout,it1,it2,isoil,iq
    INTEGER*4 day_ini
    !
    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)
    REAL factq(nq),tauq(nq)
    character*3 nomq
    
    !   Local saved variables:
    !   ----------------------
    
    print*,'OK DANS PHYPARAM'
    print*,'nq ',nq
    print*,'latitude0',ngrid,lati(1:2),lati(ngrid-1:ngrid)
    print*,'nlayer',nlayer
    print*,'size pdq ',ngrid*nlayer*4,ngrid*nlayer*nq, &
         &      size(pdq),size(lati),size(pq),size(factq)
    
    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
       
       PRINT*,'FIRSTCALL B '
       print*,'INIIO avant iophys_ini '
       call iophys_ini('phys.nc    ',presnivs)       
    ENDIF

    !-----------------------------------------------------------------------
    !    1. Initialisations :
    !    --------------------
    
    icount=icount+1
    
    pdq(:,:,:) = 0. ! we do not use tracers in this physics package    
    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 iophys_ecrit('u',nlayer,'Vent zonal moy','m/s',pu)
       call iophys_ecrit('v',nlayer,'Vent meridien moy','m/s',pv)
       call iophys_ecrit('temp',nlayer,'Temperature','K',pt)
       call iophys_ecrit('geop',nlayer,'Geopotential','m2/s2',pphi)
       call iophys_ecrit('plev',nlayer,'plev','Pa',pplev(:,1:nlayer))
       
       call iophys_ecrit('du',nlayer,'du',' ',pdu)
       call iophys_ecrit('dv',nlayer,'du',' ',pdv)
       call iophys_ecrit('dt',nlayer,'du',' ',pdt)
       call iophys_ecrit('dtsw',nlayer,'dtsw',' ',zdtsw)
       call iophys_ecrit('dtlw',nlayer,'dtlw',' ',zdtlw)
       
       do iq=1,nq
          nomq="tr."
          write(nomq(2:3),'(i1.1)') iq
          call iophys_ecrit(nomq,nlayer,nomq,' ',pq(:,:,iq))
       enddo
       
       call iophys_ecrit('ts',1,'Surface temper','K',tsurf)
       call iophys_ecrit('coslon',1,'coslon',' ',coslon)
       call iophys_ecrit('sinlon',1,'sinlon',' ',sinlon)
       call iophys_ecrit('coslat',1,'coslat',' ',coslat)
       call iophys_ecrit('sinlat',1,'sinlat',' ',sinlat)
       call iophys_ecrit('mu0',1,'mu0',' ',mu0)
       call iophys_ecrit('alb',1,'alb',' ',albedo)
       call iophys_ecrit('fract',1,'fract',' ',fract)
       call iophys_ecrit('ps',1,'Surface pressure','Pa',pplev)
       call iophys_ecrit('slp',1,'Sea level pressure','Pa',zpmer)
       call iophys_ecrit('swsurf',1,'SW surf','Pa',zfluxsw)
       call iophys_ecrit('lwsurf',1,'LW surf','Pa',zfluxlw)
       
    endif
    
    !-----------------------------------------------------------------------
    IF(lastcall) THEN
       call iotd_fin
       PRINT*,'Ecriture du fichier de reinitialiastion de la physique'
       write(75)  tsurf,tsoil
    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