C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_gencost_seaicev4.F,v 1.20 2017/04/03 23:16:38 ou.wang Exp $
C $Name:  $

#include "ECCO_OPTIONS.h"

      subroutine COST_GENCOST_SEAICEV4(mythid)

c     ==================================================================
c     SUBROUTINE cost_gencost_seaicev4
c     ==================================================================
c
c     o Evaluate cost function contributions of ice concentration.
c
c     ==================================================================
c     SUBROUTINE cost_gencost_seaicev4
c     ==================================================================

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"
#ifdef ALLOW_CAL
# include "cal.h"
#endif
#ifdef ALLOW_ECCO
# include "ecco.h"
#endif
#ifdef ALLOW_SEAICE
# include "SEAICE_SIZE.h"
# include "SEAICE_COST.h"
# include "SEAICE_PARAMS.h"
#endif

c     == routine arguments ==
      integer mythid

#ifdef ALLOW_SEAICE
#ifdef ALLOW_GENCOST_CONTRIBUTION

c     == local variables ==

      integer nnzsiv4, nnzbar
      parameter (nnzsiv4 = 1 , nnzbar = 1)
      integer nrecloc
      integer localstartdate(4)

catn changing names to make more self-explanatory
c old:sst  -> model has deficiency in iceconc -> new:deconc
c old:heff -> model has excess of iceconc     -> new:exconc

      _RL areabar    (1-olx:snx+olx,1-oly:sny+oly,1,nsx,nsy)
      _RL deconcbar  (1-olx:snx+olx,1-oly:sny+oly,1,nsx,nsy)
      _RL exconcbar  (1-olx:snx+olx,1-oly:sny+oly,1,nsx,nsy)
      _RL localweight  (1-olx:snx+olx,1-oly:sny+oly,1,nsx,nsy)
      _RL dummy
      _RL localperiod
      _RL spminloc
      _RL spmaxloc
      _RL spzeroloc

      character*(MAX_LEN_FNAM) areabarfile
      character*(MAX_LEN_FNAM) deconcbarfile
      character*(MAX_LEN_FNAM) exconcbarfile
      character*(MAX_LEN_FNAM) localobsfile

      integer igen_conc, igen_deconc, igen_exconc

      integer bi,bj
      integer i,j,k
      integer itlo,ithi
      integer jtlo,jthi
      integer jmin,jmax
      integer imin,imax
      integer irec, jrec
      integer  il, k2
      integer localrec
      integer obsrec
      logical dosumsq, dovarwei

      integer preproc_i(NGENPPROC)
      _RL preproc_r(NGENPPROC)
      character*(MAX_LEN_FNAM) preproc(NGENPPROC)
      character*(MAX_LEN_FNAM) preproc_c(NGENPPROC)

      _RL localmask  (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)

      _RL localobs   (1-olx:snx+olx,1-oly:sny+oly,nnzsiv4,nsx,nsy)
      _RL localtmp   (1-olx:snx+olx,1-oly:sny+oly,nnzsiv4,nsx,nsy)
      _RL localdif   (1-olx:snx+olx,1-oly:sny+oly,nnzsiv4,nsx,nsy)
      _RL difmask    (1-olx:snx+olx,1-oly:sny+oly,nnzsiv4,nsx,nsy)
      _RL difmask1   (1-olx:snx+olx,1-oly:sny+oly,nnzsiv4,nsx,nsy)

      character*(128) fname0, fname0w, fname1

      character*(MAX_LEN_FNAM) localobswfile
      logical exst

c     == external functions ==

      integer  ilnblnk
      external 

c     == end of interface ==

      jtlo = mybylo(mythid)
      jthi = mybyhi(mythid)
      itlo = mybxlo(mythid)
      ithi = mybxhi(mythid)
      jmin = 1
      jmax = sny
      imin = 1
      imax = snx

c=============== PART 0: initilization ===================

c-- detect the relevant gencost indices
      igen_conc=0
      igen_deconc=0
      igen_exconc=0
      do k=1,NGENCOST
        if (gencost_name(k).EQ.'siv4-conc') igen_conc=k
        if (gencost_name(k).EQ.'siv4-deconc') igen_deconc=k
        if (gencost_name(k).EQ.'siv4-exconc') igen_exconc=k
      enddo

