Back to home page

MITgcm

 
 

    


File indexing completed on 2023-02-07 06:09:54 UTC

view on githubraw file Latest commit 2b959ba3 on 2023-02-06 20:20:10 UTC
3996ba3d32 Jean*0001 #include "ECCO_OPTIONS.h"
2b959ba38e Mart*0002 #ifdef ALLOW_SEAICE
                0003 # include "SEAICE_OPTIONS.h"
                0004 #endif
                0005 
                0006 C--  File cost_gencost_seaicev4.F:
                0007 C--   Contents
                0008 C--   o COST_GENCOST_SEAICEV4
                0009 C--   o GET_EXCONC_DECONC
                0010 
                0011 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
2ac072a19d Gael*0012 
f8e779c983 antn*0013       subroutine cost_gencost_seaicev4(myThid)
2ac072a19d Gael*0014 
                0015 c     ==================================================================
                0016 c     SUBROUTINE cost_gencost_seaicev4
                0017 c     ==================================================================
                0018 c
                0019 c     o Evaluate cost function contributions of ice concentration.
                0020 c
                0021 c     ==================================================================
                0022 c     SUBROUTINE cost_gencost_seaicev4
                0023 c     ==================================================================
                0024 
                0025       implicit none
                0026 
                0027 c     == global variables ==
                0028 
                0029 #include "EEPARAMS.h"
                0030 #include "SIZE.h"
                0031 #include "PARAMS.h"
                0032 #include "GRID.h"
                0033 #ifdef ALLOW_CAL
                0034 # include "cal.h"
                0035 #endif
f09238ab8f Gael*0036 #ifdef ALLOW_ECCO
13d362b8c1 Ou W*0037 # include "ECCO_SIZE.h"
                0038 # include "ECCO.h"
f09238ab8f Gael*0039 #endif
2b959ba38e Mart*0040 c#ifdef ALLOW_SEAICE
                0041 c# include "SEAICE_SIZE.h"
                0042 c# include "SEAICE_COST.h"
                0043 c# include "SEAICE_PARAMS.h"
                0044 c#endif
2ac072a19d Gael*0045 
                0046 c     == routine arguments ==
f8e779c983 antn*0047       integer myThid
2ac072a19d Gael*0048 
f09238ab8f Gael*0049 #ifdef ALLOW_SEAICE
2ac072a19d Gael*0050 #ifdef ALLOW_GENCOST_CONTRIBUTION
                0051 
                0052 c     == local variables ==
                0053 
                0054       integer nrecloc
                0055       integer localstartdate(4)
                0056 
0a8c2c2ff2 An T*0057 catn changing names to make more self-explanatory
                0058 c old:sst  -> model has deficiency in iceconc -> new:deconc
272838bfc2 An T*0059 c old:heff -> model has excess of iceconc     -> new:exconc
0a8c2c2ff2 An T*0060 
f8e779c983 antn*0061       _RL areabar    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
                0062       _RL deconcbar  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
                0063       _RL exconcbar  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
                0064       _RL localweight  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
7e836c24c5 An T*0065       _RL dummy
2ac072a19d Gael*0066       _RL localperiod
                0067       _RL spminloc
                0068       _RL spmaxloc
                0069       _RL spzeroloc
                0070 
272838bfc2 An T*0071       character*(MAX_LEN_FNAM) areabarfile
                0072       character*(MAX_LEN_FNAM) deconcbarfile
                0073       character*(MAX_LEN_FNAM) exconcbarfile
2ac072a19d Gael*0074       character*(MAX_LEN_FNAM) localobsfile
                0075 
0a8c2c2ff2 An T*0076       integer igen_conc, igen_deconc, igen_exconc
2ac072a19d Gael*0077 
                0078       integer bi,bj
11c3150c71 Mart*0079       integer k
7e836c24c5 An T*0080       integer irec, jrec
13e29837cc An T*0081       integer  il, k2
2ac072a19d Gael*0082       integer localrec
                0083       integer obsrec
13e29837cc An T*0084       logical dosumsq, dovarwei
2ac072a19d Gael*0085 
7e836c24c5 An T*0086       integer preproc_i(NGENPPROC)
                0087       _RL preproc_r(NGENPPROC)
                0088       character*(MAX_LEN_FNAM) preproc(NGENPPROC)
                0089       character*(MAX_LEN_FNAM) preproc_c(NGENPPROC)
