program layer_field_SEK05
c@***h* OCEAN_HEAT_CONTENT/layer_field_SEK05.f
c
c NAME
c layer_field_SEK05.f
c
c PURPOSE
c Specifically for OCEAN_HEAT_CONTENT (OHC), OHC standard
c error in a given layer of the ocean is calculated using
c standard error of temperature anomaly (TA) objectively analyzed
c fields (calculated in analysis.f).
c
c DESCRIPTION
c layer_field_SEK05.f is a multipurpose program which calculates an
c standard error of an ocean variable over a volume of water. 
c A gridded field is read in for each standard depth level within the
c designated depth range, the standard error of the ocean variable
c is calculated 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_SEK05.f comes from
c 'layer-field-10ann_climatology.inf', including designated
c standard depth layers, time period, temporal span (monthly,
c quarterly, yearly, pentadal), input and output file forms,
c latitude/longitude span.
c
c OUTPUT
c Files of standard error OHC for user designated layer in
c units of 10^18 joules. Files are in the World Ocean Atlas
c (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 04/28/10 array c6mon is updated to handle monthly anomalies
C 01/17/06
C   
C 9/8/04 heat_const fixed (a0 droped)
C 8/11/04 BASED on layer-fieldK04
C         !!! doesn't work for FRESHWATER
C
C OUT_SE= OUT_UNITS * ODSQ_AREA * SQRT(VI), where
C
C VI is a vertical integral of (FKTR*SD*W*H)**2
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 For Heat: could be CpRo fixed or f(T,S,P)
C H is a thickness of elementary sublayer
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=33)
      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 (cp0=3991.868)
      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 fanom(idim,jdim)
      dimension fv(33), h(33)
      dimension iup(30),ideep(30)
      character*80 infile, infmean, outfile, junk
      character*80 klimat, klimas, klimad, insanom,  mess
      character*5 c6mon(23)
      character*4 cam(33)
      character*2 ctemp, c4
      character*1 cvar
      character*6 cstr(10)

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

      data h/0.,10.,20.,30.,50.,75.,100.,125.,150.,200.,250.,
     *      300.,400.,500.,600.,700.,800.,900.,1000.,1100.,
     *      1200.,1300.,1400.,1500.,1750.,2000.,2500.,3000.,
     *      3500.,4000.,4500.,5000.,5500./
      data cam/'0   ','10  ','20  ','30  ','50  ','75  ','100 ',
     *         '125 ','150 ','200 ','250 ','300 ','400 ','500 ',
     *         '600 ','700 ','800 ','900 ','1000','1100','1200',
     *         '1300','1400','1500','1750','2000','2500','3000',
     *         '3500','4000','4500','5000','5500'/
      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/

      open(17,file='./layer_fieldSE_K05.inf',
     *     status='old',form='formatted')

      read(17,'(a80)') junk
      read(17,'(a80)') junk
      read(17,*) itime1

      itime=itime1-1

      read(17,'(a80)') junk
      read(17,'(a80)') junk
      read(17,*) ihalo
      read(17,'(a80)') junk
      read(17,*) ihctsp
      read(17,'(a80)') junk
      read(17,*) km2m
      read(17,'(a80)') junk
      read(17,'(a80)') insanom
      read(17,'(a80)') junk
      read(17,*) idefo
      read(17,'(a80)') junk
      read(17,*) out_units
      read(17,'(a80)') junk
      read(17,*) iarea
      read(17,'(a80)') junk
      read(17,'(a80)') infile
      read(17,'(a80)') junk
      read(17,*) iflag
      read(17,'(a80)') junk
      read(17,*) idodi
      read(17,'(a80)') junk
      read(17,'(a80)') infmean
      read(17,'(a80)') junk
      read(17,'(a80)') outfile
      read(17,'(a80)') junk
      read(17,*) ifst,ilst
      read(17,'(a80)') 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,'(a80)') junk
      read(17,'(a80)') klimat
      read(17,'(a80)') klimas
      read(17,'(a80)') klimad
      read(17,'(a80)') junk
      read(17,*) meter2dbar
      read(17,'(a80)') junk
      read(17,*) saltmin

      close(17)
      write(cvar,'(i1)') iflag
      if(iflag .eq. 1) cvar='n'
      if(iflag .eq. 3) cvar=' '

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,80,0)
       call fileopen(klimat,'rb+'//CHAR(0),15)
       call addCend(klimas,80,0)
       call fileopen(klimas,'rb+'//CHAR(0),16)
       call addCend(klimad,80,0)
       call fileopen(klimad,'rb+'//CHAR(0),19)

       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

C open .SD input file
      if(mclimate.eq.0)then
       cstr(1)=ctemp
       cstr(2)=c4
       cstr(3)=c6mon(i6+i6mod)
       cstr(4)='00AM'
       cstr(5)='.s' // cvar
       call janameit(5,infile,'$',mess,cstr)
      else
       cstr(1)=ctemp
       cstr(2)='.s' // cvar
       call janameit(2,infile,'$',mess,cstr)
      endif
       call addCend(mess,80,0)
       call fileopen(mess,'rb+'//CHAR(0),10)

C open .POINTS input file if required
      if(idodi.eq.1)then
      if(mclimate.eq.0)then
       cstr(1)=ctemp
       cstr(2)=c4
       cstr(3)=c6mon(i6+i6mod)
       cstr(4)='00AM.'
       cstr(5)='points'
       call janameit(5,infmean,'$',mess,cstr)
      else
       cstr(1)=ctemp // '.'
       cstr(2)='points'
       call janameit(2,infmean,'$',mess,cstr)
      endif
       call addCend(mess,80,0)
       call fileopen(mess,'rb+'//CHAR(0),11)
      endif

C open Salinity Anomaly input file (only for freshwater)
      if(ihalo.eq.3)then
      if(mclimate.eq.0)then
       cstr(1)=ctemp
       cstr(2)=c4
       cstr(3)=c6mon(i6+i6mod)
       cstr(4)='00AM'
       call janameit(4,insanom,'$',mess,cstr)
      else
       cstr(1)=ctemp
       call janameit(1,insanom,'$',mess,cstr)
      endif
       call addCend(mess,80,0)
       call fileopen(mess,'rb+'//CHAR(0),17)
      endif

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)
          if(idodi.eq.1) call onefileclose(11)
          if(ihalo.eq.3) call onefileclose(17)
          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
       cstr(10)='se'
       call janameit(10,outfile,'$',mess,cstr)
       call addCend(mess,80,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(idodi.eq.1)then
          call readlevel(1,1,1,1,1,idim,jdim,idim,jdim,idim,jdim,
     *                   k,k,11,fm)
       else
            call rinit2D(fm,idim,jdim,1.)
       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)
       if(ihalo.eq.3)then
          call readlevel(1,1,1,1,1,idim,jdim,idim,jdim,idim,jdim,
     *                   k,k,17,fanom)
       else
            call rinit2D(fanom,idim,jdim,1.)
       endif

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

	fval=f(i,j)
        fnrm=fm(i,j)

        sanom=fanom(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. sanom.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=-1.*densSALT/densFRESH
           rkof= rkof * x2 / (sanom+x2)**2.
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 / sqrt(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 boxintegralSE(fv,n1,n2,bmiss,fout,fw)

        if(fout.gt.bmiss)then
         if(ihalo.lt.6)then
           f(i,j) = 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)
      if(idodi.eq.1) call onefileclose(11)
      if(ihalo.eq.3) call onefileclose(17)

 2000 continue
 

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