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
69361556c2 Mart*0045       IMPLICIT NONE
5001c65f45 Patr*0046
69361556c2 Mart*0047 #if (defined ALLOW_COST && defined ALLOW_COST_ICE)
5001c65f45 Patr*0048 c     == global variables ==
                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"
69361556c2 Mart*0058 #endif /* ALLOW_COST & ALLOW_COST_ICE */
5001c65f45 Patr*0059
                0060 c     == routine arguments ==
                0061       _RL     mytime
                0062       integer myiter
                0063       integer mythid
                0064
69361556c2 Mart*0065 #if (defined ALLOW_COST && defined ALLOW_COST_ICE)
                0066 c     == external functions ==
                0067       integer  ilnblnk
                0068       external ilnblnk
5001c65f45 Patr*0069
                0070 c     == local variables ==
                0071 c     msgBuf      - Informational/error message buffer
                0072       CHARACTER*(MAX_LEN_MBUF) msgBuf
ae6832360b Mart*0073       integer bi,bj,i,j,kSrf
5001c65f45 Patr*0074       _RL tempVar
                0075
                0076 c     == end of interface ==
                0077
0320e25227 Mart*0078       if ( usingPCoords ) then
                0079        kSrf = Nr
                0080       else
                0081        kSrf = 1
                0082       endif
3889aa6caa Patr*0083       if ( myTime .GT. (endTime - lastinterval) ) then
989b1a2fcf Jean*0084          tempVar = 1. _d 0/
                0085      &             ( ( 1. _d 0 + min(endTime-startTime,lastinterval) )
7c7521a1da Jean*0086      &             / deltaTClock )
5001c65f45 Patr*0087
3ad0d94cb0 Patr*0088 cph(
0320e25227 Mart*0089       write(standardMessageUnit,*) 'ph-ice B ', myiter,
344ddc3242 Mart*0090      &        theta(4,4,kSrf,1,1), area(4,4,1,1), heff(4,4,1,1)
3ad0d94cb0 Patr*0091 cph)
5001c65f45 Patr*0092          if ( cost_ice_flag .eq. 1 ) then
                0093 c     sea-ice volume
                0094             do bj=myByLo(myThid),myByHi(myThid)
                0095                do bi=myBxLo(myThid),myBxHi(myThid)
                0096                   do j = 1,sny
                0097                      do i =  1,snx
                0098                         objf_ice(bi,bj) = objf_ice(bi,bj) +
f7d3a281ce Mart*0099      &                       tempVar * rA(i,j,bi,bj) * HEFF(i,j,bi,bj)
5001c65f45 Patr*0100                      enddo
                0101                   enddo
                0102                enddo
                0103             enddo
                0104
                0105          elseif ( cost_ice_flag .eq. 2 ) then
                0106 c     sea-ice area
                0107             do bj=myByLo(myThid),myByHi(myThid)
                0108                do bi=myBxLo(myThid),myBxHi(myThid)
                0109                   do j = 1,sny
                0110                      do i =  1,snx
                0111                         objf_ice(bi,bj) = objf_ice(bi,bj) +