2ac072a19d Gael*0090 
11c3150c71 Mart*0091       _RL localmask  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
                0092       _RL localobs   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
                0093       _RL localtmp   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
                0094       _RL localdif   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
                0095       _RL difmask    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
                0096       _RL difmask1   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
2ac072a19d Gael*0097 
7e836c24c5 An T*0098       character*(128) fname0, fname0w, fname1
2ac072a19d Gael*0099 
8227586f14 An T*0100       character*(MAX_LEN_FNAM) localobswfile
2ac072a19d Gael*0101       logical exst
                0102 
                0103 c     == external functions ==
                0104 
                0105       integer  ilnblnk
                0106       external ilnblnk
                0107 
                0108 c     == end of interface ==
                0109 
272838bfc2 An T*0110 c=============== PART 0: initilization ===================
13e29837cc An T*0111 
2ac072a19d Gael*0112 c-- detect the relevant gencost indices
                0113       igen_conc=0
0a8c2c2ff2 An T*0114       igen_deconc=0
                0115       igen_exconc=0
2ac072a19d Gael*0116       do k=1,NGENCOST
                0117         if (gencost_name(k).EQ.'siv4-conc') igen_conc=k
0a8c2c2ff2 An T*0118         if (gencost_name(k).EQ.'siv4-deconc') igen_deconc=k
                0119         if (gencost_name(k).EQ.'siv4-exconc') igen_exconc=k
2ac072a19d Gael*0120       enddo
                0121 
f8e779c983 antn*0122 c-- Dependency:
272838bfc2 An T*0123 c A) igen_conc can exist on its own
                0124 c B) igen_deconc needs igen_conc
                0125 c C) igen_exconc needs both igen_conc and igen_deconc
                0126       if (igen_conc.NE.0) then
2ac072a19d Gael*0127 
7e836c24c5 An T*0128 c-- initialize objf and num:
11c3150c71 Mart*0129       do bj = myByLo(myThid), myByHi(myThid)
                0130         do bi = myBxLo(myThid), myBxHi(myThid)
7e836c24c5 An T*0131           objf_gencost(bi,bj,igen_conc) = 0. _d 0
272838bfc2 An T*0132           num_gencost(bi,bj,igen_conc)  = 0. _d 0
                0133           if(igen_deconc.ne.0) then
                0134             objf_gencost(bi,bj,igen_deconc) = 0. _d 0
                0135             num_gencost(bi,bj,igen_deconc)  = 0. _d 0
                0136           endif
                0137           if(igen_exconc.ne.0) then
                0138             objf_gencost(bi,bj,igen_exconc) = 0. _d 0
                0139             num_gencost(bi,bj,igen_exconc)  = 0. _d 0
                0140           endif
7e836c24c5 An T*0141         enddo
                0142       enddo
                0143 
2ac072a19d Gael*0144 c--   Initialise local variables.
8227586f14 An T*0145       nrecloc=0
                0146       localperiod=0.
                0147 
272838bfc2 An T*0148       areabarfile=gencost_barfile(igen_conc)
                0149       if(igen_deconc.ne.0) deconcbarfile=gencost_barfile(igen_deconc)
                0150       if(igen_exconc.ne.0) exconcbarfile=gencost_barfile(igen_exconc)
                0151 
8fed38679a Gael*0152       localobsfile=gencost_datafile(igen_conc)
8227586f14 An T*0153       localobswfile=gencost_errfile(igen_conc)
7e836c24c5 An T*0154       dummy=gencost_dummy(igen_conc)
2ac072a19d Gael*0155       localstartdate(1)=modelstartdate(1)
                0156       localstartdate(2)=modelstartdate(2)
                0157       localstartdate(3)=modelstartdate(3)
                0158       localstartdate(4)=modelstartdate(4)
                0159       spminloc=gencost_spmin(igen_conc)
                0160       spmaxloc=gencost_spmax(igen_conc)
                0161       spzeroloc=gencost_spzero(igen_conc)
272838bfc2 An T*0162 
8227586f14 An T*0163       localperiod=gencost_period(igen_conc)
                0164       nrecloc=gencost_nrec(igen_conc)