c-- Dependency: 
c A) igen_conc can exist on its own
c B) igen_deconc needs igen_conc
c C) igen_exconc needs both igen_conc and igen_deconc
      if (igen_conc.NE.0) then

c-- initialize objf and num:
      do bj = jtlo,jthi
        do bi = itlo,ithi
          objf_gencost(bi,bj,igen_conc) = 0. _d 0
          num_gencost(bi,bj,igen_conc)  = 0. _d 0
          if(igen_deconc.ne.0) then
            objf_gencost(bi,bj,igen_deconc) = 0. _d 0
            num_gencost(bi,bj,igen_deconc)  = 0. _d 0
          endif
          if(igen_exconc.ne.0) then
            objf_gencost(bi,bj,igen_exconc) = 0. _d 0
            num_gencost(bi,bj,igen_exconc)  = 0. _d 0
          endif
        enddo
      enddo

c--   Initialise local variables.
      nrecloc=0
      localperiod=0.

      areabarfile=gencost_barfile(igen_conc)
      if(igen_deconc.ne.0) deconcbarfile=gencost_barfile(igen_deconc)
      if(igen_exconc.ne.0) exconcbarfile=gencost_barfile(igen_exconc)

      localobsfile=gencost_datafile(igen_conc)
      localobswfile=gencost_errfile(igen_conc)
      dummy=gencost_dummy(igen_conc)
      localstartdate(1)=modelstartdate(1)
      localstartdate(2)=modelstartdate(2)
      localstartdate(3)=modelstartdate(3)
      localstartdate(4)=modelstartdate(4)
      spminloc=gencost_spmin(igen_conc)
      spmaxloc=gencost_spmax(igen_conc)
      spzeroloc=gencost_spzero(igen_conc)

      localperiod=gencost_period(igen_conc)
      nrecloc=gencost_nrec(igen_conc)

c-- flag to add cost: true=(obs-mod)*(obs-mod)*weight
      dosumsq=.TRUE.
      dovarwei=.FALSE.
      do k2 = 1, NGENPPROC
        preproc(k2)=gencost_preproc(k2,igen_conc)
        preproc_i(k2)=gencost_preproc_i(k2,igen_conc)
        preproc_c(k2)=gencost_preproc_c(k2,igen_conc)
        preproc_r(k2)=gencost_preproc_r(k2,igen_conc)
        if (preproc(k2).EQ.'variaweight') dovarwei=.TRUE.
        if (preproc(k2).EQ.'nosumsq') dosumsq=.FALSE.
      enddo

c-- initialize arrays, copy maskC to localmask
      call ECCO_ZERO(localobs,1,spzeroloc,myThid)
      call ECCO_ZERO(localweight,1,zeroRL,myThid)
      call ECCO_ZERO(localmask,Nr,zeroRL,myThid)
      call ECCO_CPRSRL(maskC,Nr,localmask,Nr,myThid)

c=============== PART 1: main loop ===================
      if ( .NOT. ( localobsfile.EQ.' ' ) ) then

c--   Loop over records for the second time.
      do irec = 1, nrecloc

c====================================================
c--------- PART 1.1 read weights --------------------
c====================================================
        exst=.FALSE.
        jrec=1
        if( dovarwei ) jrec = irec
        call COST_GENCAL(areabarfile,gencost_errfile(igen_conc),
     &     jrec, localstartdate, localperiod, fname1,
     &     fname0w, localrec, obsrec, exst, mythid)
        call ECCO_ZERO(localweight,nnzsiv4,zeroRL,myThid)
#ifdef SEAICECOST_JPL
       fname0w=gencost_errfile(igen_conc)
       call ECCO_READWEI(fname0w,localweight,localrec,nnzsiv4,mythid)
       call ECCO_READWEI(gencost_errfile(igen_deconc),
     &      gencost_weight(1,1,1,1,igen_deconc),localrec,nnzsiv4,mythid)
       call ECCO_READWEI(gencost_errfile(igen_exconc),
     &      gencost_weight(1,1,1,1,igen_exconc),localrec,nnzsiv4,mythid)
#else
        if ( (localrec. GT. 0).AND.(obsrec .GT. 0).AND.(exst) ) then
          call ECCO_READWEI(fname0w,localweight,localrec,nnzsiv4,mythid)
        else
          WRITE(standardMessageUnit,'(A)')
     &     'siv4cost WARNING: ALL WEIGHTS ZEROS! NO CONTRIBUTION'
        endif
