program layer_field_10ann_climatology
c@***h* OCEAN_HEAT_CONTENT/layer_field_10ann_climatology.f
c
c NAME
c layer_field_10ann_climatology.f
c
c PURPOSE
c Specifically for OCEAN_HEAT_CONTENT (OHC), OHC in a given
c layer of the ocean is calculated using temperature anomaly
c (TA) objectively analyzed fields (calculated in analysis.f)
c the heat capacity of sea water and the density of sea water
c over the area of each latitude/longitude grid box over the
c given depth range. 
c
c DESCRIPTION
c layer_field_10ann_climatology.f is a multipurpose program which
c calculates an ocean variable over a volume of water.  A gridded 
c field is read in for each standard depth level within the
c designated depth range, the ocean variable is calculated
c over the area of each grid box and a depth from halfway
c to the preceding depth level to halfway to the next
c depth level (in the case of the first depth level only
c the second half depth, in the case of the last depth only
c the first half depth) using the input value at each 
c given depth level.  Resultant values are summed (or
c in some cases averaged) over each grid box.
c
c INPUT
c User input for layer_field_10ann_climatology.f comes from
c 'layer_field_10ann_climatology.inf', including designated standard
c depth layers, time period, temporal span (monthly quarterly, yearly,
c pentadal), input and output file forms, latitude/longitude span. 
c
c OUTPUT
c Files of OHC for user designated layer in units of 10^18 joules.
c Files are in the World Ocean Atlas (WOA) binary form.
c
cAUTHOR
c John Antonov for National Centers for Environmental Information (NCEI)
c
cCOPYRIGHT
cTHIS SOFTWARE AND ITS DOCUMENTATION ARE CONSIDERED TO BE IN THE
cPUBLIC DOMAIN AND THUS ARE AVAILABLE FOR UNRESTRICTED PUBLIC USE.
cTHEY ARE FURNISHED "AS IS." THE  AUTHORS, THE UNITED STATES GOVERNMENT,
cITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES,  AND AGENTS MAKE NO
cWARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS OF THE SOFTWARE
cAND DOCUMENTATION FOR ANY PURPOSE. THEY ASSUME NO RESPONSIBILITY
c(1) FOR THE USE OF THE SOFTWARE AND DOCUMENTATION; OR (2) TO PROVIDE
cTECHNICAL SUPPORT TO USERS.
c
cHISTORY
c August 8, 2017 Initial version prepared for Reference Environmental
c Data Record (REDR) program
c
c*****
C 2/21/07 reads in monthly anomalies 
C 9/8/04 heat_const revised (a0 droped)
C 4/13/04: added in-situ density climatology
C 4/8/4 HC computation with Cp(t,s,p) and Ro(t,s,p)
C
C
C OUT_VAR= OUT_UNITS * ODSQ_AREA * VERTICAL_INTEGRAL(FKTR*ANOMALY)
C
C where FKTR=1 for any OUT_VAR exept for STERIC COMPONENTS (SCs)
C For SCs, FKTR is thermal expansion or saline contraction coeff 
C which is a function of climatological ANNUAL temperature and salinity
C 
C ihalo:
C =0 (halosteric), 1 (thermosteric), 2 (heat content),
C  3 (freshwater), 4 (oxygen), 5 (whatever),
C  6 (layer average= vertical integral / thickness of layer)
C 
C itime:
C            climatological monthly mean,
C            1yr, 3mons, 6mons or 5yr anomalies
C            up to 2009 year
C
C iup(30) and ideep(30)
C           upper and lower depths of the layer
C           (up to 30 layers)
C
C idefo:
C       DEAFAULT CONSTANTS for OUTPUT 
C STER_CONST to convert TC or SC units from meters into mm
C
C HEAT_CONST and U18J to provide Heat Content in 10**18J
C Ro = 1.02 * 10**3 kg m**-3
C Cp = 4.187 * 10**3 J kg**-1 K**-1
C
C OXY_CONST to get output in 10**15 micro mol
C (conversion factor is 0.44642857 x 10 ** 5 micromol/m**3)
C
C FRESH_CONST to get freshwater in 10**9 cubic meter (that is cubic km)
C
C SALT_CONST=1. to get salt content in 1.e+10 kg
C salt cont = salt anom * density * volume, where input
C             salt anom [pss] or [1.e-3 kg of salt per kg of water]
C             density=(dens_in-situ+1000.)*10e-3 [1.e+3 kg/m**3]
C             volume=depth[m]*area[1.e+10 m**2]
C thus 1.e-3 x 1.e+3 x 1.e+10 =>  1.e+10 kg 
C
C OUT_UNITS is a user provided factor for NON DEAFAULT computation
C
C a0 = 1.237 * 10**10 m**2 area of 1x1 square at the equator
C
C NOTES:
C 2/4/04: adopted to work under KARA Linux
C 4/7/04: convertion of depth in METER to DBAR using OCL's function
C

      parameter (idim=360, jdim=180, kdim=102, kdimall=137)
      parameter (bmiss=-1.E10, nrec=idim*jdim)

      parameter (ster_const=1.e3)
      parameter (a0=1.237, heat_const=1.02*4.187, u18j=1.e-2)
      parameter (oxy_const=0.44642857)
      parameter (fresh_const=10.)
      parameter (salt_const=1.)
      parameter (pi=3.14159265359,pirat=pi/180.)
 
      dimension f(idim,jdim), fm(idim,jdim), fh(idim,jdim,kdim)
      dimension fkt(idim,jdim), fks(idim,jdim), fkd(idim,jdim)
      dimension fv(kdim), h(kdim),hall(kdimall)
      dimension iup(30),ideep(30)
      character*200 infile, infmean, junk
      character*200 outfile
      character*200 klimat, klimas, klimad, mess
      character*5 c6mon(23)
      character*4 cam(kdim)
      character*2 ctemp, c4
      character*1 cvar
      character*6 dotend, cstr(10)