2ac072a19d Gael*0165 
13e29837cc An T*0166 c-- flag to add cost: true=(obs-mod)*(obs-mod)*weight
                0167       dosumsq=.TRUE.
                0168       dovarwei=.FALSE.
                0169       do k2 = 1, NGENPPROC
                0170         preproc(k2)=gencost_preproc(k2,igen_conc)
                0171         preproc_i(k2)=gencost_preproc_i(k2,igen_conc)
                0172         preproc_c(k2)=gencost_preproc_c(k2,igen_conc)
                0173         preproc_r(k2)=gencost_preproc_r(k2,igen_conc)
                0174         if (preproc(k2).EQ.'variaweight') dovarwei=.TRUE.
                0175         if (preproc(k2).EQ.'nosumsq') dosumsq=.FALSE.
                0176       enddo
272838bfc2 An T*0177 
11c3150c71 Mart*0178 c--   initialize arrays, copy 2D maskInC to 2D localmask; this means
                0179 c     that sea ice cost function contributions on open boundary points
                0180 c     are masked out, probably something that eccov4 can live
                0181 c     with. Alternatively, on has to copy the first level of maskC to
                0182 c     localmask explicitly.
7e836c24c5 An T*0183       call ecco_zero(localobs,1,spzeroloc,myThid)
                0184       call ecco_zero(localweight,1,zeroRL,myThid)
11c3150c71 Mart*0185       call ecco_zero(localmask,1,zeroRL,myThid)
                0186       call ecco_cprsrl(maskInC,localmask,1,1,myThid)
2ac072a19d Gael*0187 
272838bfc2 An T*0188 c=============== PART 1: main loop ===================
2ac072a19d Gael*0189       if ( .NOT. ( localobsfile.EQ.' ' ) ) then
                0190 
                0191 c--   Loop over records for the second time.
                0192       do irec = 1, nrecloc
                0193 
272838bfc2 An T*0194 c====================================================
                0195 c--------- PART 1.1 read weights --------------------
                0196 c====================================================
7e836c24c5 An T*0197         exst=.FALSE.
                0198         jrec=1
13e29837cc An T*0199         if( dovarwei ) jrec = irec
272838bfc2 An T*0200         call cost_gencal(areabarfile,gencost_errfile(igen_conc),
7e836c24c5 An T*0201      &     jrec, localstartdate, localperiod, fname1,
f8e779c983 antn*0202      &     fname0w, localrec, obsrec, exst, myThid)
11c3150c71 Mart*0203         call ecco_zero(localweight,1,zeroRL,myThid)
6b2230d510 Ou W*0204 #ifdef SEAICECOST_JPL
                0205        fname0w=gencost_errfile(igen_conc)
b9b1c4d04c Timo*0206        call ecco_readwei(fname0w,localweight,localrec,
11c3150c71 Mart*0207      &      1,1,dosumsq,myThid)
6b2230d510 Ou W*0208        call ecco_readwei(gencost_errfile(igen_deconc),
ebcc29af97 Jean*0209      &      gencost_weight(1-OLx,1-OLy,1,1,igen_deconc),localrec,
11c3150c71 Mart*0210      &      1,1,dosumsq,myThid)
6b2230d510 Ou W*0211        call ecco_readwei(gencost_errfile(igen_exconc),
ebcc29af97 Jean*0212      &      gencost_weight(1-OLx,1-OLy,1,1,igen_exconc),localrec,
11c3150c71 Mart*0213      &      1,1,dosumsq,myThid)
6b2230d510 Ou W*0214 #else
7e836c24c5 An T*0215         if ( (localrec. GT. 0).AND.(obsrec .GT. 0).AND.(exst) ) then
b9b1c4d04c Timo*0216           call ecco_readwei(fname0w,localweight,localrec,
11c3150c71 Mart*0217      &      1,1,dosumsq,myThid)
2ac072a19d Gael*0218         else
7e836c24c5 An T*0219           WRITE(standardMessageUnit,'(A)')
                0220      &     'siv4cost WARNING: ALL WEIGHTS ZEROS! NO CONTRIBUTION'
2ac072a19d Gael*0221         endif
6b2230d510 Ou W*0222 #endif
2ac072a19d Gael*0223 
272838bfc2 An T*0224 c====================================================
                0225 c--------- PART 1.2 read barfiles ------------------
                0226 c====================================================
                0227 c-- set all bars to zeros:
11c3150c71 Mart*0228         call ecco_zero(areabar,1,zeroRL,myThid)
                0229         call ecco_zero(deconcbar,1,zeroRL,myThid)
                0230         call ecco_zero(exconcbar,1,zeroRL,myThid)
