C $Header: /u/gcmpack/MITgcm/pkg/seaice/cost_ice_test.F,v 1.14 2014/12/03 16:47:45 mlosch Exp $
C $Name: $
#include "SEAICE_OPTIONS.h"
subroutine COST_ICE_TEST( mytime, myiter, mythid )
c ==================================================================
c SUBROUTINE cost_ice_test
c ==================================================================
c
c o Compute sea-ice cost function. The following options can be
c selected with data.seaice (SEAICE_PARM02) variable cost_ice_flag
c
c cost_ice_flag = 1
c - compute mean sea-ice volume
c costIceStart < mytime < costIceEnd
c
c cost_ice_flag = 2
c - compute mean sea-ice area
c costIceStart < mytime < costIceEnd
c
c cost_ice_flag = 3
c - heat content of top level plus latent heat of sea-ice
c costIceStart < mytime < costIceEnd
c
c cost_ice_flag = 4
c - heat content of top level
c costIceStart < mytime < costIceEnd
c
c cost_ice_flag = 5
c - heat content of top level plus sea-ice plus latent heat of snow
c costIceStart < mytime < costIceEnd
c
c cost_ice_flag = 6
c - quadratic cost function measuring difference between pkg/seaice
c AREA variable and simulated sea-ice measurements at every time
c step.
c
c ==================================================================
c
c started: menemenlis@jpl.nasa.gov 26-Feb-2003
c
c ==================================================================
c SUBROUTINE cost_ice_test
c ==================================================================
implicit none
c == global variables ==
#ifdef ALLOW_COST_ICE
#include "EEPARAMS.h"
#include "SIZE.h"
#include "GRID.h"
#include "PARAMS.h"
#include "SEAICE_SIZE.h"
#include "SEAICE_COST.h"
#include "SEAICE.h"
#include "DYNVARS.h"
#include "cost.h"
#endif /* ALLOW_COST_ICE */
c == routine arguments ==
_RL mytime
integer myiter
integer mythid
#ifdef ALLOW_COST_ICE
c == local variables ==
c msgBuf - Informational/error message buffer
CHARACTER*(MAX_LEN_MBUF) msgBuf
integer bi,bj,i,j,kSrf
_RL tempVar
c == external functions ==
integer ilnblnk
external
c == end of interface ==
if ( myTime .GT. (endTime - lastinterval) ) then
tempVar = 1. _d 0/
& ( ( 1. _d 0 + min(endTime-startTime,lastinterval) )
& / deltaTClock )
kSrf = 1
cph(
write(standardMessageUnit,*) 'ph-ice B ', myiter,
& theta(4,4,kSrf,1,1), area(4,4,1,1), heff(4,4,1,1)
cph)
if ( cost_ice_flag .eq. 1 ) then
c sea-ice volume
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do j = 1,sny
do i = 1,snx
objf_ice(bi,bj) = objf_ice(bi,bj) +
& tempVar * rA(i,j,bi,bj) * HEFF(i,j,bi,bj)
enddo
enddo
enddo
enddo
elseif ( cost_ice_flag .eq. 2 ) then
c sea-ice area
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do j = 1,sny
do i = 1,snx
objf_ice(bi,bj) = objf_ice(bi,bj) +
& tempVar * rA(i,j,bi,bj) * AREA(i,j,bi,bj)
enddo
enddo
enddo
enddo
c heat content of top level:
c theta * delZ * (sea water heat capacity = 3996 J/kg/K)
c * (density of sea-water = 1026 kg/m^3)
c
c heat content of sea-ice:
c tice * heff * (sea ice heat capacity = 2090 J/kg/K)
c * (density of sea-ice = 910 kg/m^3)
c
c note: to remove mass contribution to heat content,
c which is not properly accounted for by volume converving
c ocean model, theta and tice are referenced to freezing
c temperature of sea-ice, here -1.96 deg C
c
c latent heat content of sea-ice:
c - heff * (latent heat of fusion = 334000 J/kg)
c * (density of sea-ice = 910 kg/m^3)
c
c latent heat content of snow:
c - hsnow * (latent heat of fusion = 334000 J/kg)
c * (density of snow = 330 kg/m^3)
elseif ( cost_ice_flag .eq. 3 ) then
c heat content of top level plus latent heat of sea-ice
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do j = 1,sny
do i = 1,snx
objf_ice(bi,bj) = objf_ice(bi,bj) +
& tempVar * rA(i,j,bi,bj) * (
& (THETA(i,j,kSrf,bi,bj) + 1.96 _d 0 ) *
& drF(1) * 3996. _d 0 * 1026. _d 0 -
& HEFF(i,j,bi,bj) * 334000. _d 0 * 910. _d 0 )
enddo
enddo
enddo
enddo
elseif ( cost_ice_flag .eq. 4 ) then
c heat content of top level
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do j = 1,sny
do i = 1,snx
objf_ice(bi,bj) = objf_ice(bi,bj) +
& tempVar * rA(i,j,bi,bj) * (
& (THETA(i,j,kSrf,bi,bj) + 1.96 _d 0 ) *
& drF(1) * 3996. _d 0 * 1026. _d 0 )
enddo
enddo
enddo
enddo
elseif ( cost_ice_flag .eq. 5 ) then
c heat content of top level plus sea-ice plus latent heat of snow
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do j = 1,sny
do i = 1,snx
objf_ice(bi,bj) = objf_ice(bi,bj) +
& tempVar * rA(i,j,bi,bj) * (
& (THETA(i,j,kSrf,bi,bj) + 1.96 _d 0 ) *
& drF(1) * 3996. _d 0 * 1026. _d 0 +
& (TICES(i,j,1,bi,bj) - 273.15 _d 0 + 1.96 _d 0 ) *
& HEFF(i,j,bi,bj) * 2090. _d 0 * 910. _d 0 -
& HEFF(i,j,bi,bj) * 334000. _d 0 * 910. _d 0 -
& HSNOW(i,j,bi,bj) * 334000. _d 0 * 330. _d 0 )
enddo
enddo
enddo
enddo
elseif ( cost_ice_flag .eq. 6 ) then
c Qadratic cost function measuring difference between pkg/seaice
c AREA variable and simulated sea-ice measurements at every time
c step. For time being no measurements are read-in. It is
c assumed that measurements are AREA=0.5 at all times everywhere.
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do j = 1,sny
do i = 1,snx
objf_ice(bi,bj) = objf_ice(bi,bj) +
& ( AREA(i,j,bi,bj) - 0.5 _d 0 ) *
& ( AREA(i,j,bi,bj) - 0.5 _d 0 )
enddo
enddo
enddo
enddo
elseif ( cost_ice_flag .eq. 7 ) then
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do j = 1,sny
do i = 1,snx
objf_ice(bi,bj) = objf_ice(bi,bj) +
& UICE(i,j,bi,bj) * UICE(i,j,bi,bj) +
& VICE(i,j,bi,bj) * VICE(i,j,bi,bj)
enddo
enddo
enddo
enddo
else
WRITE(msgBuf,'(A)')
& 'COST_ICE: invalid cost_ice_flag'
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT , myThid )
STOP 'ABNORMAL END: S/R COST_ICE'
endif
endif
cph(
write(standardMessageUnit,*) 'ph-ice C ', myiter, objf_ice(1,1)
cph)
#endif /* ALLOW_COST_ICE */
return
end