C ODSQA is a function that returns an absolute value of cosine of LAT
      odsqa(alat)=abs(cos(alat))

      data c6mon/' ',
     *           '1-6 ','7-12',
     *           '1-3','4-6','7-9','10-12',
     *           ' ',' ',' ',' ',
     *   '1-1','2-2','3-3','4-4','5-5','6-6','7-7','8-8',
     *   '9-9','10-10','11-11','12-12'/
      data f/nrec*bmiss/

c     Get depths
      call setstandepsfromfile(numdeps,hall)
      if ( numdeps .gt. kdim ) numdeps=kdim
      do 35 n=1,numdeps
       h(n)=hall(n)
       cam(n)='XXXX'
       if ( h(n) .lt. 10. ) then
        cam(n)='X   '
       elseif ( h(n) .lt. 100. ) then
        cam(n)='XX  '
       elseif ( h(n) .lt. 1000. ) then
        cam(n)='XXX '
       endif
       inum=h(n)
       call insertnum('X',inum,cam(n),4)
35    continue

      open(17,file='./layer_field_10ann_climatology.inf',
     *     status='old',form='formatted')
      read(17,'(a)') junk
      read(17,'(a)') junk
      read(17,*) itime1

      itime=itime1-1

      read(17,'(a)') junk
      read(17,'(a)') junk
      read(17,*) ihalo
      read(17,'(a)') junk
      read(17,*) ihctsp
      read(17,'(a)') junk
      read(17,*) km2m
      read(17,'(a)') junk
      read(17,*) idefo
      read(17,'(a)') junk
      read(17,*) out_units
      read(17,'(a)') junk
      read(17,*) iarea
      read(17,'(a)') junk
      read(17,'(a)') infile
      read(17,'(a)') junk
      read(17,'(a6)')      dotend
      read(17,'(a)') junk
      read(17,*) iflag
      read(17,'(a)') junk
      read(17,'(a)') infmean
      read(17,'(a)') junk
      read(17,'(a)') outfile
      read(17,'(a)') junk
      read(17,*) ifst,ilst
      read(17,'(a)') junk
      read(17,*) (iup(i),ideep(i), i=1,10)
      read(17,*) (iup(i),ideep(i), i=11,20)
      read(17,*) (iup(i),ideep(i), i=21,30)
      read(17,'(a)') junk
      read(17,'(a)') klimat
      read(17,'(a)') klimas
      read(17,'(a)') klimad
      read(17,'(a)') junk
      read(17,*) meter2dbar
      read(17,'(a)') junk
      read(17,*) saltmin

      close(17)
      write(cvar,'(i1)') iflag

