pkg/seaice/seaice_cost_test.F

869864d4b6 Patr*0001 #include "SEAICE_OPTIONS.h"
5001c65f45 Patr*0002
2b959ba38e Mart*0003       subroutine seaice_cost_test( mytime, myiter, mythid )
5001c65f45 Patr*0004
                0005 c     ==================================================================
2b959ba38e Mart*0006 c     SUBROUTINE seaice_cost_test
5001c65f45 Patr*0007 c     ==================================================================
                0008 c
465da1ecf8 Dimi*0009 c     o Compute sea-ice cost function.  The following options can be
                0010 c       selected with data.seaice (SEAICE_PARM02) variable cost_ice_flag
5001c65f45 Patr*0011 c
                0012 c     cost_ice_flag = 1
                0013 c     - compute mean sea-ice volume
                0014 c       costIceStart < mytime < costIceEnd
                0015 c
                0016 c     cost_ice_flag = 2
                0017 c     - compute mean sea-ice area
                0018 c       costIceStart < mytime < costIceEnd
                0019 c
                0020 c     cost_ice_flag = 3
                0021 c     - heat content of top level plus latent heat of sea-ice
                0022 c       costIceStart < mytime < costIceEnd
                0023 c
                0024 c     cost_ice_flag = 4
                0025 c     - heat content of top level
                0026 c       costIceStart < mytime < costIceEnd
                0027 c
                0028 c     cost_ice_flag = 5
                0029 c     - heat content of top level plus sea-ice plus latent heat of snow
                0030 c       costIceStart < mytime < costIceEnd
                0031 c
                0032 c     cost_ice_flag = 6
                0033 c     - quadratic cost function measuring difference between pkg/seaice
                0034 c       AREA variable and simulated sea-ice measurements at every time
                0035 c       step.
                0036 c
                0037 c     ==================================================================
                0038 c
                0039 c     started: menemenlis@jpl.nasa.gov 26-Feb-2003
                0040 c
                0041 c     ==================================================================
2b959ba38e Mart*0042 c     SUBROUTINE seaice_cost_test
5001c65f45 Patr*0043 c     ==================================================================
                0044
                0045       implicit none
                0046
                0047 c     == global variables ==
                0048 #ifdef ALLOW_COST_ICE
                0049 #include "EEPARAMS.h"
                0050 #include "SIZE.h"
                0051 #include "GRID.h"
                0052 #include "PARAMS.h"
03c669d1ab Jean*0053 #include "SEAICE_SIZE.h"
869864d4b6 Patr*0054 #include "SEAICE_COST.h"
5001c65f45 Patr*0055 #include "SEAICE.h"
                0056 #include "DYNVARS.h"
869864d4b6 Patr*0057 #include "cost.h"
5001c65f45 Patr*0058 #endif /* ALLOW_COST_ICE */
                0059
                0060 c     == routine arguments ==
                0061
                0062       _RL     mytime
                0063       integer myiter
                0064       integer mythid
                0065
e4775240e5 Dimi*0066 #ifdef ALLOW_COST_ICE
5001c65f45 Patr*0067
                0068 c     == local variables ==
                0069
                0070 c     msgBuf      - Informational/error message buffer
                0071       CHARACTER*(MAX_LEN_MBUF) msgBuf
ae6832360b Mart*0072       integer bi,bj,i,j,kSrf
5001c65f45 Patr*0073       _RL tempVar
                0074
                0075 c     == external functions ==
                0076
                0077       integer  ilnblnk
                0078       external ilnblnk
                0079
                0080 c     == end of interface ==
                0081
0320e25227 Mart*0082       if ( usingPCoords ) then
                0083        kSrf = Nr
                0084       else
                0085        kSrf = 1
                0086       endif