f7d3a281ce Mart*0112      &                       tempVar * rA(i,j,bi,bj) * AREA(i,j,bi,bj)
5001c65f45 Patr*0113                      enddo
                0114                   enddo
                0115                enddo
                0116             enddo
                0117
                0118 c     heat content of top level:
                0119 c     theta * delZ * (sea water heat capacity = 3996 J/kg/K)
                0120 c                  * (density of sea-water = 1026 kg/m^3)
                0121 c
                0122 c     heat content of sea-ice:
                0123 c     tice * heff * (sea ice heat capacity = 2090 J/kg/K)
                0124 c                 * (density of sea-ice = 910 kg/m^3)
                0125 c
                0126 c     note: to remove mass contribution to heat content,
                0127 c     which is not properly accounted for by volume converving
                0128 c     ocean model, theta and tice are referenced to freezing
                0129 c     temperature of sea-ice, here -1.96 deg C
                0130 c
                0131 c     latent heat content of sea-ice:
                0132 c     - heff * (latent heat of fusion = 334000 J/kg)
                0133 c            * (density of sea-ice = 910 kg/m^3)
                0134 c
                0135 c     latent heat content of snow:
                0136 c     - hsnow * (latent heat of fusion = 334000 J/kg)
                0137 c             * (density of snow = 330 kg/m^3)
                0138
                0139          elseif ( cost_ice_flag .eq. 3 ) then
                0140 c     heat content of top level plus latent heat of sea-ice
                0141             do bj=myByLo(myThid),myByHi(myThid)
                0142              do bi=myBxLo(myThid),myBxHi(myThid)
                0143               do j = 1,sny
                0144                do i =  1,snx
                0145                 objf_ice(bi,bj) = objf_ice(bi,bj) +
                0146      &                 tempVar * rA(i,j,bi,bj) * (
989b1a2fcf Jean*0147      &                 (THETA(i,j,kSrf,bi,bj) + 1.96 _d 0 ) *
0320e25227 Mart*0148      &                 drF(kSrf) * 3996. _d 0 * 1026. _d 0 -
989b1a2fcf Jean*0149      &                 HEFF(i,j,bi,bj) * 334000. _d 0 * 910. _d 0 )
5001c65f45 Patr*0150                enddo
                0151               enddo
                0152              enddo
                0153             enddo
                0154
                0155          elseif ( cost_ice_flag .eq. 4 ) then
                0156 c     heat content of top level
                0157             do bj=myByLo(myThid),myByHi(myThid)
                0158              do bi=myBxLo(myThid),myBxHi(myThid)
                0159               do j = 1,sny
                0160                do i =  1,snx
                0161                 objf_ice(bi,bj) = objf_ice(bi,bj) +
                0162      &                 tempVar * rA(i,j,bi,bj) * (
989b1a2fcf Jean*0163      &                 (THETA(i,j,kSrf,bi,bj) + 1.96 _d 0 ) *
0320e25227 Mart*0164      &                 drF(kSrf) * 3996. _d 0 * 1026. _d 0 )
5001c65f45 Patr*0165                enddo
                0166               enddo
                0167              enddo
                0168             enddo
                0169
                0170          elseif ( cost_ice_flag .eq. 5 ) then
                0171 c     heat content of top level plus sea-ice plus latent heat of snow
                0172             do bj=myByLo(myThid),myByHi(myThid)
                0173              do bi=myBxLo(myThid),myBxHi(myThid)
                0174               do j = 1,sny
                0175                do i =  1,snx
                0176                 objf_ice(bi,bj) = objf_ice(bi,bj) +
                0177      &                 tempVar * rA(i,j,bi,bj) * (
989b1a2fcf Jean*0178      &                 (THETA(i,j,kSrf,bi,bj) + 1.96 _d 0 ) *
0320e25227 Mart*0179      &                 drF(kSrf) * 3996. _d 0 * 1026. _d 0 +
6e5facdf0e Mart*0180      &                 (TICES(i,j,1,bi,bj) - 273.15 _d 0 + 1.96 _d 0 ) *
989b1a2fcf Jean*0181      &                 HEFF(i,j,bi,bj) * 2090. _d 0 * 910. _d 0 -
                0182      &                 HEFF(i,j,bi,bj) * 334000. _d 0 * 910. _d 0 -
                0183      &                 HSNOW(i,j,bi,bj) * 334000. _d 0 * 330. _d 0 )
5001c65f45 Patr*0184                enddo
                0185               enddo
                0186              enddo
                0187             enddo
                0188
                0189          elseif ( cost_ice_flag .eq. 6 ) then
                0190 c     Qadratic cost function measuring difference between pkg/seaice
                0191 c     AREA variable and simulated sea-ice measurements at every time
                0192 c     step.  For time being no measurements are read-in.  It is
                0193 c     assumed that measurements are AREA=0.5 at all times everywhere.
                0194             do bj=myByLo(myThid),myByHi(myThid)
                0195                do bi=myBxLo(myThid),myBxHi(myThid)
                0196                   do j = 1,sny
                0197                      do i =  1,snx
                0198                         objf_ice(bi,bj) = objf_ice(bi,bj) +
989b1a2fcf Jean*0199      &                       ( AREA(i,j,bi,bj) - 0.5 _d 0 ) *
                0200      &                       ( AREA(i,j,bi,bj) - 0.5 _d 0 )
5001c65f45 Patr*0201                      enddo
                0202                   enddo
                0203                enddo
                0204             enddo
                0205
d877a5eaeb Patr*0206          elseif ( cost_ice_flag .eq. 7 ) then
                0207
                0208             do bj=myByLo(myThid),myByHi(myThid)
                0209                do bi=myBxLo(myThid),myBxHi(myThid)
                0210                   do j = 1,sny
                0211                      do i =  1,snx
                0212                         objf_ice(bi,bj) = objf_ice(bi,bj) +
                0213      &                       UICE(i,j,bi,bj) * UICE(i,j,bi,bj) +
                0214      &                       VICE(i,j,bi,bj) * VICE(i,j,bi,bj)
                0215
                0216                      enddo
                0217                   enddo
                0218                enddo
                0219             enddo
                0220
5001c65f45 Patr*0221          else
                0222             WRITE(msgBuf,'(A)')
                0223      &           'COST_ICE: invalid cost_ice_flag'
                0224             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0225      &           SQUEEZE_RIGHT , myThid )
                0226             STOP 'ABNORMAL END: S/R COST_ICE'
                0227          endif
                0228       endif
                0229
3874013cca Patr*0230 cph(
344ddc3242 Mart*0231       write(standardMessageUnit,*) 'ph-ice C ', myiter, objf_ice(1,1)
3874013cca Patr*0232 cph)
                0233
69361556c2 Mart*0234 #endif /* ALLOW_COST & ALLOW_COST_ICE */
5001c65f45 Patr*0235
69361556c2 Mart*0236       RETURN
                0237       END