File indexing completed on 2023-07-14 05:10:23 UTC
view on githubraw file Latest commit de57a2ec on 2023-07-13 16:55:13 UTC
8f7d13d0c9 Jean*0001 #include "ECCO_OPTIONS.h"
f1e8c1d7ee Gael*0002
0003 subroutine cost_gencost_sstv4(
f8e779c983 antn*0004 I myThid
f1e8c1d7ee Gael*0005 & )
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 implicit none
0021
0022
0023
0024 #include "EEPARAMS.h"
0025 #include "SIZE.h"
0026 #include "PARAMS.h"
0027 #include "GRID.h"
f09238ab8f Gael*0028 #ifdef ALLOW_CAL
0029 # include "cal.h"
0030 #endif
0031 #ifdef ALLOW_ECCO
13d362b8c1 Ou W*0032 # include "ECCO_SIZE.h"
0033 # include "ECCO.h"
f09238ab8f Gael*0034 #endif
ebac42f582 An T*0035 #ifdef ALLOW_SMOOTH
f09238ab8f Gael*0036 # include "SMOOTH.h"
ebac42f582 An T*0037 #endif
f1e8c1d7ee Gael*0038
0039
0040
f8e779c983 antn*0041 integer myThid
f1e8c1d7ee Gael*0042
0043 #ifdef ALLOW_SMOOTH
0044 #ifdef ALLOW_GENCOST_CONTRIBUTION
0045
0046
0047 integer bi,bj
0048 integer i,j,k
0049 integer itlo,ithi
0050 integer jtlo,jthi
0051 integer jmin,jmax
0052 integer imin,imax
0053 integer irec,jrec,krec
0054 integer ilps
0055
0056 logical doglobalread
0057 logical ladinit
0058
f09238ab8f Gael*0059 _RL mydummy
f8e779c983 antn*0060 _RL mybar(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0061 _RL anom_sst(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0062 _RL obs_sst (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0063 _RL nb_sst (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0064 _RL msk_sst (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
0065 _RL tmp_sst (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
f1e8c1d7ee Gael*0066 _RL spval
0067
0068 _RL junk,junkweight
0069
0070 integer ndaysave
0071 _RL ndaysaveRL
0072
3248fb8b5a Gael*0073 integer k2, k2_lsc
0074
de57a2ec4b Mart*0075 character*(MAX_LEN_FNAM) fname
0076 character*(MAX_LEN_FNAM) fname2
e5310f7c13 Gael*0077
0078 #ifdef ALLOW_ECCO_DEBUG
0079 character*(MAX_LEN_MBUF) msgBuf
0080 INTEGER ioUnit
0081 #endif
0082 logical exst
f1e8c1d7ee Gael*0083
0084 _RL daytime
0085 _RL diffsecs
0086 integer il, localrec
0087 integer dayiter
0088 integer daydate(4)
0089 integer difftime(4)
e36a586433 Jean*0090 integer tempDate_1
f1e8c1d7ee Gael*0091 integer middate(4)
e5310f7c13 Gael*0092 integer locstartdate(4)
f1e8c1d7ee Gael*0093 integer yday, ymod
0094 integer md, dd, sd, ld, wd
0095
f09238ab8f Gael*0096 integer kgen, kgen_lsc
b9b1c4d04c Timo*0097 logical dosumsq
f1e8c1d7ee Gael*0098
0099
0100
0101 integer ilnblnk
0102 external ilnblnk
0103
0104
0105
f8e779c983 antn*0106 jtlo = myByLo(myThid)
0107 jthi = myByHi(myThid)
0108 itlo = myBxLo(myThid)
0109 ithi = myBxHi(myThid)
f1e8c1d7ee Gael*0110 jmin = 1
f8e779c983 antn*0111 jmax = sNy
f1e8c1d7ee Gael*0112 imin = 1
f8e779c983 antn*0113 imax = sNx
f1e8c1d7ee Gael*0114
e5310f7c13 Gael*0115 #ifdef ALLOW_ECCO_DEBUG
0116 ioUnit=standardMessageUnit
0117 #endif
0118
f1e8c1d7ee Gael*0119
f09238ab8f Gael*0120 kgen=0
0121 kgen_lsc=0
3248fb8b5a Gael*0122 k2_lsc=0
f1e8c1d7ee Gael*0123 do k=1,NGENCOST
3248fb8b5a Gael*0124 if (gencost_name(k).EQ.'sstv4-amsre') kgen=k
0125 if (gencost_name(k).EQ.'sstv4-amsre-lsc') then
0126 kgen_lsc=k
0127 do k2 = 1, NGENPPROC
0128 if (gencost_posproc(k2,kgen_lsc).EQ.'smooth') k2_lsc=k2
0129 enddo
0130 endif
f1e8c1d7ee Gael*0131 enddo
0132
f09238ab8f Gael*0133 if (kgen.NE.0) then
f8e779c983 antn*0134
f09238ab8f Gael*0135
b9b1c4d04c Timo*0136 dosumsq=.TRUE.
11c3150c71 Mart*0137 call ecco_zero( gencost_weight(1-OLx,1-OLy,1,1,kgen),
0138 & 1, zeroRL, myThid )
0139 call ecco_zero( gencost_weight(1-OLx,1-OLy,1,1,kgen_lsc),
0140 & 1, zeroRL, myThid )
cbd85e4123 Gael*0141 if ( gencost_errfile(kgen) .NE. ' ' )
11c3150c71 Mart*0142 & call ecco_readwei( gencost_errfile(kgen),
0143 & gencost_weight(1-OLx,1-OLy,1,1,kgen),
0144 & 1, 1, 1, dosumsq, myThid )
cbd85e4123 Gael*0145 if ( gencost_errfile(kgen_lsc) .NE. ' ' )
11c3150c71 Mart*0146 & call ecco_readwei( gencost_errfile(kgen_lsc),
0147 & gencost_weight(1-OLx,1-OLy,1,1,kgen_lsc),
0148 & 1, 1, 1, dosumsq, myThid )
cbd85e4123 Gael*0149
f8e779c983 antn*0150 call cal_FullDate(19920101,0,locstartdate,myThid)
f09238ab8f Gael*0151
f1e8c1d7ee Gael*0152
0153 doglobalread = .false.
0154 ladinit = .false.
0155
f09238ab8f Gael*0156 ilps=ilnblnk( gencost_barfile(kgen) )
de57a2ec4b Mart*0157 write(fname,'(2a,i10.10)')
f09238ab8f Gael*0158 & gencost_barfile(kgen)(1:ilps),'.',eccoiter
f1e8c1d7ee Gael*0159
f09238ab8f Gael*0160 spval = gencost_spmin(kgen)
0161 mydummy = gencost_dummy(kgen)
f8e779c983 antn*0162
f1e8c1d7ee Gael*0163
0164
0165
0166
0167 ndaysave=7
0168 ndaysaveRL=ndaysave
0169
0170 do irec = 1, ndaysrec-ndaysave+1, 7
0171
0172 do bj = jtlo,jthi
0173 do bi = itlo,ithi
0174 do j = jmin,jmax
0175 do i = imin,imax
0176 anom_sst(i,j,bi,bj) = 0. _d 0
0177 obs_sst(i,j,bi,bj) = 0. _d 0
0178 nb_sst(i,j,bi,bj) = 0. _d 0
0179 msk_sst(i,j,bi,bj) = 0. _d 0
0180 enddo
0181 enddo
0182 enddo
0183 enddo
0184
0185
0186
0187
0188 do jrec=1,ndaysave
0189
0190 krec=irec+jrec-1
0191
0192
101f75e5cd Gael*0193 #ifdef ALLOW_AUTODIFF
f09238ab8f Gael*0194 call active_read_xy( fname, mybar, krec, doglobalread,
f8e779c983 antn*0195 & ladinit, eccoiter, myThid,
f09238ab8f Gael*0196 & mydummy )
101f75e5cd Gael*0197 #else
0198 CALL READ_REC_XY_RL( fname, mybar, kRec, 1, myThid )
0199 #endif
f1e8c1d7ee Gael*0200
0201
e5310f7c13 Gael*0202 daytime = FLOAT(secondsperday*(krec-1)) + modelstart
0203 dayiter = hoursperday*(krec-1)+modeliter0
f8e779c983 antn*0204 call cal_getdate( dayiter, daytime, daydate, myThid )
0205 call cal_convdate( daydate,yday,md,dd,sd,ld,wd,myThid )
e5310f7c13 Gael*0206 ymod = locstartdate(1)/10000
0207 if ( ymod .GE. yday ) then
e244d7989d Gael*0208 middate(1)=1
f8e779c983 antn*0209 call cal_FullDate(locstartdate(1),0,middate,myThid)
f1e8c1d7ee Gael*0210 else
e244d7989d Gael*0211 middate(1)=1
e36a586433 Jean*0212 tempDate_1 = yday*10000+100+1
f8e779c983 antn*0213 call cal_FullDate( tempDate_1, 0, middate, myThid)
f1e8c1d7ee Gael*0214 endif
f8e779c983 antn*0215 call cal_TimePassed( middate, daydate, difftime, myThid )
0216 call cal_ToSeconds( difftime, diffsecs, myThid )
e5310f7c13 Gael*0217
f1e8c1d7ee Gael*0218 localrec = int(diffsecs/86400.) + 1
0219
f09238ab8f Gael*0220 il=ilnblnk(gencost_datafile(kgen))
de57a2ec4b Mart*0221 write(fname2,'(2a,i4)')
f09238ab8f Gael*0222 & gencost_datafile(kgen)(1:il), '_', yday
e5310f7c13 Gael*0223 inquire( file=fname2, exist=exst )
0224
0225 #ifdef ALLOW_ECCO_DEBUG
0226 WRITE(msgBuf,'(A,I4,A,I4,A,I10,A,1PE15.2)') 'sstv4 reading ',
0227 & yday,' ',ymod,' ',localrec,' ',diffsecs
0228 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
0229
f8e779c983 antn*0230 CALL CAL_PRINTDATE(middate,myThid)
0231 CALL CAL_PRINTDATE(daydate,myThid)
0232 CALL CAL_PRINTDATE(difftime,myThid)
e36a586433 Jean*0233 #endif
f1e8c1d7ee Gael*0234
e244d7989d Gael*0235 if ( ( localrec .GT. 0 ).AND.(diffsecs .GT. 0.d0) ) then
f8e779c983 antn*0236 CALL READ_REC_3D_RL( fname2, cost_iprec, 1,
0237 & tmp_sst, localrec, 1, myThid )
f1e8c1d7ee Gael*0238 else
0239 do bj = jtlo,jthi
0240 do bi = itlo,ithi
0241 do j = jmin,jmax
0242 do i = imin,imax
0243 tmp_sst(i,j,bi,bj) = spval
0244 enddo
0245 enddo
0246 enddo
0247 enddo
0248 endif
0249
0250
0251 do bj = jtlo,jthi
0252 do bi = itlo,ithi
0253 do j = jmin,jmax
0254 do i = imin,imax
0255 if ( (tmp_sst(i,j,bi,bj).GT.spval).AND.
0256 & (maskc(i,j,1,bi,bj).EQ.1.) ) then
0257 anom_sst(i,j,bi,bj)= anom_sst(i,j,bi,bj)+
f09238ab8f Gael*0258 & mybar(i,j,bi,bj)-tmp_sst(i,j,bi,bj)
f1e8c1d7ee Gael*0259 obs_sst(i,j,bi,bj)= obs_sst(i,j,bi,bj)+
e36a586433 Jean*0260 & tmp_sst(i,j,bi,bj)
f1e8c1d7ee Gael*0261 nb_sst(i,j,bi,bj)=nb_sst(i,j,bi,bj)+1. _d 0
0262 endif
0263 enddo
0264 enddo
0265 enddo
0266 enddo
0267
0268 enddo
0269
0270
0271 do bj = jtlo,jthi
0272 do bi = itlo,ithi
0273 do j = jmin,jmax
0274 do i = imin,imax
0275 if ( nb_sst(i,j,bi,bj) .NE. 0. ) then
e36a586433 Jean*0276 obs_sst(i,j,bi,bj) =
f1e8c1d7ee Gael*0277 & obs_sst(i,j,bi,bj)/nb_sst(i,j,bi,bj)
e36a586433 Jean*0278 anom_sst(i,j,bi,bj) =
f1e8c1d7ee Gael*0279 & anom_sst(i,j,bi,bj)/nb_sst(i,j,bi,bj)
0280 msk_sst(i,j,bi,bj) = 1. _d 0
0281 endif
0282 enddo
0283 enddo
0284 enddo
0285 enddo
0286
0287
0288
0289
0290 #ifdef ALLOW_GENCOST_SSTV4_OUTPUT
de57a2ec4b Mart*0291 write(fname2,'(1a)') 'sstdiff_raw'
f8e779c983 antn*0292 CALL WRITE_REC_3D_RL( fname2, precFloat32, 1,
0293 & anom_sst, irec, 1, myThid )
de57a2ec4b Mart*0294 write(fname2,'(1a)') 'sstobs_raw'
f8e779c983 antn*0295 CALL WRITE_REC_3D_RL( fname2, precFloat32, 1,
0296 & obs_sst, irec, 1, myThid )
f1e8c1d7ee Gael*0297 #endif
0298
3248fb8b5a Gael*0299 if ( useSMOOTH.AND.(k2_lsc.GT.0) )
7b8b86ab99 Timo*0300 & call smooth_hetero2d(anom_sst,maskInC,
3248fb8b5a Gael*0301 & gencost_posproc_c(k2_lsc,kgen_lsc),
f8e779c983 antn*0302 & gencost_posproc_i(k2_lsc,kgen_lsc),myThid)
f1e8c1d7ee Gael*0303
0304 #ifdef ALLOW_GENCOST_SSTV4_OUTPUT
3248fb8b5a Gael*0305 if ( useSMOOTH.AND.(k2_lsc.GT.0) )
7b8b86ab99 Timo*0306 & call smooth_hetero2d(obs_sst,maskInC,
3248fb8b5a Gael*0307 & gencost_posproc_c(k2_lsc,kgen_lsc),
f8e779c983 antn*0308 & gencost_posproc_i(k2_lsc,kgen_lsc),myThid)
f1e8c1d7ee Gael*0309
de57a2ec4b Mart*0310 write(fname2,'(1a)') 'sstdiff_smooth'
f8e779c983 antn*0311 CALL WRITE_REC_3D_RL( fname2, precFloat32, 1,
0312 & anom_sst, irec, 1, myThid )
de57a2ec4b Mart*0313 write(fname2,'(1a)') 'sstobs_smooth'
f8e779c983 antn*0314 CALL WRITE_REC_3D_RL( fname2, precFloat32, 1,
0315 & obs_sst, irec, 1, myThid )
e36a586433 Jean*0316 #endif
f1e8c1d7ee Gael*0317
0318
0319
0320
0321 do bj = jtlo,jthi
0322 do bi = itlo,ithi
0323 do j = jmin,jmax
0324 do i = imin,imax
0325 junk = anom_sst(i,j,bi,bj)
f09238ab8f Gael*0326 junkweight = gencost_weight(i,j,bi,bj,kgen_lsc)*
f1e8c1d7ee Gael*0327 & maskc(i,j,1,bi,bj)
f4cfad5b95 Ou W*0328 objf_gencost(bi,bj,kgen_lsc) =
0329 & objf_gencost(bi,bj,kgen_lsc)
f1e8c1d7ee Gael*0330 & +junk*junk*junkweight/ndaysaveRL
0331 if ( (junkweight.GT.0.).AND.(nb_sst(i,j,bi,bj).GT.0.) )
f4cfad5b95 Ou W*0332 & num_gencost(bi,bj,kgen_lsc) =
0333 & num_gencost(bi,bj,kgen_lsc) + 1. _d 0 /ndaysaveRL
f1e8c1d7ee Gael*0334 enddo
0335 enddo
0336 enddo
0337 enddo
0338
0339 enddo
0340
0341
0342
0343
0344
0345 do irec = 1, ndaysrec
0346
0347
101f75e5cd Gael*0348 #ifdef ALLOW_AUTODIFF
f09238ab8f Gael*0349 call active_read_xy( fname, mybar, irec, doglobalread,
f8e779c983 antn*0350 & ladinit, eccoiter, myThid,
f09238ab8f Gael*0351 & mydummy )
101f75e5cd Gael*0352 #else
0353 CALL READ_REC_XY_RL( fname, mybar, iRec, 1, myThid )
0354 #endif
f1e8c1d7ee Gael*0355
0356
e5310f7c13 Gael*0357 daytime = FLOAT(secondsperday*(irec-1)) + modelstart
0358 dayiter = hoursperday*(irec-1)+modeliter0
f8e779c983 antn*0359 call cal_getdate( dayiter, daytime, daydate, myThid )
0360 call cal_convdate( daydate,yday,md,dd,sd,ld,wd,myThid )
e5310f7c13 Gael*0361 ymod = locstartdate(1)/10000
0362 if ( ymod .GE. yday ) then
e244d7989d Gael*0363 middate(1)=1
f8e779c983 antn*0364 call cal_FullDate(locstartdate(1),0,middate,myThid)
f1e8c1d7ee Gael*0365 else
e244d7989d Gael*0366 middate(1)=1
0367 tempDate_1 = yday*10000+100+1
f8e779c983 antn*0368 call cal_FullDate( tempDate_1, 0, middate, myThid)
f1e8c1d7ee Gael*0369 endif
f8e779c983 antn*0370 call cal_TimePassed( middate, daydate, difftime, myThid )
0371 call cal_ToSeconds( difftime, diffsecs, myThid )
e5310f7c13 Gael*0372
f1e8c1d7ee Gael*0373 localrec = int(diffsecs/86400.) + 1
0374
f09238ab8f Gael*0375 il=ilnblnk(gencost_datafile(kgen))
de57a2ec4b Mart*0376 write(fname2,'(2a,i4)')
f09238ab8f Gael*0377 & gencost_datafile(kgen)(1:il), '_', yday
e5310f7c13 Gael*0378 inquire( file=fname2, exist=exst )
0379
0380 #ifdef ALLOW_ECCO_DEBUG
0381 WRITE(msgBuf,'(A,I4,A,I4,A,I10,A,1PE15.2)') 'sstv4 reading ',
0382 & yday,' ',ymod,' ',localrec,' ',diffsecs
0383 CALL PRINT_MESSAGE( msgBuf, ioUnit, SQUEEZE_RIGHT, myThid )
0384
f8e779c983 antn*0385 CALL CAL_PRINTDATE(middate,myThid)
0386 CALL CAL_PRINTDATE(daydate,myThid)
0387 CALL CAL_PRINTDATE(difftime,myThid)
e5310f7c13 Gael*0388 #endif
f1e8c1d7ee Gael*0389
e244d7989d Gael*0390 if ( ( localrec .GT. 0 ).AND.(diffsecs .GT. 0.d0) ) then
f8e779c983 antn*0391 call READ_REC_3D_RL( fname2, cost_iprec, 1,
0392 & tmp_sst, localrec, 1, myThid )
f1e8c1d7ee Gael*0393 else
0394 do bj = jtlo,jthi
0395 do bi = itlo,ithi
0396 do j = jmin,jmax
0397 do i = imin,imax
0398 tmp_sst(i,j,bi,bj) = spval
0399 enddo
0400 enddo
0401 enddo
0402 enddo
0403 endif
0404
0405
0406 do bj = jtlo,jthi
0407 do bi = itlo,ithi
0408 do j = jmin,jmax
0409 do i = imin,imax
0410 if ( (tmp_sst(i,j,bi,bj).GT.spval).AND.
0411 & (maskc(i,j,1,bi,bj).EQ.1.) ) then
0412 anom_sst(i,j,bi,bj) =
f09238ab8f Gael*0413 & mybar(i,j,bi,bj)-tmp_sst(i,j,bi,bj)
f1e8c1d7ee Gael*0414 msk_sst(i,j,bi,bj) = 1. _d 0
0415 else
0416 anom_sst(i,j,bi,bj) = 0. _d 0
0417 msk_sst(i,j,bi,bj) = 0. _d 0
0418 endif
0419 enddo
0420 enddo
0421 enddo
0422 enddo
0423
0424 #ifdef ALLOW_GENCOST_SSTV4_OUTPUT
de57a2ec4b Mart*0425 write(fname2,'(1a)') 'sstdiff_point'
f8e779c983 antn*0426 CALL WRITE_REC_3D_RL( fname2, precFloat32, 1,
0427 & anom_sst, irec, 1, myThid )
e36a586433 Jean*0428 #endif
f1e8c1d7ee Gael*0429
0430
0431
0432 do bj = jtlo,jthi
0433 do bi = itlo,ithi
0434 do j = jmin,jmax
0435 do i = imin,imax
0436 junk = anom_sst(i,j,bi,bj)
f09238ab8f Gael*0437 junkweight = gencost_weight(i,j,bi,bj,kgen)*
f1e8c1d7ee Gael*0438 & maskc(i,j,1,bi,bj)*msk_sst(i,j,bi,bj)
f4cfad5b95 Ou W*0439 objf_gencost(bi,bj,kgen) =
0440 & objf_gencost(bi,bj,kgen)+junk*junk*junkweight
f1e8c1d7ee Gael*0441 if (junkweight.GT.0.)
f4cfad5b95 Ou W*0442 & num_gencost(bi,bj,kgen) =
0443 & num_gencost(bi,bj,kgen) + 1. _d 0
f1e8c1d7ee Gael*0444 enddo
0445 enddo
0446 enddo
0447 enddo
0448
0449 enddo
0450
f09238ab8f Gael*0451
0452 endif
0453
f1e8c1d7ee Gael*0454 #endif /* ifdef ALLOW_GENCOST_CONTRIBUTION */
0455 #endif /* ifdef ALLOW_SMOOTH */
0456
f8e779c983 antn*0457 RETURN
0458 END