3889aa6caa Patr*0087       if ( myTime .GT. (endTime - lastinterval) ) then
989b1a2fcf Jean*0088          tempVar = 1. _d 0/
                0089      &             ( ( 1. _d 0 + min(endTime-startTime,lastinterval) )
7c7521a1da Jean*0090      &             / deltaTClock )
5001c65f45 Patr*0091
3ad0d94cb0 Patr*0092 cph(
0320e25227 Mart*0093       write(standardMessageUnit,*) 'ph-ice B ', myiter,
344ddc3242 Mart*0094      &        theta(4,4,kSrf,1,1), area(4,4,1,1), heff(4,4,1,1)
3ad0d94cb0 Patr*0095 cph)
5001c65f45 Patr*0096          if ( cost_ice_flag .eq. 1 ) then
                0097 c     sea-ice volume
                0098             do bj=myByLo(myThid),myByHi(myThid)
                0099                do bi=myBxLo(myThid),myBxHi(myThid)
                0100                   do j = 1,sny
                0101                      do i =  1,snx
                0102                         objf_ice(bi,bj) = objf_ice(bi,bj) +
f7d3a281ce Mart*0103      &                       tempVar * rA(i,j,bi,bj) * HEFF(i,j,bi,bj)
5001c65f45 Patr*0104                      enddo
                0105                   enddo
                0106                enddo
                0107             enddo
                0108
                0109          elseif ( cost_ice_flag .eq. 2 ) then
                0110 c     sea-ice area
                0111             do bj=myByLo(myThid),myByHi(myThid)
                0112                do bi=myBxLo(myThid),myBxHi(myThid)
                0113                   do j = 1,sny
                0114                      do i =  1,snx
                0115                         objf_ice(bi,bj) = objf_ice(bi,bj) +
f7d3a281ce Mart*0116      &                       tempVar * rA(i,j,bi,bj) * AREA(i,j,bi,bj)
5001c65f45 Patr*0117                      enddo
                0118                   enddo
                0119                enddo
                0120             enddo
                0121
                0122 c     heat content of top level:
                0123 c     theta * delZ * (sea water heat capacity = 3996 J/kg/K)
                0124 c                  * (density of sea-water = 1026 kg/m^3)
                0125 c
                0126 c     heat content of sea-ice:
                0127 c     tice * heff * (sea ice heat capacity = 2090 J/kg/K)
                0128 c                 * (density of sea-ice = 910 kg/m^3)
                0129 c
                0130 c     note: to remove mass contribution to heat content,
                0131 c     which is not properly accounted for by volume converving
                0132 c     ocean model, theta and tice are referenced to freezing
                0133 c     temperature of sea-ice, here -1.96 deg C
                0134 c
                0135 c     latent heat content of sea-ice:
                0136 c     - heff * (latent heat of fusion = 334000 J/kg)
                0137 c            * (density of sea-ice = 910 kg/m^3)
                0138 c
                0139 c     latent heat content of snow:
                0140 c     - hsnow * (latent heat of fusion = 334000 J/kg)
                0141 c             * (density of snow = 330 kg/m^3)
                0142
                0143          elseif ( cost_ice_flag .eq. 3 ) then
                0144 c     heat content of top level plus latent heat of sea-ice
                0145             do bj=myByLo(myThid),myByHi(myThid)
                0146              do bi=myBxLo(myThid),myBxHi(myThid)
                0147               do j = 1,sny
                0148                do i =  1,snx
                0149                 objf_ice(bi,bj) = objf_ice(bi,bj) +
                0150      &                 tempVar * rA(i,j,bi,bj) * (
989b1a2fcf Jean*0151      &                 (THETA(i,j,kSrf,bi,bj) + 1.96 _d 0 ) *
0320e25227 Mart*0152      &                 drF(kSrf) * 3996. _d 0 * 1026. _d 0 -
989b1a2fcf Jean*0153      &                 HEFF(i,j,bi,bj) * 334000. _d 0 * 910. _d 0 )
5001c65f45 Patr*0154                enddo
                0155               enddo
                0156              enddo
                0157             enddo
                0158
                0159          elseif ( cost_ice_flag .eq. 4 ) then
                0160 c     heat content of top level
                0161             do bj=myByLo(myThid),myByHi(myThid)
                0162              do bi=myBxLo(myThid),myBxHi(myThid)
                0163               do j = 1,sny
                0164                do i =  1,snx
                0165                 objf_ice(bi,bj) = objf_ice(bi,bj) +
                0166      &                 tempVar * rA(i,j,bi,bj) * (
989b1a2fcf Jean*0167      &                 (THETA(i,j,kSrf,bi,bj) + 1.96 _d 0 ) *
0320e25227 Mart*0168      &                 drF(kSrf) * 3996. _d 0 * 1026. _d 0 )
5001c65f45 Patr*0169                enddo
                0170               enddo
                0171              enddo
                0172             enddo
                0173
                0174          elseif ( cost_ice_flag .eq. 5 ) then
                0175 c     heat content of top level plus sea-ice plus latent heat of snow
                0176             do bj=myByLo(myThid),myByHi(myThid)
                0177              do bi=myBxLo(myThid),myBxHi(myThid)
                0178               do j = 1,sny
                0179                do i =  1,snx
                0180                 objf_ice(bi,bj) = objf_ice(bi,bj) +
                0181      &                 tempVar * rA(i,j,bi,bj) * (
989b1a2fcf Jean*0182      &                 (THETA(i,j,kSrf,bi,bj) + 1.96 _d 0 ) *
0320e25227 Mart*0183      &                 drF(kSrf) * 3996. _d 0 * 1026. _d 0 +
6e5facdf0e Mart*0184      &                 (TICES(i,j,1,bi,bj) - 273.15 _d 0 + 1.96 _d 0 ) *
989b1a2fcf Jean*0185      &                 HEFF(i,j,bi,bj) * 2090. _d 0 * 910. _d 0 -
                0186      &                 HEFF(i,j,bi,bj) * 334000. _d 0 * 910. _d 0 -
                0187      &                 HSNOW(i,j,bi,bj) * 334000. _d 0 * 330. _d 0 )
