C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/verification/fizhi-gridalt-hs/code/fizhi_init_vars.F,v 1.3 2005/05/05 21:27:39 molod Exp $
C $Name:  $

#include "FIZHI_OPTIONS.h"
       subroutine fizhi_init_vars (myThid)
c-----------------------------------------------------------------------
c  Routine to initialise the fizhi state.
c  
c  Input: myThid       - Process number calling this routine
c
c  Notes: 
c   1) For a Cold Start - 
c      This routine takes the initial condition on the dynamics grid
c      and interpolates to the physics grid to initialize the state
c      variables that are on both grids. It initializes the variables
c      of the turbulence scheme to 0., and the land state from a model
c      climatology.
c   2) For a Restart, read the fizhi pickup file
c   3) The velocity component physics fields are on an A-Grid
c
c Calls: dyn2phys (x4)
c-----------------------------------------------------------------------
       implicit none
#include "SIZE.h"
#include "fizhi_SIZE.h"
#include "fizhi_land_SIZE.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "gridalt_mapping.h"
#include "fizhi_coms.h"
#include "fizhi_land_coms.h"
#include "fizhi_earth_coms.h"
#include "EEPARAMS.h"
#include "SURFACE.h"
#include "PARAMS.h"
#include "chronos.h"

       integer myThid

c pe on dynamics and physics grid refers to bottom edge
       _RL pephy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nrphys+1,nSx,nSy)
       _RL pedyn(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1,nSx,nSy)
       _RL windphy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nrphys,nSx,nSy)
       _RL udyntemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
       _RL vdyntemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
       _RL tempphy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nrphys,nSx,nSy)

       integer i, j, L, bi, bj, Lbotij
       integer im1, im2, jm1, jm2, idim1, idim2, jdim1, jdim2
       logical alarm
       external alarm

       im1 = 1-OLx
       im2 = sNx+OLx
       jm1 = 1-OLy
       jm2 = sNy+OLy
       idim1 = 1
       idim2 = sNx
       jdim1 = 1
       jdim2 = sNy

c   First Check to see if we can start a fizhi experiment at current time
c    All Fizhi alarms must be on for the first time step of a segment

      if( .not.alarm('moist') .or. .not.alarm('turb')   .or.
     .    .not.alarm('radsw') .or. .not.alarm('radlw') ) then
       print *,' Cant Start Fizhi experiment at ',nymd,' ',nhms
       stop
      endif

C Deal Here with Variables that are on a Fizhi Pickup or need Initialization

      IF ( startTime.EQ.baseTime .AND. nIter0.EQ.0 ) THEN
      print *,' In fizhi_init_vars: Beginning of New Experiment '

       do bj = myByLo(myThid), myByHi(myThid)
       do bi = myBxLo(myThid), myBxHi(myThid)

C Build pressures on dynamics grid
        do j = 1,sNy
        do i = 1,sNx
         do L = 1,Nr
          pedyn(i,j,L,bi,bj) = 0.
         enddo
        enddo
        enddo
        do j = 1,sNy
        do i = 1,sNx
         Lbotij = ksurfC(i,j,bi,bj)
         if(Lbotij.ne.0.) 
     .    pedyn(i,j,Lbotij,bi,bj) = Ro_surf(i,j,bi,bj) + etaH(i,j,bi,bj)
        enddo
        enddo
        do j = 1,sNy
        do i = 1,sNx
         Lbotij = ksurfC(i,j,bi,bj)
         do L = Lbotij+1,Nr+1
          pedyn(i,j,L,bi,bj) = pedyn(i,j,L-1,bi,bj) -
     .                        drF(L-1)*hfacC(i,j,L-1,bi,bj)
         enddo
c Do not use a zero field as the top edge pressure for interpolation
         if(pedyn(i,j,Nr+1,bi,bj).lt.1.e-5)
     .                               pedyn(i,j,Nr+1,bi,bj) = 1.e-5
        enddo
        enddo
C Build pressures on physics grid
        do j = 1,sNy
        do i = 1,sNx
         pephy(i,j,1,bi,bj)=Ro_surf(i,j,bi,bj) + etaH(i,j,bi,bj)
         do L = 2,Nrphys+1
          pephy(i,j,L,bi,bj)=pephy(i,j,L-1,bi,bj)-dpphys0(i,j,L-1,bi,bj)
         enddo
c Do not use a zero field as the top edge pressure for interpolation
         if(pephy(i,j,Nrphys+1,bi,bj).lt.1.e-5)
     .                               pephy(i,j,Nrphys+1,bi,bj) = 1.e-5
        enddo
        enddo
c
c Create an initial wind magnitude field on the physics grid -
c   Use a log wind law with z0=1cm, u*=1 cm/sec,
c   do units and get u = .025*ln(dP*10), with dP in pa.
        do L = 1,Nrphys
        do j = 1,sNy
        do i = 1,sNx
         windphy(i,j,L,bi,bj) = 0.025 *
     .             log((pephy(i,j,1,bi,bj)-pephy(i,j,L+1,bi,bj))*10.)
        enddo
        enddo
        enddo

       enddo
       enddo
                                                                                   
c Create initial fields on phys. grid - Move Dynamics u and v to A-Grid
       call CtoA(myThid,uvel,vvel,maskW,maskS,im1,im2,jm1,jm2,Nr,
     .                     Nsx,Nsy,1,sNx,1,sNy,udyntemp,vdyntemp)

       do bj = myByLo(myThid), myByHi(myThid)
       do bi = myBxLo(myThid), myBxHi(myThid)