C to use default constant for output or not
      if(idefo.eq.1)then
         if(ihalo.le.1)then
C steric components
            out_units=ster_const
            iarea=0
         elseif(ihalo.eq.2)then
C heat content
            out_units=heat_const * u18j
            iarea=1
         elseif(ihalo.eq.3)then
C freshwater
            out_units=fresh_const
            iarea=1
         elseif(ihalo.eq.4)then
C oxy
            out_units=oxy_const
            iarea=1
         elseif(ihalo.eq.5)then
C salt content 
            out_units=salt_const
            iarea=1
         elseif(ihalo.eq.6)then
c average anomaly of layer
            out_units=1.
            iarea=0
         endif
      endif

C Open Climatological Temperature, Salinit, and Density
       call addCend(klimat,200,0)
       call fileopen(klimat,'rb+'//CHAR(0),15)
       call addCend(klimas,200,0)
       call fileopen(klimas,'rb+'//CHAR(0),16)
       call addCend(klimad,200,0)
       call fileopen(klimad,'rb+'//CHAR(0),19)

C Open grand mean anomaly file
       if(iflag.eq.0)then
       call addCend(infmean,200,0)
       call fileopen(infmean,'rb+'//CHAR(0),11)
       endif
 
       mclimate=0
C initialize time loops
       if(itime.eq.5) then
C        MONTHLY CLIMATOLOGY
         mclimate=1
         ilast=ilst
         naTE=0
         itime1=1
         i6mod=0
       elseif(itime.eq.4) then
C        PENTADAL ANOMALIES
         ilast=ilst-4
         naTE=4
         itime1=1
         i6mod=0
       else
C        YEARLY, 6mons, or 3mons ANOMALIES
         ilast=ilst
         naTE=0
         i6mod=itime
       endif

C time LOOP
 
      do 2000 infl=ifst, ilast

      call nameTIME(naTE,infl,ctemp,c4)

c "subtime" loop

      do 2000 i6=1,itime1

      if(mclimate.eq.0)then
       cstr(1)=ctemp
       cstr(2)=c4
       cstr(3)=c6mon(i6+i6mod)
c      cstr(4)='00AM'
       cstr(4)=dotend
       call janameit(4,infile,'$',mess,cstr)
      else
       cstr(1)=ctemp
       cstr(2)=dotend
       call janameit(2,infile,'$',mess,cstr)
      endif
       call addCend(mess,200,0)
       call fileopen(mess,'rb+'//CHAR(0),10)
 
C Start reading data from a file LOOP 
 
      do 1999 itr=1,30

      n1=iup(itr)
      n2=ideep(itr)
      if (n1.eq.0) then
          call onefileclose(10)
          goto 2000
      endif

       cstr(1)=cam(n1)
       cstr(2)='-'
       cstr(3)=cam(n2)
       cstr(4)='-'
       cstr(5)=ctemp
       cstr(6)=c4
       cstr(7)=c6mon(i6+i6mod)
       cstr(8)='.'
       cstr(9)=cvar
       call janameit(9,outfile,'$',mess,cstr)
       call addCend(mess,200,0)
       call fileopen(mess,'rb+'//CHAR(0),12)

      do 999 k=n1,n2

       call readlevel(1,1,1,1,1,idim,jdim,idim,jdim,idim,jdim,
     *                k,k,10,f)
       if(iflag.eq.0)then
          call readlevel(1,1,1,1,1,idim,jdim,idim,jdim,idim,jdim,
     *                   k,k,11,fm)
       else
            call rinit2D(fm,idim,jdim,0.)
       endif
       call readlevel(1,1,1,1,1,idim,jdim,idim,jdim,idim,jdim,
     *                k,k,15,fkt)
       call readlevel(1,1,1,1,1,idim,jdim,idim,jdim,idim,jdim,
     *                k,k,16,fks)
       call readlevel(1,1,1,1,1,idim,jdim,idim,jdim,idim,jdim,
     *                k,k,19,fkd)

       do i=1,idim
       do j=1,jdim

	fval=f(i,j)
        fnrm=fm(i,j)
        saltref=fks(i,j)
        densitu=fkd(i,j)

C for freshwater
C apply restriction for minimum climatological salinity
C UPDATE below to use with Steric and HC based on NON-CONSTANT density
        if(ihalo.eq.3)then
        if(saltref.gt.bmiss)then
        if(saltref.lt.saltmin)then
           saltref=saltmin
        endif
        endif
        endif

	if((fval.le.bmiss) .or. (fnrm.le.bmiss)) then
            fh(i,j,k)=bmiss
	else

           x1=fkt(i,j)
           x2=saltref
           x3degree=j-90.5
           x3= x3degree * pirat

           if(meter2dbar.eq.1)then
c convert METER to DBAR using
c OCL's function
              if(k.eq.1)then
                 x4=0.
              else
                 pressPrev=fpress(h(k-1),x3degree,0.,0.,3)
                 x4=fpress(h(k),x3degree,h(k-1),pressPrev,3)
              endif
           else
c or Leroy(1969) I used before April 7, 2004
                 x4=deplat(h(k),x3)
           endif

         if(ihalo.eq.0)then
           call thermik(x1,x2,x4,ralf,rbet)
           rkof=-1.0*rbet
         elseif(ihalo.eq.1)then
           call thermik(x1,x2,x4,ralf,rbet)
           rkof=ralf
         elseif(ihalo.eq.3)then
           densSALT=fdense(x1,x2,x4)+1000.
           densFRESH=fdense(x1,0.,x4)+1000.
           rkof=(densSALT*(fval-fnrm))/(densFRESH*(x2+fval-fnrm))
           rkof=-1.*rkof
           fval=1.
           fnrm=0.
C fresh output in METER or CUBIC KM
           if(km2m.lt.1)then
              iarea=0
              out_units=1.
           endif
         elseif(ihalo.eq.2 .and. ihctsp.eq.1)then
           dTSP  = 1.e-3 * (fdense(x1,x2,x4)+1000.)
           CpTSP = 1.e-3 * cpsw(x2,x1,x4)
           rkof=dTSP*CpTSP
           out_units=u18j
         elseif(ihalo.eq.5) then
           rkof= 1.e-3 * (densitu+1000.)
         else
           rkof=1.
         endif

	   fh(i,j,k)= rkof * (fval-fnrm)
           if ( i .eq. 180 .and. j .eq. 90 ) 
     *     write(6,*) 'FH',k,rkof,fval,fnrm

	endif

       enddo
       enddo
 999  continue
 
C Compute vertical integral or average anomaly 
C of current ITR-th layer and
C write it as the 1st record into output file
 
      do 900 j=1,jdim

       glat=(j-90.5)*pirat
       if(iarea.eq.1)then
          a1x1sq = a0 * odsqa(glat)
       else
          a1x1sq = 1.
       endif

       do 900 i=1,idim
        do k=n1,n2
           fv(k)=fh(i,j,k)
        enddo

        call boxintegral01(fv,n1,n2,bmiss,fout,fw,h,kdim)

        if(fout.gt.bmiss)then
         if(ihalo.lt.6)then
           f(i,j) = out_units * a1x1sq * fout
           if ( i .eq. 180 .and. j .eq. 90 ) 
     * write(6,*) 'LAYER',out_units,a1x1sq,fout
         else
           f(i,j) = out_units * a1x1sq * fout/fw
         endif
        else
           f(i,j) = bmiss
        endif

  900 continue
  
      call writelevel(1,1,1,1,1,idim,jdim,idim,jdim,idim,jdim,
     *                1,1,12,f)
      call onefileclose(12)

 1999 continue

      call onefileclose(10)

 2000 continue
 
      if(iflag.eq.0) call onefileclose(11)

      stop
      end
c     include 'sub_Cp.f'
c     include './antonov/sub_alphabeta.f'
c     include './antonov/sub_nameTIME.f'
c     include './antonov/sub_fullname.f'
c     include './antonov/sub_init_array.f'
c     include './antonov/sub_depth_integral_tmod.f'