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