c Create initial fields on phys. grid - interpolate from dyn. grid
        call dyn2phys(udyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,
     . 1,sNx,1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,1,tempphy)
c   Note: Interpolation gives bottom-up arrays (level 1 is bottom),
c         Physics works top-down. so -> need to flip arrays
        do L = 1,Nrphys
        do j = 1,sNy
        do i = 1,sNx
         uphy(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
        enddo
        enddo
        enddo
        call dyn2phys(vdyntemp,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,
     . 1,sNx,1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,1,tempphy)
        do L = 1,Nrphys
        do j = 1,sNy
        do i = 1,sNx
         vphy(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
        enddo
        enddo
        enddo
        call dyn2phys(theta,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,
     . 1,sNx,1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,tempphy)
        do L = 1,Nrphys
        do j = 1,sNy
        do i = 1,sNx
         thphy(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
        enddo
        enddo
        enddo

        call dyn2phys(salt,pedyn,im1,im2,jm1,jm2,Nr,Nsx,Nsy,
     . 1,sNx,1,sNy,bi,bj,windphy,pephy,ksurfC,Nrphys,nlperdyn,0,tempphy)
        do L = 1,Nrphys
        do j = 1,sNy
        do i = 1,sNx
         sphy(i,j,Nrphys+1-L,bi,bj) = tempphy(i,j,L,bi,bj)
        enddo
        enddo
        enddo

c Zero out fizhi tendency arrays on the fizhi grid
        do L = 1,Nrphys
        do j = 1,sNy
        do i = 1,sNx
         duphy(i,j,L,bi,bj) = 0.
         dvphy(i,j,L,bi,bj) = 0.
         dthphy(i,j,L,bi,bj) = 0.
         dsphy(i,j,L,bi,bj) = 0.
        enddo
        enddo
        enddo

c Zero out fizhi tendency arrays on the dynamics grid
        do L = 1,Nr
        do j = jm1,jm2
        do i = im1,im2
         guphy(i,j,L,bi,bj) = 0.
         gvphy(i,j,L,bi,bj) = 0.
         gthphy(i,j,L,bi,bj) = 0.
         gsphy(i,j,L,bi,bj) = 0.
        enddo
        enddo
        enddo

c Initialize vegetation tile tke, xlmt, khmt, xxmt, yymt, ctmt, zetamt,
        if( (nhms.eq.nhms0) .and. (nymd.eq.nymd0) ) then
         print *,' Cold Start: Zero out Turb second moments '
         do i = 1,nchp
          ctmt(i,bi,bj) = 0.
          xxmt(i,bi,bj) = 0.
          yymt(i,bi,bj) = 0.
          zetamt(i,bi,bj) = 0.
         enddo
         do L = 1,Nrphys
         do i = 1,nchp
          tke(i,L,bi,bj) = 0.
          xlmt(i,L,bi,bj) = 0.
          khmt(i,L,bi,bj) = 0.
         enddo
         enddo
        else
         print *,' Need initial Values for TKE - dont have them! '
         stop
        endif

c Now initialize vegetation tile land state too - tcanopy, etc... 
c       call fizhi_init_vegsurftiles( nymd,nhms, 'D', myThid )
c Now initialize land state too - tcanopy, etc... SET FOR NOW,
c                                              READ CLIM FOR REAL
        do i = 1,nchp
        tcanopy(i,bi,bj) = 283.
        tdeep(i,bi,bj) = 282.5
        ecanopy(i,bi,bj) = 2.e-2
        swetshal(i,bi,bj) = 0.6
        swetroot(i,bi,bj) = 0.5
        swetdeep(i,bi,bj) = 0.5
        capac(i,bi,bj) = 0.
        snodep(i,bi,bj) = 0.
        enddo

c Now initialize fizhi arrays that will be on a pickup
        print *,' Initialize fizhi arrays that will be on pickup '
        imstturblw(bi,bj) = 0
        imstturbsw(bi,bj) = 0
        iras(bi,bj) = 0
        nlwcld(bi,bj) = 0
        nlwlz(bi,bj) = 0
        nswcld(bi,bj) = 0
        nswlz(bi,bj) = 0
        do L = 1,Nrphys
        do j = 1,sNy
        do i = 1,sNx
         swlz(i,j,L,bi,bj) = 0.
         lwlz(i,j,L,bi,bj) = 0.
         qliqavesw(i,j,L,bi,bj) = 0.
         qliqavelw(i,j,L,bi,bj) = 0.
         fccavesw(i,j,L,bi,bj) = 0.
         fccavelw(i,j,L,bi,bj) = 0.
         cldtot_sw(i,j,L,bi,bj) = 0.
         cldras_sw(i,j,L,bi,bj) = 0.
         cldlsp_sw(i,j,L,bi,bj) = 0.
         cldtot_lw(i,j,L,bi,bj) = 0.
         cldras_lw(i,j,L,bi,bj) = 0.
         cldlsp_lw(i,j,L,bi,bj) = 0.
        enddo
        enddo
        enddo
        do j = 1,sNy
        do i = 1,sNx
         rainlsp(i,j,bi,bj) = 0.
         raincon(i,j,bi,bj) = 0.
         snowfall(i,j,bi,bj) = 0.
        enddo
        enddo

       enddo
       enddo

      ELSE
      print *,' In fizhi_init_vars: Read from restart '
                                                                                   
C--   Read fizhi package state variables from pickup file

       call fizhi_read_pickup( nIter0, myThid )
       CALL FIZHI_READ_VEGTILES( nIter0, 'D', myThid )

      ENDIF

       return
       end