block data blkdat
c
c=======================================================================
c                                                                    ===
c  Sub-program to initialize some variables and arrays in named      ===
c  common blocks when the ANALYTICAL option is activated.            ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <param.h>
#include <pconst.h>
#include <scalar.h>
#include <fullwd.h>
#include <filtdat.h>
#include <rhomean.h>
#include <cvmix.h>
#include <runid.h>
#include <version.h>
#ifdef analytical
# include <hybrid.h>
# include <moddat.h>
# include <options.h>
#endif
#ifdef ldrifters
# include <tracks.h>
#endif
#include <pe_netcdf.h>
#include <iounits.h>
#ifdef peprf
# include <hydro.h>
# include <prf.h>
#endif
#if defined oias & defined fcsterr
# include <ioi.h>
#endif
#if defined bioMcGillic | defined bioFasham | defined bioAnder | defined bioDuse
# include <cbiopnh.h>
#endif
#if defined analytical | defined peprf | (defined oias & defined fcsterr ) 
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
#endif
#ifdef analytical
      integer profpts
      parameter (profpts=mprof*3)
#endif
#ifdef peprf
      integer prfpts
      parameter (prfpts=(mnbprf+1)*4)
#endif
#if defined oias & defined fcsterr
      integer trcpts,volpts
      parameter (volpts=imtjmt*km, trcpts=volpts*nt)
#endif
c
c=======================================================================
c  Initialize some common block data.
c=======================================================================
c
c-----------------------------------------------------------------------
c  Set the model name, closing date of current version & version number.
c-----------------------------------------------------------------------
c
      data model /'PE MODEL  '/
      data vdate /'08/21/2001'/
      data vnum  /'9.6       '/
c
c-----------------------------------------------------------------------
c  Set the type of mixing timestep.  If EB=.TRUE., an Euler
c  backward step is done; if EB=.FALSE., a forward step is done.
c-----------------------------------------------------------------------
c
      data eb /.false./
c
c-----------------------------------------------------------------------
c  Initialize various quantities.
c-----------------------------------------------------------------------
c
      data ifprnt,nwrite /0,1000000/
      data ntsi,ncon /1,1/
      data na,nb,nc /0,0,0/
      data omega /earthrot/
      data radius /re/
      data radian /rad2deg/
      data grav /980.6/
      data icntv,icntt,icntz,icntp /0,0,0,0/
c
c-----------------------------------------------------------------------
c  Set I/O units.  Units 101 & 103 are reserved for nested runs.
c-----------------------------------------------------------------------
c
      data stdinp,stdout /5,6/
      data kontrl,kflds /11,12/
      data labs /13,14,15/
      data trkinp /34/
      data ioiout,ioipar /35,81/
      data bioinp /43/
#ifdef bioDuse
      data absinp /44/
#endif
      data tsrcinp /45/
      data inpprf,outprf /47,48/
      data usrinp /50/
      data xtrinp /60/
      data tdbxin /90/
#ifdef ext_tide
      data tdphin /97/
#endif
c
c-----------------------------------------------------------------------
c  Initialize netCDF output states.
c-----------------------------------------------------------------------
c
      data ncnrgst,ncoutst,nctrkst /3*0/
#ifdef pressbias
      data ncpbst /0/
#endif 
c
c-----------------------------------------------------------------------
c  Parameters associated with equation of state.
c-----------------------------------------------------------------------
c
c  Mean salinity value subtracted from salinity to improve accuracy.
c
      data smean /35.0/
c
c  Mean density (kg/m^3).
c
      data rho0 /1025.0/
c
c  Mean temperature (TBAR) and salinity (SBAR) used in linear equation
c  of state and their associate mean density (RBAR).
c
      data tbar /10.0/
      data sbar /35.0/
      data rbar /1026.9524/
c
c  Thermal expansion and saline contraction coefficients.
c
      data tcoef /1.67e-04/
      data scoef /7.62e-04/
c
#ifdef ldrifters
c-----------------------------------------------------------------------
c  Parameters associated with Lagrangian drifters.
c-----------------------------------------------------------------------
c
      data nfloats /0/
c
#endif
#ifdef peprf
c-----------------------------------------------------------------------
c  Parameters associated with profile sampling.
c-----------------------------------------------------------------------
c
      data cm_bgn,ctd_bgn,xbt_bgn /3*0/
      data cm_emt,ctd_emt,xbt_emt /3*1/
      data nbcm,nbctd,nbxbt /3*0/