5001c65f45 Patr*0188                enddo
                0189               enddo
                0190              enddo
                0191             enddo
                0192
                0193          elseif ( cost_ice_flag .eq. 6 ) then
                0194 c     Qadratic cost function measuring difference between pkg/seaice
                0195 c     AREA variable and simulated sea-ice measurements at every time
                0196 c     step.  For time being no measurements are read-in.  It is
                0197 c     assumed that measurements are AREA=0.5 at all times everywhere.
                0198             do bj=myByLo(myThid),myByHi(myThid)
                0199                do bi=myBxLo(myThid),myBxHi(myThid)
                0200                   do j = 1,sny
                0201                      do i =  1,snx
                0202                         objf_ice(bi,bj) = objf_ice(bi,bj) +
989b1a2fcf Jean*0203      &                       ( AREA(i,j,bi,bj) - 0.5 _d 0 ) *
                0204      &                       ( AREA(i,j,bi,bj) - 0.5 _d 0 )
5001c65f45 Patr*0205                      enddo
                0206                   enddo
                0207                enddo
                0208             enddo
                0209
d877a5eaeb Patr*0210          elseif ( cost_ice_flag .eq. 7 ) then
                0211
                0212             do bj=myByLo(myThid),myByHi(myThid)
                0213                do bi=myBxLo(myThid),myBxHi(myThid)
                0214                   do j = 1,sny
                0215                      do i =  1,snx
                0216                         objf_ice(bi,bj) = objf_ice(bi,bj) +
                0217      &                       UICE(i,j,bi,bj) * UICE(i,j,bi,bj) +
                0218      &                       VICE(i,j,bi,bj) * VICE(i,j,bi,bj)
                0219
                0220                      enddo
                0221                   enddo
                0222                enddo
                0223             enddo
                0224
5001c65f45 Patr*0225          else
                0226             WRITE(msgBuf,'(A)')
                0227      &           'COST_ICE: invalid cost_ice_flag'
                0228             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0229      &           SQUEEZE_RIGHT , myThid )
                0230             STOP 'ABNORMAL END: S/R COST_ICE'
                0231          endif
                0232       endif
                0233
3874013cca Patr*0234 cph(
344ddc3242 Mart*0235       write(standardMessageUnit,*) 'ph-ice C ', myiter, objf_ice(1,1)
3874013cca Patr*0236 cph)
                0237
5001c65f45 Patr*0238 #endif /* ALLOW_COST_ICE */
                0239
989b1a2fcf Jean*0240       return
5001c65f45 Patr*0241       end