#endif

c====================================================
c--------- PART 1.2 read barfiles ------------------
c====================================================
c-- set all bars to zeros:
        call ECCO_ZERO(areabar,nnzbar,zeroRL,myThid)
        call ECCO_ZERO(deconcbar,nnzbar,zeroRL,myThid)
        call ECCO_ZERO(exconcbar,nnzbar,zeroRL,myThid)

c--1.2.A sea-ice concentration barfile
        exst=.FALSE.
        call COST_GENCAL(areabarfile,gencost_datafile(igen_conc),
     &   irec,localstartdate,localperiod,fname1,
     &   fname0,localrec,obsrec,exst,mythid)
        call COST_GENREAD(fname1,areabar,localtmp,irec,nnzbar,
     &       nrecloc,preproc,preproc_c,preproc_i,preproc_r, 
     &       dummy,mythid)

c--1.2.B sst as proxy for deconc barfile, needs igen_conc
        if(igen_deconc.ne.0) then
         exst=.FALSE.
         call COST_GENCAL(deconcbarfile,gencost_datafile(igen_conc),
     &    irec,localstartdate,localperiod,fname1,
     &    fname0,localrec,obsrec,exst,mythid)
         call COST_GENREAD(fname1,deconcbar,localtmp,irec,nnzbar,
     &       nrecloc,preproc,preproc_c,preproc_i,preproc_r, 
     &       dummy,mythid)
        endif

c--1.2.C heff as proxy for exconc barfile, need igen_conc and igen_exconc
        if(igen_deconc.ne.0 .and. igen_exconc.ne.0) then
         exst=.FALSE.
         call COST_GENCAL(exconcbarfile,gencost_datafile(igen_conc),
     &    irec,localstartdate,localperiod,fname1,
     &    fname0,localrec,obsrec,exst,mythid)
         call COST_GENREAD(fname1,exconcbar,localtmp,irec,nnzbar,
     &       nrecloc,preproc,preproc_c,preproc_i,preproc_r, 
     &       dummy,mythid)
        endif

c====================================================
c--------- PART 1.3 read data --------------------
c====================================================
c-- initialize to spzerloc = -9999.
        call ECCO_ZERO(localobs,nnzsiv4,spzeroloc,mythid)
        if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) ) then
          call MDSREADFIELD( fname0, cost_iprec, cost_yftype, nnzsiv4,
     &         localobs, localrec, mythid )
        else
          il=ilnblnk( fname0 )
          WRITE(standardMessageUnit,'(2A)')
     &     'siv4cost WARNING: DATA MISING! NO CONTRIBUTION, ',
     &     fname0(1:il)
        endif

c====================================================
c--------- PART 1.4 Cost calculation -------------
c====================================================
c compute obs minus bar (localdif) and mask (difmask) 
        call ECCO_ZERO(localdif,nnzsiv4,zeroRL,mythid)
        call ECCO_ZERO(difmask,nnzsiv4,zeroRL,mythid)
        call ECCO_DIFFMSK( 
     I     areabar, nnzbar, localobs, nnzsiv4, localmask,
     I     spminloc, spmaxloc, spzeroloc,
     O     localdif, difmask,
     I     myThid )

c---1.4.A area term:
        call ECCO_ADDCOST(
     I     localdif,localweight,difmask,nnzsiv4,dosumsq,
     O     objf_gencost(1,1,igen_conc),num_gencost(1,1,igen_conc),
     I     mythid)

c---1.4.B defficient ice term: (old: sst term, new: deconc)
c Add ice: model_A==0 but obs_A > 0, calc enthalpy E:
        if(igen_deconc.ne.0) then
         call ECCO_ZERO(difmask1,nnzsiv4,zeroRL,mythid)
         call ECCO_ZERO(localdif,nnzsiv4,zeroRL,mythid)
         call ECCO_ZERO(localtmp,nnzsiv4,zeroRL,mythid)
 
         call GET_EXCONC_DECONC(
     I    localobs,nnzsiv4,areabar,exconcbar,deconcbar,nnzbar,
     I    difmask,'de',
     O    localdif,difmask1,localtmp,
     I    myThid )
#ifdef SEAICECOST_JPL
        call ECCO_CP(gencost_weight(1,1,1,1,igen_deconc),nnzsiv4,
     O       localtmp,nnzsiv4,myThid)