c
      data icmpos  /prfpts*0/
      data ictdpos /prfpts*0/
      data ixbtpos /prfpts*0/
c
      data castid,hflag,nhpts,jdayoff /0,0,kmp1,2440000/
c
      data axbt,cm,ctd,eof,xbt,xbts,xctd /7*.false./
c
#endif
c-----------------------------------------------------------------------
c  Parameters associated with vertical mixing.
c-----------------------------------------------------------------------
c
c  Specific heat (ergs/gram/^C).  It is approximately 4.0e+7, and
c  varies only slightly (see Gill, 1982, Appendix 3).
c
      data cp /4.0e+7/
c
c  Minimum wind stress (dynes/cm^2) used in computation of mixed-layer
c  depth:
c          TAUMIN = RHOair  * Cd * Wmin^2
c  where (Pond and Picard, p 108)
c          RHOair = 1.3 kg/m^3 ; Cd = 1.4e-3
c          TAUMIN = 1.64e-3 dynes/cm^2 (using Wmin = 0.3 m/s)
c
      data taumin /1.64e-3/
c
c  Threshold latitude used to compute Coriolis (F0).  Since Coriolis is
c  zero at the Equator and there is a division by F0 in the algorithm
c  for Ekman depth, this latitude value is used in the computation for
c  Ekman depth at latitudes less than LATMF0.
c
      data latmf0 /10.0/
c
c  Wind stress coefficient MCOEF for estimating mixed-layer depth. It is
c  the "m" parameter from Eq. 10.30 in Niiler and Kraus (1977; in:
c  Modeling and Prediction of the Upper Layers of the Ocean, Ed by E.B.
c  Kraus, p 143-172).  The larger MCOEF, the larger the influence of
c  wind stress on  mixed-layer depth (relative to buoyancy). It is
c  generally observed that MCOEF=1.25, although it can be > 8  during
c  storm events (ibid.). Actually, MCOEF is not constant but increases
c  as the Richardson Number decreases (ibid). Presuming storm events of
c  MCOEF=8 to occur 1/7 or 1/30 of the time would suggest a time-mean
c  MCOEF of 2.21 or 1.48, respectively. Thus, generally use
c  1.25 =< MCOEF =< 2.21.
c
      data mcoef /1.25/
c
c  Buoyancy coefficient NCOEF for estimating mixed-layer depth. It is
c  the "n" parameter from Eq. 10.30 in Niiler and Kraus (1977). It is a
c  weight for the effects of positive buoyancy; 0 =< NCOEF =< 1:
c      NCOEF=0  all convectively produced energy is dissipated i.e.
c               a positive buoyancy flux (surface cooling) does not
c               influence MLD.
c      NCOEF=1  convectively-produced eddies are not dissipated i.e.
c               positive buoyancy influences mixed-layer depth in an
c               equal and opposite manner as negative buoyancy).
c  Observations suggest 0.0<NCOEF<0.113, and perhaps typically
c  NCOEF=0.02 or 0.03 (ibid.). Actually, NCOEF is not constant but
c  decreases as mixed-layer depth increases (ibid.).  However, the
c  original algorithm is actually poorly defined (i.e., can blow up)
c  unless NCOEF=0 ! So use NCOEF=0 ! (Gill and Turner (1976) also
c  support NCOEF=0).
c
      data ncoef /c0/
c
c-----------------------------------------------------------------------
c  Parameters associated with irradiance.
c-----------------------------------------------------------------------
c
c  Light extinction coefficient (1/cm). Parsons et al. (1984,p95) say
c  ATTH2O=3.84e-4 for clear ocean water (Jerlov,1968), but round up due
c  to POM and DOM to 4.0e-4  1/cm (about 100m = 2% light level).
c  Note: coastal water can be much thicker! (see Parsons et al., p 70).
c
      data atth2o /4.0e-4/
c
#ifdef bioDuse
c  Lower and upper limits of the spectral range for use in productivity
c  models.
c
      data lmdamin /400.0/
      data lmdamax /700.0/
c
c  Conversion factor for watts to quanta per sec, for 400-700 nm
c  range.  From Morel and Smith (1974), units are quanta/sec/watt,
c  and accuracy is plus or minus a few percent for most meteorological
c  conditions.
c
      data ht2qnta /2.77e+18/
