File indexing completed on 2026-01-28 06:08:56 UTC
view on githubraw file Latest commit 6a9e386e on 2026-01-27 18:10:16 UTC
96b006450c dngo*0001 #include "STREAMICE_OPTIONS.h"
0002 #ifdef ALLOW_COST
0003 #include "COST_OPTIONS.h"
0004 #endif
0005
0006 subroutine streamice_cost_final( myThid )
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 IMPLICIT NONE
0018
0019
0020 #include "SIZE.h"
0021 #include "EEPARAMS.h"
0022 #include "PARAMS.h"
0023 #include "DYNVARS.h"
0024 #include "GRID.h"
0025 #ifdef ALLOW_STREAMICE
0026 # include "STREAMICE.h"
0027 #endif
0028
0029 #ifdef ALLOW_COST
0030 #include "cost.h"
0031 #endif
0032
0033
0034
0035 integer myThid
0036
0037 #ifdef ALLOW_COST_STREAMICE
0038
0039 integer bi, bj
6a9e386e2e dngo*0040 integer i, j
96b006450c dngo*0041 integer itlo,ithi
0042 integer jtlo,jthi
6a9e386e2e dngo*0043
0044
0045
0046 _RL i_numcells, dCdx, dCdy, dBdx, dBdy, idt
0047 _RL utmp, vtmp, uotmp, votmp, cotmp
96b006450c dngo*0048 _RL dRdx, dRdy, cfricval, bglenval, cfricvalp1, bglenvalp1
0049 INTEGER ILNBLNK
0050 EXTERNAL ILNBLNK
0051 CHARACTER*(MAX_LEN_FNAM) STREAMICExvelOptimFile
0052 CHARACTER*(MAX_LEN_FNAM) STREAMICEyvelOptimFile
0053 CHARACTER*(MAX_LEN_FNAM) STREAMICEerrvelOptimFile
0054 LOGICAL calc_prior_cost_bglen
0055 _RL U_obs (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0056 _RL V_obs (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0057 _RL U_err (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0058 _RL R0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0059
0060 _RL smooth_cost_fric (nSx,nSy)
feb7fa5d1e dngo*0061 _RL smooth_cost_maxmelt (nSx,nSy)
96b006450c dngo*0062 _RL smooth_cost_prior (nSx,nSy)
0063 _RL smooth_cost_bglen (nSx,nSy)
0064 _RL cost_misfit (nSx,nSy)
0065 _RL cost_misfit_norm (nSx,nSy)
0066
0067 _RL smooth_cost_fric_gl, smooth_cost_prior_gl,
0068 & smooth_cost_bglen_gl, cost_misfit_gl,
0069 & cost_thin_gl, cost_misfit_norm_gl,
feb7fa5d1e dngo*0070 & cost_vel_td_gl, smooth_cost_maxmelt_gl
96b006450c dngo*0071 _RL rhoi, rhow, r, i_r, h, hf, r_low_tmp, haf
0072
0073 jtlo = myByLo(myThid)
0074 jthi = myByHi(myThid)
0075 itlo = myBxLo(myThid)
0076 ithi = myBxHi(myThid)
0077
0078 i_numcells = 1.0/(Nx*Ny)
0079 idt = 1.0/deltaT*365.0*86400.0
0080
0081 _BARRIER
0082
0083 STREAMICExvelOptimFile=
0084 & STREAMICEvelOptimSnapBasename(1:
0085 & ILNBLNK(STREAMICEvelOptimSnapBasename))//"u.bin"
0086 STREAMICEyvelOptimFile=
0087 & STREAMICEvelOptimSnapBasename(1:
0088 & ILNBLNK(STREAMICEvelOptimSnapBasename))//"v.bin"
0089 STREAMICEerrvelOptimFile=
0090 & STREAMICEvelOptimSnapBasename(1:
0091 & ILNBLNK(STREAMICEvelOptimSnapBasename))//"err.bin"
0092
0093 if (streamice_do_snapshot_cost) then
0094 CALL READ_FLD_XY_RL( STREAMICExvelOptimFile, ' ',
0095 & U_obs, 0, myThid )
0096 print *, 'GOT HERE read velobsu'
0097 CALL READ_FLD_XY_RL( STREAMICEyvelOptimFile, ' ',
0098 & V_obs, 0, myThid )
0099 print *, 'GOT HERE read velobsv'
0100 CALL READ_FLD_XY_RL( STREAMICEerrvelOptimFile, ' ',
0101 & U_err, 0, myThid )
0102 print *, 'GOT HERE read uerr'
0103 endif
0104 if (STREAMICEBglenCostMaskFile .ne. ' ') then
0105 calc_prior_cost_bglen = .true.
0106 else
0107 calc_prior_cost_bglen = .false.
0108 endif
0109 CALL READ_FLD_XY_RL(STREAMICEtopogFile, ' ',
0110 & R0, 0, myThid )
0111 print *, 'GOT HERE read topo file'
0112
6a9e386e2e dngo*0113 _EXCH_XY_RL(R0, myThid)
96b006450c dngo*0114
0115
0116
0117
0118
0119
0120 rhoi = streamice_density
0121 rhow = streamice_density_ocean_avg
0122 r=rhoi/rhow
0123 i_r = 1/r
0124
0125 do bj = jtlo,jthi
0126 do bi = itlo,ithi
0127
0128 smooth_cost_fric (bi,bj) = 0.0
feb7fa5d1e dngo*0129 smooth_cost_maxmelt (bi,bj) = 0.0
96b006450c dngo*0130 smooth_cost_bglen (bi,bj) = 0.0
0131 smooth_cost_prior (bi,bj) = 0.0
0132 cost_misfit (bi,bj) = 0.0
0133 cost_misfit_norm (bi,bj) = 0.0
0134
0135 do j=1,sNy
0136 do i=1,sNx
0137
0138 #ifdef ALLOW_STREAMICE_TC_COST
0139 IF ( .not. STREAMICE_do_timedep_cost ) THEN
0140 #endif
0141 cfricval = C_basal_friction(i,j,bi,bj)
0142 bglenval = B_glen(i,j,bi,bj)
0143
0144 cfricvalp1 = C_basal_friction(i+1,j,bi,bj)
0145 dCdx = (Cfricvalp1-cfricval)/
0146 & (dxC(i+1,j,bi,bj))
0147 cfricvalp1 = C_basal_friction(i,j+1,bi,bj)
0148 dCdy = (Cfricvalp1-cfricval)/
0149 & (dxC(i,j+1,bi,bj))
0150
0151 bglenvalp1 = B_glen(i+1,j,bi,bj)
0152 dBdx = (bglenvalp1-
0153 & bglenval) /
0154 & (dxC(i+1,j,bi,bj))
0155 bglenvalp1 = B_glen(i,j+1,bi,bj)
0156 dBdy = (bglenvalp1-
0157 & bglenval) /
0158 & (dyC(i,j+1,bi,bj))
0159
0160 dRdx = (R_low_si(i+1,j,bi,bj)-
0161 & R_low_si(i,j,bi,bj)) /
0162 & (dxC(i+1,j,bi,bj))
0163 dRdy = (R_low_si(i,j+1,bi,bj)-
0164 & R_low_si(i,j,bi,bj)) /
0165 & (dyC(i,j+1,bi,bj))
0166
0167 #ifdef ALLOW_STREAMICE_TC_COST
0168 endif
0169 #endif
0170
0171 if (STREAMICE_do_snapshot_cost) then
0172
0173 utmp = streamice_u_surf(i,j,bi,bj)
0174 vtmp = streamice_v_surf(i,j,bi,bj)
0175 uotmp = U_obs(i,j,bi,bj)
0176 votmp = V_obs(i,j,bi,bj)
0177 cotmp = sqrt(uotmp**2+votmp**2)
0178
0179 IF((cotmp.ne.0.0) .and.
0180 & streamice_hmask(i,j,bi,bj).eq.1.0) THEN
0181
0182 tile_fc (bi,bj) = tile_fc (bi,bj) +
0183 & streamice_wgt_vel *
0184 & (
0185 & 0.5 * (streamice_u_surf(i,j,bi,bj)-
0186 & U_obs(i,j,bi,bj))**2 +
0187 & 0.5 * (streamice_v_surf(i,j,bi,bj)-
0188 & V_obs(i,j,bi,bj))**2
0189 & ) * I_numcells
0190 & / (1.0+U_err(i,j,bi,bj)**2)
0191
0192 cost_misfit (bi,bj) = cost_misfit (bi,bj) +
0193 & streamice_wgt_vel *
0194 & (
0195 & 0.5 * (streamice_u_surf(i,j,bi,bj)-
0196 & U_obs(i,j,bi,bj))**2 +
0197 & 0.5 * (streamice_v_surf(i,j,bi,bj)-
0198 & V_obs(i,j,bi,bj))**2
0199 & ) * I_numcells
0200 & / (1.0+U_err(i,j,bi,bj)**2)
0201
0202 tile_fc (bi,bj) = tile_fc (bi,bj) +
0203 & streamice_wgt_vel_norm *
0204 & (
0205 & 0.5 * (streamice_u_surf(i,j,bi,bj)-
0206 & U_obs(i,j,bi,bj))**2 +
0207 & 0.5 * (streamice_v_surf(i,j,bi,bj)-
0208 & V_obs(i,j,bi,bj))**2
0209 & ) * I_numcells
0210 & / (U_err(i,j,bi,bj)**2+
0211 & U_obs(i,j,bi,bj)**2+V_obs(i,j,bi,bj)**2)
0212
0213 cost_misfit_norm (bi,bj) = cost_misfit_norm (bi,bj) +
0214 & streamice_wgt_vel_norm *
0215 & (
0216 & 0.5 * (streamice_u_surf(i,j,bi,bj)-
0217 & U_obs(i,j,bi,bj))**2 +
0218 & 0.5 * (streamice_v_surf(i,j,bi,bj)-
0219 & V_obs(i,j,bi,bj))**2
0220 & ) * I_numcells
0221 & / (U_err(i,j,bi,bj)**2+
0222 & U_obs(i,j,bi,bj)**2+V_obs(i,j,bi,bj)**2)
0223
0224 ENDIF
0225 endif
0226
0227 #ifdef ALLOW_STREAMICE_TC_COST
0228 IF ( .not. STREAMICE_do_timedep_cost ) THEN
0229 #endif
0230
0231 if(streamice_hmask(i,j,bi,bj).eq.1.0) then
0232 tile_fc (bi,bj) = tile_fc (bi,bj) +
0233 & streamice_wgt_tikh_beta * (dCdx**2+dCdy**2) * I_numcells
0234 smooth_cost_fric (bi,bj) = smooth_cost_fric (bi,bj) +
0235 & streamice_wgt_tikh_beta * (dCdx**2+dCdy**2) * I_numcells
0236
0237 tile_fc (bi,bj) = tile_fc (bi,bj) +
0238 & streamice_wgt_tikh_bglen * (dBdx**2+dBdy**2) * I_numcells
0239 smooth_cost_bglen (bi,bj) = smooth_cost_bglen (bi,bj) +
0240 & streamice_wgt_tikh_bglen * (dBdx**2+dBdy**2) * I_numcells
0241
0242 tile_fc (bi,bj) = tile_fc (bi,bj) +
0243 & 0.e5 * (dRdx**2+dRdy**2) * I_numcells
0244
0245 h = H_streamice(i,j,bi,bj)
0246 hf = -1.0 * i_r * R_low_si (i,j,bi,bj)
0247
0248 IF ((h-hf) .gt. 5. .AND. B_glen0(i,j,bi,bj).gt.0.0 .and.
0249 & calc_prior_cost_bglen) then
0250
0251 tile_fc (bi,bj) = tile_fc (bi,bj) +
0252 & streamice_wgt_prior_bglen * I_numcells *
0253 & (B_glen(i,j,bi,bj)-B_glen0(i,j,bi,bj))**2
0254
0255 smooth_cost_prior (bi,bj) = smooth_cost_prior (bi,bj) +
0256 & streamice_wgt_prior_bglen * I_numcells *
0257 & (B_glen(i,j,bi,bj)-B_glen0(i,j,bi,bj))**2
0258
0259 ENDIF
0260 endif
0261
0262 #ifdef ALLOW_STREAMICE_TC_COST
0263 ENDIF
0264 #endif
0265
0266 if (streamice_do_verification_cost) then
0267 if (streamice_hmask(i,j,bi,bj).eq.1.0) then
0268 tile_fc (bi,bj) = tile_fc (bi,bj) +
0269 & u_streamice(i,j,bi,bj)**2+v_streamice(i,j,bi,bj)**2+
0270 & h_streamice(i,j,bi,bj)**2
0271 endif
0272 endif
0273
0274 if (streamice_do_vaf_cost) then
0275 if (streamice_hmask(i,j,bi,bj).eq.1.0) then
0276
0277 r_low_tmp = R_low_si (i,j,bi,bj)
0278
0279 if (r_low_tmp .ge. 0.0) then
0280 HAF = h_streamice (i,j,bi,bj)
0281 else
0282 HAF = h_streamice (i,j,bi,bj) +
0283 & (streamice_density_ocean_avg / streamice_density
0284 & * r_low_tmp)
0285 endif
0286
0287 if (HAF .gt. 0.0) then
0288 tile_fc(bi,bj) = tile_fc(bi,bj) +
0289 & HAF * rA (i,j,bi,bj)
0290 endif
0291
0292 endif
0293 endif
0294
0295 end do
0296 end do
0297
0298 tile_fc (bi,bj) =
0299 & tile_fc (bi,bj) + cost_func1_streamice(bi,bj)
0300
0301 #ifdef ALLOW_STREAMICE_TC_COST
0302 if ( STREAMICE_do_timedep_cost ) then
0303 smooth_cost_prior(bi,bj) = cost_prior_streamice(bi,bj)
0304 smooth_cost_fric(bi,bj) = cost_smooth_fric_streamice(bi,bj)
0305 smooth_cost_bglen(bi,bj) = cost_smooth_glen_streamice(bi,bj)
feb7fa5d1e dngo*0306 smooth_cost_maxmelt(bi,bj) =
0307 & cost_smooth_bdotmelt_streamice(bi,bj)
96b006450c dngo*0308 endif
0309 #endif
0310
0311 end do
0312 end do
0313
0314 CALL GLOBAL_SUM_TILE_RL
0315 & ( smooth_cost_fric, smooth_cost_fric_gl, myThid )
0316 CALL GLOBAL_SUM_TILE_RL
0317 & ( smooth_cost_bglen, smooth_cost_bglen_gl, myThid )
0318 CALL GLOBAL_SUM_TILE_RL
0319 & ( smooth_cost_prior, smooth_cost_prior_gl, myThid )
0320 CALL GLOBAL_SUM_TILE_RL
0321 & ( cost_misfit, cost_misfit_gl, myThid )
0322 CALL GLOBAL_SUM_TILE_RL
0323 & ( cost_misfit_norm, cost_misfit_norm_gl, myThid )
0324 CALL GLOBAL_SUM_TILE_RL
0325 & ( cost_surf_streamice, cost_thin_gl, myThid )
0326 CALL GLOBAL_SUM_TILE_RL
0327 & ( cost_vel_streamice, cost_vel_td_gl, myThid )
feb7fa5d1e dngo*0328 CALL GLOBAL_SUM_TILE_RL
0329 & ( smooth_cost_maxmelt, smooth_cost_maxmelt_gl, myThid )
0330 write(standardMessageUnit,'(A,D22.15)') 'maxmlt smooth contr = ',
0331 & smooth_cost_maxmelt_gl
96b006450c dngo*0332
0333 if (STREAMICE_do_timedep_cost .or.
0334 & STREAMICE_do_snapshot_cost) then
0335
0336 write(standardMessageUnit,'(A,D22.15)') 'fric smooth contr = ',
0337 & smooth_cost_fric_gl
0338 write(standardMessageUnit,'(A,D22.15)') 'bglen smooth contr = ',
0339 & smooth_cost_bglen_gl
0340 write(standardMessageUnit,'(A,D22.15)') 'prior smooth contr = ',
0341 & smooth_cost_prior_gl
0342 write(standardMessageUnit,'(A,D22.15)') 'vel misfit = ',
0343 & cost_misfit_gl
0344 write(standardMessageUnit,'(A,D22.15)') 'vel misfit norm = ',
0345 & cost_misfit_norm_gl
0346
0347 endif
0348
0349 if (STREAMICE_do_timedep_cost) then
0350
0351 write(standardMessageUnit,'(A,D22.15)') 'thinning contr = ',
0352 & cost_thin_gl
0353 write(standardMessageUnit,'(A,D22.15)') 'td vel misfit = ',
0354 & cost_vel_td_gl
0355
0356 endif
0357
0358 #endif
0359
0360 RETURN
0361 END