#endif
        call ECCO_ADDCOST(
     I      localdif,localtmp,difmask1,nnzsiv4,dosumsq,
     O      objf_gencost(1,1,igen_deconc),num_gencost(1,1,igen_deconc),
     I      mythid)
        endif

c---1.4.C excessive ice term:  (old: heff and sst term, new: exconc)
c Removing ice: model_A > 0 but obs_A==0, calc enthalpy E:
        if(igen_deconc.ne.0 .and. igen_exconc.ne.0) then
         call ECCO_ZERO(difmask1,nnzsiv4,zeroRL,mythid)
         call ECCO_ZERO(localdif,nnzsiv4,zeroRL,mythid)
         call ECCO_ZERO(localtmp,nnzsiv4,zeroRL,mythid)

         call GET_EXCONC_DECONC(
     I    localobs,nnzsiv4,areabar,exconcbar,deconcbar,nnzbar,
     I    difmask,'ex',
     O    localdif,difmask1,localtmp,
     I    myThid )
#ifdef SEAICECOST_JPL
        call ECCO_CP(gencost_weight(1,1,1,1,igen_exconc),nnzsiv4,
     O       localtmp,nnzsiv4,myThid)
#endif
        call ECCO_ADDCOST(
     I      localdif,localtmp,difmask1,nnzsiv4,dosumsq,
     O      objf_gencost(1,1,igen_exconc),num_gencost(1,1,igen_exconc),
     I      mythid)
        endif

      enddo

      endif !if ( .NOT. ( localobsfile.EQ.' ' ) ) then
      endif !if (igen_conc.NE.0)

#endif /* ALLOW_GENCOST_CONTRIBUTION */
#endif /* ALLOW_SEAICE */

      RETURN
      end