c
c  Spectral attenuation coefficients (m^-1) from Morel and Antoine,
c  1994 (JPO, 24:1652-1665) and Morel 1988 (...)
c  the 400 to 700 nm range only is used here at this time
c  in 5 nm increments
c
      data atth2os /0.021,0.020,0.020,0.019,0.018,0.018,0.017,0.017,
     *              0.017,0.017,0.017,0.017,0.017,0.017,0.018,0.018,
     *              0.019,0.020,0.022,0.024,0.027,0.032,0.038,0.045,
     *              0.049,0.051,0.052,0.054,0.057,0.062,0.064,0.068,
     *              0.072,0.076,0.081,0.094,0.107,0.128,0.157,0.200,
     *              0.253,0.279,0.296,0.303,0.310,0.315,0.320,0.325,
     *              0.330,0.340,0.350,0.370,0.405,0.418,0.430,0.440,
     *              0.450,0.470,0.500,0.550,0.650/
c
c  Chi(lambda) values also from Morel and Antoine, 1994.
c  Same range as above.  Used with e(lambda) values below
c  to estimate attenuation due to chlorophyll (phytoplankton).
c
      data chil /0.110,0.111,0.112,0.113,0.113,0.110,0.108,0.106,
     *           0.104,0.100,0.097,0.094,0.090,0.086,0.082,0.079,
     *           0.075,0.073,0.069,0.066,0.064,0.060,0.058,0.054,
     *           0.050,0.047,0.047,0.045,0.044,0.043,0.041,0.040,
     *           0.039,0.038,0.036,0.034,0.033,0.033,0.032,0.033,
     *           0.034,0.035,0.036,0.038,0.038,0.040,0.042,0.043,
     *           0.044,0.045,0.045,0.046,0.047,0.049,0.052,0.052,
     *           0.051,0.044,0.039,0.034,0.030/
c
c  e(lambda) values, paired with above chi values.
c
      data elam /0.668,0.672,0.680,0.687,0.693,0.701,0.707,0.708,
     *           0.707,0.704,0.701,0.699,0.700,0.703,0.703,0.703,
     *           0.703,0.704,0.702,0.700,0.700,0.695,0.690,0.685,
     *           0.680,0.675,0.670,0.665,0.660,0.655,0.650,0.645,
     *           0.640,0.630,0.623,0.615,0.610,0.614,0.618,0.622,
     *           0.626,0.630,0.634,0.638,0.642,0.647,0.653,0.658,
     *           0.663,0.667,0.672,0.677,0.682,0.687,0.695,0.697,
     *           0.693,0.665,0.640,0.620,0.600/
c
c  wsj(lambda) values, taken from Anderson, 1993 (L&O,38:1403-1419).
c
      data wsj  /0.0053,0.0103,0.0108,0.0113,0.0014,0.0014,0.0107,
     *           0.0126,0.0127,0.0140,0.0154,0.0154,0.0157,0.0156,
     *           0.0158,0.0163,0.0169,0.0156,0.0161,0.0169,0.0164,
     *           0.0170,0.0171,0.0163,0.0164,0.0172,0.0178,0.0177,
     *           0.0174,0.0180,0.0180,0.0181,0.0179,0.0182,0.0180,
     *           0.0185,0.0186,0.0186,0.0169,0.0178,0.0181,0.0187,
     *           0.0187,0.0186,0.0191,0.0188,0.0187,0.0190,0.0190,
     *           0.0189,0.0183,0.0173,0.0182,0.0194,0.0191,0.0189,
     *           0.0188,0.0180,0.0164,0.0170,0.0089/
c
c     photoacclimation rate for chl:c (?) ratio
c
      data chlcrate /3.89e-5/
#endif
c
c  Initialize run titles.
c
      data titlmdl /'pe c','ode:','    ','    ','filt','er:(','    ',
     *              '    ',')   ','    ','    ','    ','    ','    ',
     *              '    ','    ','    ','    ','    ','    '/
      data titlebc /'bdy ','cond',':   ','    ','    ','    ','    ',
     *              '    ','    ','    ','    ','    ','    ','    ',
     *              '    ','    ','    ','    ','    ','    '/
c
c  Set names of tracer variables.
c
      data tname(1,1) /'temp                                        '/
      data tname(2,1) /'temperature                                 '/
      data tname(3,1) /'Celsius                                     '/
      data tname(4,1) /'temperature, scalar, series                 '/
      data tname(1,2) /'salt                                        '/
      data tname(2,2) /'salinity                                    '/
      data tname(3,2) /'PSU                                         '/
      data tname(4,2) /'salinity, scalar, series                    '/