272838bfc2 An T*0231 
                0232 c--1.2.A sea-ice concentration barfile
7e836c24c5 An T*0233         exst=.FALSE.
272838bfc2 An T*0234         call cost_gencal(areabarfile,gencost_datafile(igen_conc),
                0235      &   irec,localstartdate,localperiod,fname1,
f8e779c983 antn*0236      &   fname0,localrec,obsrec,exst,myThid)
11c3150c71 Mart*0237         call cost_genread(fname1,areabar,localtmp,irec,1,1,
f8e779c983 antn*0238      &       nrecloc,preproc,preproc_c,preproc_i,preproc_r,
                0239      &       dummy,myThid)
272838bfc2 An T*0240 
                0241 c--1.2.B sst as proxy for deconc barfile, needs igen_conc
                0242         if(igen_deconc.ne.0) then
                0243          exst=.FALSE.
                0244          call cost_gencal(deconcbarfile,gencost_datafile(igen_conc),
                0245      &    irec,localstartdate,localperiod,fname1,
f8e779c983 antn*0246      &    fname0,localrec,obsrec,exst,myThid)
11c3150c71 Mart*0247          call cost_genread(fname1,deconcbar,localtmp,
                0248      &        irec,1,1,
                0249      &        nrecloc,preproc,preproc_c,preproc_i,preproc_r,
                0250      &        dummy,myThid)
272838bfc2 An T*0251         endif
                0252 
                0253 c--1.2.C heff as proxy for exconc barfile, need igen_conc and igen_exconc
                0254         if(igen_deconc.ne.0 .and. igen_exconc.ne.0) then
                0255          exst=.FALSE.
                0256          call cost_gencal(exconcbarfile,gencost_datafile(igen_conc),
                0257      &    irec,localstartdate,localperiod,fname1,
f8e779c983 antn*0258      &    fname0,localrec,obsrec,exst,myThid)
11c3150c71 Mart*0259          call cost_genread(fname1,exconcbar,localtmp,
                0260      &        irec,1,1,
                0261      &        nrecloc,preproc,preproc_c,preproc_i,preproc_r,
                0262      &        dummy,myThid)
272838bfc2 An T*0263         endif
                0264 
                0265 c====================================================
                0266 c--------- PART 1.3 read data --------------------
                0267 c====================================================
                0268 c-- initialize to spzerloc = -9999.
11c3150c71 Mart*0269         call ecco_zero(localobs,1,spzeroloc,myThid)
7e836c24c5 An T*0270         if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) ) then
11c3150c71 Mart*0271          CALL READ_REC_3D_RL( fname0, cost_iprec, 1,
f8e779c983 antn*0272      &                        localobs, localrec, 0, myThid )
2ac072a19d Gael*0273         else
7e836c24c5 An T*0274           il=ilnblnk( fname0 )
                0275           WRITE(standardMessageUnit,'(2A)')
                0276      &     'siv4cost WARNING: DATA MISING! NO CONTRIBUTION, ',
                0277      &     fname0(1:il)
2ac072a19d Gael*0278         endif
                0279 
272838bfc2 An T*0280 c====================================================
                0281 c--------- PART 1.4 Cost calculation -------------
                0282 c====================================================
f8e779c983 antn*0283 c compute obs minus bar (localdif) and mask (difmask)
11c3150c71 Mart*0284         call ecco_zero(localdif,1,zeroRL,myThid)
                0285         call ecco_zero(difmask,1,zeroRL,myThid)
f8e779c983 antn*0286         call ecco_diffmsk(
11c3150c71 Mart*0287      I     areabar, localobs, localmask,
                0288      I     1, 1, spminloc, spmaxloc, spzeroloc,
7e836c24c5 An T*0289      O     localdif, difmask,
                0290      I     myThid )
8227586f14 An T*0291 
272838bfc2 An T*0292 c---1.4.A area term:
05f114eba6 An T*0293         call ecco_addcost(
11c3150c71 Mart*0294      I     localdif,localweight,difmask,1,1,dosumsq,
05f114eba6 An T*0295      O     objf_gencost(1,1,igen_conc),num_gencost(1,1,igen_conc),
f8e779c983 antn*0296      I     myThid)
05f114eba6 An T*0297 
272838bfc2 An T*0298 c---1.4.B defficient ice term: (old: sst term, new: deconc)
                0299 c Add ice: model_A==0 but obs_A > 0, calc enthalpy E:
                0300         if(igen_deconc.ne.0) then
