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-|--+----|