#if defined bioAnder | defined bioFasham | defined bioMcGillic
      data tname(1,ino3)/'NO3                                         '/
      data tname(2,ino3)/'nitrate concentration                       '/
      data tname(3,ino3)/'millimoles nitrogen meter-3                 '/
      data tname(4,ino3)/'nitrate, scalar, series                     '/
      data tname(1,iphy)/'phyto                                       '/
      data tname(2,iphy)/'phytoplankton concentration                 '/
      data tname(3,iphy)/'millimoles nitrogen meter-3                 '/
      data tname(4,iphy)/'phytoplankton, scalar, series               '/
      data tname(1,izoo)/'zoo                                         '/
      data tname(2,izoo)/'zooplankton concentration                   '/
      data tname(3,izoo)/'millimoles nitrogen meter-3                 '/
      data tname(4,izoo)/'zooplankton, scalar, series                 '/
      data tname(1,inh4)/'NH4                                         '/
      data tname(2,inh4)/'ammonium concentration                      '/
      data tname(3,inh4)/'millimoles nitrogen meter-3                 '/
      data tname(4,inh4)/'ammonium, scalar, series                    '/
# if defined bioAnder
      data tname(1,idet)/'detritus                                    '/
      data tname(2,idet)/'detritus concentration                      '/
      data tname(3,idet)/'millimoles nitrogen meter-3                 '/
      data tname(4,idet)/'detritus, scalar, series                    '/
# elif defined bioFasham
      data tname(1,ipon)/'PON                                         '/
      data tname(2,ipon)/'particulate organic nitrogen concentration  '/
      data tname(3,ipon)/'millimoles nitrogen meter-3                 '/
      data tname(4,ipon)/'PON, scalar, series                         '/
      data tname(1,idon)/'DON                                         '/
      data tname(2,idon)/'dissolved organic nitrogen concentration    '/
      data tname(3,idon)/'millimoles nitrogen meter-3                 '/
      data tname(4,idon)/'DON, scalar, series                         '/
      data tname(1,ibac)/'bacteria                                    '/
      data tname(2,ibac)/'bacteria concentration                      '/
      data tname(3,ibac)/'millimoles nitrogen meter-3                 '/
      data tname(4,ibac)/'bacteria, scalar, series                    '/
# endif
#endif
#ifdef bioDuse
      data tname(1,izoo)/'zoo                                         '/
      data tname(2,izoo)/'zooplankton concentration                   '/
      data tname(3,izoo)/'millimoles nitrogen meter-3                 '/
      data tname(4,izoo)/'zooplankton, scalar, series                 '/
      data tname(1,ino3)/'NO3                                         '/
      data tname(2,ino3)/'nitrate concentration                       '/
      data tname(3,ino3)/'millimoles nitrogen meter-3                 '/
      data tname(4,ino3)/'nitrate, scalar, series                     '/
      data tname(1,inh4)/'NH4                                         '/
      data tname(2,inh4)/'ammonium concentration                      '/
      data tname(3,inh4)/'millimoles nitrogen meter-3                 '/
      data tname(4,inh4)/'ammonium, scalar, series                    '/ 
      data tname(1,idet)/'detritus                                    '/
      data tname(2,idet)/'detritus concentration                      '/
      data tname(3,idet)/'millimoles nitrogen meter-3                 '/
      data tname(4,idet)/'detritus, scalar, series                    '/  
      data tname(1,iqn3)/'CELLNO3                                     '/
      data tname(2,iqn3)/'cellular nitrogen (from nitrate)            '/
      data tname(3,iqn3)/'millimoles nitrogen meter-3                 '/
      data tname(4,iqn3)/'cellnitrate, scalar, series                 '/
      data tname(1,iqn4)/'CELLNH4                                     '/
      data tname(2,iqn4)/'cellular nitrogen (from ammonia)            '/
      data tname(3,iqn4)/'millimoles nitrogen meter-3                 '/
      data tname(4,iqn4)/'cellammonia, scalar, series                 '/
      data tname(1,ichl)/'CHL                                         '/
      data tname(2,ichl)/'chlorophyll concentration                   '/
      data tname(3,ichl)/'milligrams meter-3                          '/
      data tname(4,ichl)/'chlorophyll, scalar, series                 '/