11c3150c71 Mart*0301          call ecco_zero(difmask1,1,zeroRL,myThid)
                0302          call ecco_zero(localdif,1,zeroRL,myThid)
                0303          call ecco_zero(localtmp,1,zeroRL,myThid)
f8e779c983 antn*0304 
272838bfc2 An T*0305          call get_exconc_deconc(
11c3150c71 Mart*0306      I    localobs,1,areabar,exconcbar,deconcbar,1,
272838bfc2 An T*0307      I    difmask,'de',
65756801a8 An T*0308      O    localdif,difmask1,localtmp,
                0309      I    myThid )
6b2230d510 Ou W*0310 #ifdef SEAICECOST_JPL
ebcc29af97 Jean*0311         call ecco_cp( gencost_weight(1-OLx,1-OLy,1,1,igen_deconc),
11c3150c71 Mart*0312      O       localtmp,1,1,myThid)
6b2230d510 Ou W*0313 #endif
05f114eba6 An T*0314         call ecco_addcost(
11c3150c71 Mart*0315      I      localdif,localtmp,difmask1,1,1,dosumsq,
272838bfc2 An T*0316      O      objf_gencost(1,1,igen_deconc),num_gencost(1,1,igen_deconc),
f8e779c983 antn*0317      I      myThid)
272838bfc2 An T*0318         endif
1bedc5345f An T*0319 
272838bfc2 An T*0320 c---1.4.C excessive ice term:  (old: heff and sst term, new: exconc)
                0321 c Removing ice: model_A > 0 but obs_A==0, calc enthalpy E:
                0322         if(igen_deconc.ne.0 .and. igen_exconc.ne.0) then
11c3150c71 Mart*0323          call ecco_zero(difmask1,1,zeroRL,myThid)
                0324          call ecco_zero(localdif,1,zeroRL,myThid)
                0325          call ecco_zero(localtmp,1,zeroRL,myThid)
272838bfc2 An T*0326 
                0327          call get_exconc_deconc(
11c3150c71 Mart*0328      I    localobs,1,areabar,exconcbar,deconcbar,1,
272838bfc2 An T*0329      I    difmask,'ex',
65756801a8 An T*0330      O    localdif,difmask1,localtmp,
                0331      I    myThid )
6b2230d510 Ou W*0332 #ifdef SEAICECOST_JPL
ebcc29af97 Jean*0333         call ecco_cp( gencost_weight(1-OLx,1-OLy,1,1,igen_exconc),
11c3150c71 Mart*0334      O       localtmp,1,1,myThid)
6b2230d510 Ou W*0335 #endif
05f114eba6 An T*0336         call ecco_addcost(
11c3150c71 Mart*0337      I      localdif,localtmp,difmask1,1,1,dosumsq,
272838bfc2 An T*0338      O      objf_gencost(1,1,igen_exconc),num_gencost(1,1,igen_exconc),
f8e779c983 antn*0339      I      myThid)
272838bfc2 An T*0340         endif
2ac072a19d Gael*0341 
                0342       enddo
                0343 
                0344       endif !if ( .NOT. ( localobsfile.EQ.' ' ) ) then
272838bfc2 An T*0345       endif !if (igen_conc.NE.0)
2ac072a19d Gael*0346 
                0347 #endif /* ALLOW_GENCOST_CONTRIBUTION */
f09238ab8f Gael*0348 #endif /* ALLOW_SEAICE */
2ac072a19d Gael*0349 
65756801a8 An T*0350       RETURN
2b959ba38e Mart*0351       END
65756801a8 An T*0352 
                0353 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0354 
                0355       subroutine get_exconc_deconc(
                0356      I    localobs,nnzobs,concbar,exconcbar,deconcbar,nnzbar,
                0357      I    localmask,flag_exconc_deconc,
                0358      O    localfld,localfldmsk,localfldweight,
                0359      I    myThid )
                0360 
                0361 C     !DESCRIPTION: \bv