C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| #include "ECCO_OPTIONS.h" subroutine GET_EXCONC_DECONC( I localobs,nnzobs,concbar,exconcbar,deconcbar,nnzbar, I localmask,flag_exconc_deconc, O localfld,localfldmsk,localfldweight, I myThid ) C !DESCRIPTION: \bv c Routine to calculate Enthalpy for the case of c defficient/excessive model seaice C \ev C !USES: implicit none c == global variables == #include "EEPARAMS.h" #include "SIZE.h" #include "PARAMS.h" #include "GRID.h" #ifdef ALLOW_SEAICE # include "SEAICE_SIZE.h" # include "SEAICE_COST.h" # include "SEAICE_PARAMS.h" #endif #ifdef ALLOW_ECCO # include "ecco.h" #endif c == routine arguments == integer mythid, nnzbar, nnzobs _RL localmask (1-olx:snx+olx,1-oly:sny+oly,1,nsx,nsy) _RL localobs (1-olx:snx+olx,1-oly:sny+oly,1,nsx,nsy) _RL concbar (1-olx:snx+olx,1-oly:sny+oly,1,nsx,nsy) _RL deconcbar (1-olx:snx+olx,1-oly:sny+oly,1,nsx,nsy) _RL exconcbar (1-olx:snx+olx,1-oly:sny+oly,1,nsx,nsy) _RL localfld (1-olx:snx+olx,1-oly:sny+oly,1,nsx,nsy) _RL localfldmsk (1-olx:snx+olx,1-oly:sny+oly,1,nsx,nsy) _RL localfldweight(1-olx:snx+olx,1-oly:sny+oly,1,nsx,nsy) character*2 flag_exconc_deconc #ifdef ALLOW_GENCOST_CONTRIBUTION #ifdef ALLOW_SEAICE c == local variables == integer bi,bj integer itlo,ithi integer jtlo,jthi integer jmin,jmax integer imin,imax integer i,j,k C- jmc: params SEAICE_freeze has been retired; set it locally until someone C who knows what this cost-cointribution means fix it. C- atn: also adding 1 normalizing factor same order of magnitude as C rhoNil*HeatCapacity_cp*dz = SEAICE_rhoice*SEAICE_lhFusion*heff C = 1e3*1e3*1e1=1e7 C- atn: lastly, define 2 cutoff values for cost to be read in from data.seaice C and initialized in seaice_readparms: SEAICE_cutoff_[area,heff] C Reason: some iceconc data set have "bogus" mask with area>0 C at climatological max locations -> not real data. So either need C to clean up the data or take SEAICE_cutoff_area>=0.15 for example. C Might need to migrate into pkg/ecco instead of in pkg/seaice. _RL SEAICE_freeze, epsilonTemp, epsilonHEFF _RL localnorm, localnormsq _RL const1,const2 CEOP SEAICE_freeze = -1.96 _d 0 epsilonTemp = 0.0001 _d 0 #ifdef SEAICECOST_JPL epsilonHEFF = 0.3 _d 0 #else epsilonHEFF = 0.01 _d 0 #endif localnorm = 1. _d -07 localnormsq=localnorm*localnorm const1=HeatCapacity_Cp*rhoNil*drF(1) const2=SEAICE_lhFusion*SEAICE_rhoIce jtlo = mybylo(mythid) jthi = mybyhi(mythid) itlo = mybxlo(mythid) ithi = mybxhi(mythid) jmin = 1-oly jmax = sny+oly imin = 1-olx imax = snx+olx c intialize call ECCO_ZERO(localfld,nnzobs,zeroRL,myThid) call ECCO_ZERO(localfldmsk,nnzobs,zeroRL,myThid) call ECCO_ZERO(localfldweight,nnzobs,zeroRL,myThid) c----------------------DECONC------------------------------- catn-- old: sst term, new: deconc c needs localconcbar and localsstbar c Add ice: model_A==0 but obs_A > 0, calc enthalpy E: c E_current = (deconcbar(i,j,k,bi,bj)-Tfreeze) c *HeatCapacity_Cp*rhoNil*drF(1) c HEFF_target = epsilon_HEFF [m] c E_target = -(HEFF_target*SEAICE_lhFusion*SEAICE_rhoIce) c cost=(Model-data)^2 if(flag_exconc_deconc.EQ.'de') then do bj = jtlo,jthi do bi = itlo,ithi do k = 1,nnzobs do j = jmin,jmax do i = imin,imax if ( (concbar(i,j,k,bi,bj) .LE. 0.).AND. & (localobs(i,j,k,bi,bj) .GT. 0.) ) then localfldmsk(i,j,k,bi,bj) = localmask(i,j,k,bi,bj) #ifndef SEAICECOST_JPL localfldweight(i,j,k,bi,bj) = localnormsq #endif localfld(i,j,k,bi,bj) = & (deconcbar(i,j,k,bi,bj)-SEAICE_freeze)*const1 & - (-1. _d 0 *epsilonHEFF*const2) endif enddo enddo enddo enddo enddo endif c----------------------EXCONC------------------------------- catn-- old: heff and sst term, new: exconc c needs localconcbar, localsstbar, and localheffbar c Removing ice: model_A > 0 but obs_A==0, calc enthalpy E: c E_current = [(deconcbar-SEAICE_freeze)*HeatCapacity_Cp*rhoNil*drF(1) c - (exconcbar * SEAICE_lhFusion * SEAICE_rhoIce) c - (HSNOW * SEAICE_lhFusion * SEAICE_rhoSnow)] c E_target = (epsilonTemp) * HeatCapacity_Cp * rhoNil * drF(1) c cost(Model-data)^2 if(flag_exconc_deconc.EQ.'ex') then do bj = jtlo,jthi do bi = itlo,ithi do k = 1,nnzobs do j = jmin,jmax do i = imin,imax if ((localobs(i,j,k,bi,bj) .LE. SEAICE_cutoff_area).AND. & (exconcbar(i,j,k,bi,bj) .GT. SEAICE_cutoff_heff)) then localfldmsk(i,j,k,bi,bj) = localmask(i,j,k,bi,bj) #ifndef SEAICECOST_JPL localfldweight(i,j,k,bi,bj) = localnormsq #endif localfld(i,j,k,bi,bj) = & ( (deconcbar(i,j,k,bi,bj)-SEAICE_freeze)*const1 #ifdef SEAICECOST_JPL & - max(exconcbar(i,j,k,bi,bj),epsilonHEFF)* #else & - exconcbar(i,j,k,bi,bj)* #endif & const2 ) - (epsilonTemp*const1) endif enddo enddo enddo enddo enddo endif #endif /* ALLOW_GENCOST_CONTRIBUTION */ #endif /* ALLOW_SEAICE */ RETURN END


C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|