#endif
#if defined codunlim | defined codlim
      data tname(1,icod)/'cod                                         '/
      data tname(2,icod)/'cod concentration                           '/
      data tname(3,icod)/'milligrams meter-3                          '/
      data tname(4,icod)/'cod, scalar, series                         '/
#endif
#ifdef analytical
c
c-----------------------------------------------------------------------
c  Set up internal parameters during Analytical initialization.
c-----------------------------------------------------------------------
c
c  Set run identification labels.
c
      data runic   / 'ANAF','LDS:',' IFF','    ','18/0','9/95' /
      data titleic / 'Anal','ytic',' Ini','tial','izat','ion:','    ',
     *               '    ','    ','.   ','    ','    ','    ','    ',
     *               '    ','    ','    ','    ','    ','    ' /
c
c  Set initialization flags.
c
c  IFLAG(1)   diagnostic initialization switch: [0] no, [1] yes
c  IFLAG(2)   interpolation switch: [0] linear, [1] cubic-spline
c  IFLAG(3)   persistent boundary conditions switch: [0] no, [1] yes
c  IFLAG(4)   boundary condition on PSI: [1] PSI, [2] d(vort)/dt
c  IFLAG(5)   not used
c  IFLAG(6)   used in preprocessing program "PE_initial.f"
c  IFLAG(7)   used in preprocessing program "PE_initial.f"
c  IFLAG(8)   hybrid coordinates interface level, KC
c  IFLAG(9)   hybrid coordinates transformation: [0] single, [1] double
c  IFLAG(10)  coastal setup switch: [0] no, [1] coast(s) only,
c                                           [2] coast(s) and island(s)
c
c                  1  2  3  4  5  6  7  8  9 10
# ifdef dblsigma
      data iflag / 0, 0, 0, 2, 0, 1, 1,12, 0, 0 /
# else
      data iflag / 0, 0, 0, 2, 0, 1, 1, 1, 0, 0 /
# endif
c
c  Set domain configuration parameters:
c     COORD  - flag for horizontal coordinate
c                 [0] cartesian (tangent plane)
c                 [1] unrotated spherical
c                 [2] rotated spherical
c     RLATD  - transformation centroid latitude (degrees; North is positive)
c     RLNGD  - transformation centroid longitude  (degrees; WEST is negative)
c     DELX   - X-offset between model and transform centers.  (cm or degrees) 
c     DELX   - Y-offset between model and transform centers.  (cm or degrees)
c     THETAD - domain rotation angle  (degrees)
c     GRIDX  - model grid spacing in the x-direction.  (cm or degrees) 
c     GRIDY  - model grid spacing in the y-direction. (cm or degrees)
c
      data coord  /0/
      data rlatd  /64.25/
      data rlngd  /-11.5/
      data delx   /c0/
      data dely   /c0/
      data thetad /c0/
      data gridx  /500000.0/
      data gridy  /500000.0/
c
c  Set double hybrid vertical transformation parameters: minimum
c  depth (cm) of coordinate interface ZC1, maximum depth (cm) of
c  coordinate interface ZC2, reference depth (cm) for coordinate
c  interface ZREF, and slope parameter (positive) of the coordinate
c  interface ZSLOPE.
c
# ifdef dblsigma
      data zc1  / 5000.0/
      data zc2  /15000.0/
      data zref /35000.0/
      data zslope /0.5/
# else
      data zc1  /c0/
      data zc2  /c0/
      data zref /c0/
      data zslope /c0/
# endif
c
c  Set KM flat level depths REFZ (cm) at the center of the vertical box.
c
      data refz / 100.0, 226.5, 288.5, 372.5, 486.0, 639.0, 846.0,
     &           1125.5,1503.0,2013.0,2702.0 /
c
c  Set in mean profile (T,S,Z).  Two extra set of values need to be
c  provided in addition to the model number of levels KM.  This set of
c  values are needed at the top and bottom of the profile in order to
c  bound the interpolation in routine MEANRHO.  These values may or
c  may not be used in ANAMRHO.
c
      data tinit /profpts*c0/
      data nprof /mprof/
#endif
#if defined oias & defined fcsterr
c
c  Initialize forecast error fields to "no data".
c
      data ufcterr /volpts*c1/
      data vfcterr /volpts*c1/
      data pfcterr /imtjmt*c1/
      data tfcterr /trcpts*c1/
#endif
      end