f8e779c983 antn*0362 c     Routine to calculate Enthalpy for the case of
65756801a8 An T*0363 c     defficient/excessive model seaice
                0364 C     \ev
                0365 
                0366 C     !USES:
                0367       implicit none
                0368 
                0369 c     == global variables ==
                0370 #include "EEPARAMS.h"
                0371 #include "SIZE.h"
                0372 #include "PARAMS.h"
                0373 #include "GRID.h"
                0374 #ifdef ALLOW_SEAICE
                0375 # include "SEAICE_SIZE.h"
                0376 # include "SEAICE_COST.h"
                0377 # include "SEAICE_PARAMS.h"
                0378 #endif
                0379 #ifdef ALLOW_ECCO
13d362b8c1 Ou W*0380 # include "ECCO_SIZE.h"
                0381 # include "ECCO.h"
65756801a8 An T*0382 #endif
                0383 
                0384 c     == routine arguments ==
                0385 
f8e779c983 antn*0386       integer myThid, nnzbar, nnzobs
                0387       _RL localmask     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
65756801a8 An T*0388 
f8e779c983 antn*0389       _RL localobs      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
                0390       _RL concbar       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
                0391       _RL deconcbar     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
                0392       _RL exconcbar     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
                0393       _RL localfld      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
                0394       _RL localfldmsk   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
                0395       _RL localfldweight(1-OLx:sNx+OLx,1-OLy:sNy+OLy,1,nSx,nSy)
65756801a8 An T*0396 
                0397       character*2 flag_exconc_deconc
                0398 
                0399 #ifdef ALLOW_SEAICE
2b959ba38e Mart*0400 #ifdef ALLOW_GENCOST_CONTRIBUTION
65756801a8 An T*0401 
                0402 c    == local variables ==
                0403       integer bi,bj
                0404       integer i,j,k
                0405 
                0406 C- jmc: params SEAICE_freeze has been retired; set it locally until someone
                0407 C       who knows what this cost-cointribution means fix it.
                0408 C- atn: also adding 1 normalizing factor same order of magnitude as
f8e779c983 antn*0409 C       rhoNil*HeatCapacity_cp*dz = SEAICE_rhoice*SEAICE_lhFusion*heff
65756801a8 An T*0410 C       = 1e3*1e3*1e1=1e7
                0411 C- atn: lastly, define 2 cutoff values for cost to be read in from data.seaice
                0412 C      and initialized in seaice_readparms: SEAICE_cutoff_[area,heff]
                0413 C      Reason: some iceconc data set have "bogus" mask with area>0
                0414 C      at climatological max locations -> not real data.  So either need
                0415 C      to clean up the data or take SEAICE_cutoff_area>=0.15 for example.
                0416 C      Might need to migrate into pkg/ecco instead of in pkg/seaice.
                0417       _RL SEAICE_freeze, epsilonTemp, epsilonHEFF
                0418       _RL localnorm, localnormsq
                0419       _RL const1,const2
                0420 CEOP
                0421 
                0422       SEAICE_freeze  = -1.96  _d 0
                0423       epsilonTemp = 0.0001 _d 0
6b2230d510 Ou W*0424 #ifdef SEAICECOST_JPL
                0425       epsilonHEFF = 0.3 _d 0
                0426 #else
65756801a8 An T*0427       epsilonHEFF = 0.01 _d 0
6b2230d510 Ou W*0428 #endif
65756801a8 An T*0429       localnorm = 1. _d -07
                0430       localnormsq=localnorm*localnorm
                0431 
                0432       const1=HeatCapacity_Cp*rhoNil*drF(1)
                0433       const2=SEAICE_lhFusion*SEAICE_rhoIce
                0434 
                0435 c intialize
                0436       call ecco_zero(localfld,nnzobs,zeroRL,myThid)
                0437       call ecco_zero(localfldmsk,nnzobs,zeroRL,myThid)
                0438       call ecco_zero(localfldweight,nnzobs,zeroRL,myThid)
                0439 
272838bfc2 An T*0440 c----------------------DECONC-------------------------------
                0441 catn-- old: sst term, new: deconc
                0442 c needs localconcbar and localsstbar
                0443 c Add ice: model_A==0 but obs_A > 0, calc enthalpy E:
                0444 c E_current = (deconcbar(i,j,k,bi,bj)-Tfreeze)
                0445 c             *HeatCapacity_Cp*rhoNil*drF(1)
                0446 c HEFF_target = epsilon_HEFF [m]
                0447 c E_target  = -(HEFF_target*SEAICE_lhFusion*SEAICE_rhoIce)
