pkg/seaice/seaice_cost_test.F

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