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