f8e779c983 antn*0448 c cost=(Model-data)^2
272838bfc2 An T*0449       if(flag_exconc_deconc.EQ.'de') then
11c3150c71 Mart*0450         do bj = myByLo(myThid), myByHi(myThid)
                0451           do bi = myBxLo(myThid), myBxHi(myThid)
65756801a8 An T*0452            do k = 1,nnzobs
11c3150c71 Mart*0453             do j = 1-OLy,sNy+OLy
                0454              do i = 1-OLx,sNx+OLx
65756801a8 An T*0455 
272838bfc2 An T*0456               if ( (concbar(i,j,k,bi,bj) .LE. 0.).AND.
                0457      &               (localobs(i,j,k,bi,bj) .GT. 0.) ) then
65756801a8 An T*0458 
                0459                localfldmsk(i,j,k,bi,bj) = localmask(i,j,k,bi,bj)
272838bfc2 An T*0460 
6b2230d510 Ou W*0461 #ifndef SEAICECOST_JPL
                0462                localfldweight(i,j,k,bi,bj) = localnormsq
                0463 #endif
272838bfc2 An T*0464                localfld(i,j,k,bi,bj) =
                0465      &          (deconcbar(i,j,k,bi,bj)-SEAICE_freeze)*const1
                0466      &               - (-1. _d 0 *epsilonHEFF*const2)
65756801a8 An T*0467               endif
                0468              enddo
                0469             enddo
                0470            enddo
                0471           enddo
                0472         enddo
                0473       endif
                0474 
272838bfc2 An T*0475 c----------------------EXCONC-------------------------------
                0476 catn-- old: heff and sst term, new: exconc
                0477 c needs localconcbar, localsstbar, and localheffbar
                0478 c Removing ice: model_A > 0 but obs_A==0, calc enthalpy E:
                0479 c E_current = [(deconcbar-SEAICE_freeze)*HeatCapacity_Cp*rhoNil*drF(1)
f8e779c983 antn*0480 c            - (exconcbar * SEAICE_lhFusion * SEAICE_rhoIce)
272838bfc2 An T*0481 c            - (HSNOW * SEAICE_lhFusion * SEAICE_rhoSnow)]
                0482 c E_target = (epsilonTemp) * HeatCapacity_Cp * rhoNil * drF(1)
                0483 c cost(Model-data)^2
                0484 
                0485       if(flag_exconc_deconc.EQ.'ex') then
11c3150c71 Mart*0486         do bj = myByLo(myThid), myByHi(myThid)
                0487           do bi = myBxLo(myThid), myBxHi(myThid)
65756801a8 An T*0488            do k = 1,nnzobs
11c3150c71 Mart*0489             do j = 1-OLy,sNy+OLy
                0490              do i = 1-OLx,sNx+OLx
65756801a8 An T*0491 
272838bfc2 An T*0492               if ((localobs(i,j,k,bi,bj) .LE. SEAICE_cutoff_area).AND.
                0493      &        (exconcbar(i,j,k,bi,bj) .GT. SEAICE_cutoff_heff)) then
65756801a8 An T*0494 
                0495                localfldmsk(i,j,k,bi,bj) = localmask(i,j,k,bi,bj)
6b2230d510 Ou W*0496 #ifndef SEAICECOST_JPL
65756801a8 An T*0497                localfldweight(i,j,k,bi,bj) = localnormsq
6b2230d510 Ou W*0498 #endif
f8e779c983 antn*0499                localfld(i,j,k,bi,bj) =
272838bfc2 An T*0500      &        ( (deconcbar(i,j,k,bi,bj)-SEAICE_freeze)*const1
6b2230d510 Ou W*0501 #ifdef SEAICECOST_JPL
                0502      &         - max(exconcbar(i,j,k,bi,bj),epsilonHEFF)*
f8e779c983 antn*0503 #else
6b2230d510 Ou W*0504      &         - exconcbar(i,j,k,bi,bj)*
                0505 #endif
                0506      &           const2 ) - (epsilonTemp*const1)
65756801a8 An T*0507               endif
                0508              enddo
                0509             enddo
                0510            enddo
                0511           enddo
                0512         enddo
                0513       endif
                0514 
                0515 #endif /* ALLOW_GENCOST_CONTRIBUTION */
                0516 #endif /* ALLOW_SEAICE */
                0517       RETURN
                0518       END