Changeset 12523
- Timestamp:
- 2020-03-09T11:59:47+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/releases/release-4.0-HEAD/src/OCE/DIA/diaharm.F90
r11536 r12523 24 24 25 25 INTEGER, PARAMETER :: jpincomax = 2.*jpmax_harmo 26 INTEGER, PARAMETER :: jpdimsparse = jpincomax*3 00*2426 INTEGER, PARAMETER :: jpdimsparse = jpincomax*366*24*2 ! 30min for a 1yr-long run 27 27 28 28 ! !!** namelist variables ** … … 35 35 INTEGER , ALLOCATABLE, DIMENSION(:) :: name 36 36 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp 37 REAL(wp), ALLOCATABLE, DIMENSION(:) :: ana_freq, ut , vt , ft 38 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: out_eta , out_u, out_v 39 40 INTEGER :: ninco, nsparse 37 REAL(wp), ALLOCATABLE, DIMENSION(:) :: ana_freq, ut, vt, ft 38 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: out_eta, out_u, out_v 39 40 INTEGER :: ninco, nsparse 41 REAL(wp) :: z1_tmp3 41 42 INTEGER , DIMENSION(jpdimsparse) :: njsparse, nisparse 42 43 INTEGER , SAVE, DIMENSION(jpincomax) :: ipos1 43 44 REAL(wp), DIMENSION(jpdimsparse) :: valuesparse 44 REAL(wp), DIMENSION(jpincomax) :: ztmp4 , ztmp7 45 REAL(wp), DIMENSION(jpincomax) :: ztmp4 , ztmp7, z1_pivot 45 46 REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) :: ztmp3 , zpilier 46 REAL(wp), SAVE, DIMENSION(jpincomax) :: zpivot47 47 48 48 CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: tname ! Names of tidal constituents ('M2', 'K1',...) … … 76 76 WRITE(numout,*) 77 77 WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization' 78 WRITE(numout,*) '~~~~~~~ '78 WRITE(numout,*) '~~~~~~~~~~~~~ ' 79 79 ENDIF 80 80 ! … … 131 131 ENDIF 132 132 133 ALLOCATE(name 133 ALLOCATE(name(nb_ana)) 134 134 DO jh=1,nb_ana 135 135 DO ji=1,jpmax_harmo … … 163 163 164 164 165 SUBROUTINE dia_harm 165 SUBROUTINE dia_harm( kt ) 166 166 !!---------------------------------------------------------------------- 167 167 !! *** ROUTINE dia_harm *** … … 172 172 !! 173 173 !!-------------------------------------------------------------------- 174 INTEGER, INTENT( IN ) ::kt175 ! 176 INTEGER :: ji, jj, jh, jc, nhc177 REAL(wp) :: ztime, ztemp174 INTEGER, INTENT( in ) :: kt 175 ! 176 INTEGER :: ji, jj, jh, jc, nhc 177 REAL(wp) :: ztime, ztemp 178 178 !!-------------------------------------------------------------------- 179 179 IF( ln_timing ) CALL timing_start('dia_harm') … … 187 187 DO jc = 1, 2 188 188 nhc = nhc+1 189 ztemp = (MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) &189 ztemp = ( MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) & 190 190 & +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh))) 191 191 ! 192 DO jj = 1,jpj 193 DO ji = 1,jpi 194 ! Elevation 195 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj) 196 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj) 197 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj) 192 DO jj = 2, jpjm1 193 DO ji = 2, jpim1 194 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp * sshn(ji,jj) * ssmask (ji,jj) ! elevation 195 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp * un_b(ji,jj) * ssumask(ji,jj) ! u-vel 196 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp * vn_b(ji,jj) * ssvmask(ji,jj) ! v-vel 198 197 END DO 199 198 END DO 200 ! 201 END DO 202 END DO 203 ! 199 END DO 200 END DO 204 201 END IF 205 202 ! … … 220 217 !! 221 218 !!-------------------------------------------------------------------- 222 INTEGER :: ji, jj, jh, jc, jn, nhan, jl 223 INTEGER :: ksp, kun, keq 224 REAL(wp) :: ztime, ztime_ini, ztime_end 225 REAL(wp) :: X1, X2 226 REAL(wp), DIMENSION(jpi,jpj,jpmax_harmo,2) :: ana_amp ! workspace 219 INTEGER :: ji, jj, jh, jc, jn, nhan 220 INTEGER :: ksp, kun, keq 221 REAL(wp) :: ztime, ztime_ini, ztime_end, z1_han 227 222 !!-------------------------------------------------------------------- 228 223 ! 229 224 IF(lwp) WRITE(numout,*) 230 IF(lwp) WRITE(numout,*) ' anharmo_end: kt=nitend_han: Perform harmonic analysis'225 IF(lwp) WRITE(numout,*) 'dia_harm_end: kt=nitend_han: Perform harmonic analysis' 231 226 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 227 228 ALLOCATE( out_eta(jpi,jpj,2*nb_ana), out_u(jpi,jpj,2*nb_ana), out_v(jpi,jpj,2*nb_ana) ) 232 229 233 230 ztime_ini = nit000_han*rdt ! Initial time in seconds at the beginning of analysis 234 231 ztime_end = nitend_han*rdt ! Final time in seconds at the end of analysis 235 232 nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis 236 233 z1_han = 1._wp / REAL(nhan-1) 234 237 235 ninco = 2*nb_ana 238 236 … … 240 238 keq = 0 241 239 DO jn = 1, nhan 242 ztime=( (nhan-jn)*ztime_ini + (jn-1)*ztime_end ) /FLOAT(nhan-1)240 ztime=( (nhan-jn)*ztime_ini + (jn-1)*ztime_end ) * z1_han 243 241 keq = keq + 1 244 242 kun = 0 … … 257 255 nsparse = ksp 258 256 257 IF( nsparse > jpdimsparse ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse') 258 IF( ninco > jpincomax ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax') 259 260 CALL SUR_DETERMINE_INIT 261 259 262 ! Elevation: 260 DO jj = 1, jpj 261 DO ji = 1, jpi 263 DO jj = 2, jpjm1 264 DO ji = 2, jpim1 265 262 266 ! Fill input array 263 kun = 0 264 DO jh = 1, nb_ana 265 DO jc = 1, 2 266 kun = kun + 1 267 ztmp4(kun)=ana_temp(ji,jj,kun,1) 268 END DO 269 END DO 270 271 CALL SUR_DETERMINE(jj) 272 267 ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,1) 268 CALL SUR_DETERMINE 269 273 270 ! Fill output array 274 271 DO jh = 1, nb_ana 275 ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1) 276 ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2) 277 END DO 278 END DO 279 END DO 280 281 ALLOCATE( out_eta(jpi,jpj,2*nb_ana), & 282 & out_u (jpi,jpj,2*nb_ana), & 283 & out_v (jpi,jpj,2*nb_ana) ) 284 285 DO jj = 1, jpj 286 DO ji = 1, jpi 287 DO jh = 1, nb_ana 288 X1 = ana_amp(ji,jj,jh,1) 289 X2 =-ana_amp(ji,jj,jh,2) 290 out_eta(ji,jj,jh ) = X1 * tmask_i(ji,jj) 291 out_eta(ji,jj,jh+nb_ana) = X2 * tmask_i(ji,jj) 272 out_eta(ji,jj,jh ) = ztmp7((jh-1)*2+1) * ssmask(ji,jj) 273 out_eta(ji,jj,jh+nb_ana) = -ztmp7((jh-1)*2+2) * ssmask(ji,jj) 292 274 END DO 293 275 END DO … … 295 277 296 278 ! ubar: 297 DO jj = 1, jpj 298 DO ji = 1, jpi 279 DO jj = 2, jpjm1 280 DO ji = 2, jpim1 281 299 282 ! Fill input array 300 kun=0 301 DO jh = 1,nb_ana 302 DO jc = 1,2 303 kun = kun + 1 304 ztmp4(kun)=ana_temp(ji,jj,kun,2) 305 END DO 306 END DO 307 308 CALL SUR_DETERMINE(jj+1) 283 ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,2) 284 CALL SUR_DETERMINE 309 285 310 286 ! Fill output array 311 287 DO jh = 1, nb_ana 312 ana_amp(ji,jj,jh,1) = ztmp7((jh-1)*2+1) 313 ana_amp(ji,jj,jh,2) = ztmp7((jh-1)*2+2) 314 END DO 315 316 END DO 317 END DO 318 319 DO jj = 1, jpj 320 DO ji = 1, jpi 321 DO jh = 1, nb_ana 322 X1= ana_amp(ji,jj,jh,1) 323 X2=-ana_amp(ji,jj,jh,2) 324 out_u(ji,jj, jh) = X1 * ssumask(ji,jj) 325 out_u(ji,jj,nb_ana+jh) = X2 * ssumask(ji,jj) 326 ENDDO 327 ENDDO 328 ENDDO 288 out_u(ji,jj, jh) = ztmp7((jh-1)*2+1) * ssumask(ji,jj) 289 out_u(ji,jj,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssumask(ji,jj) 290 END DO 291 292 END DO 293 END DO 329 294 330 295 ! vbar: 331 DO jj = 1, jpj 332 DO ji = 1, jpi 296 DO jj = 2, jpjm1 297 DO ji = 2, jpim1 298 333 299 ! Fill input array 334 kun=0 335 DO jh = 1,nb_ana 336 DO jc = 1,2 337 kun = kun + 1 338 ztmp4(kun)=ana_temp(ji,jj,kun,3) 339 END DO 340 END DO 341 342 CALL SUR_DETERMINE(jj+1) 300 ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,3) 301 CALL SUR_DETERMINE 343 302 344 303 ! Fill output array 345 304 DO jh = 1, nb_ana 346 ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1) 347 ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2) 348 END DO 349 350 END DO 351 END DO 352 353 DO jj = 1, jpj 354 DO ji = 1, jpi 355 DO jh = 1, nb_ana 356 X1=ana_amp(ji,jj,jh,1) 357 X2=-ana_amp(ji,jj,jh,2) 358 out_v(ji,jj, jh)=X1 * ssvmask(ji,jj) 359 out_v(ji,jj,nb_ana+jh)=X2 * ssvmask(ji,jj) 360 END DO 361 END DO 362 END DO 305 out_v(ji,jj, jh) = ztmp7((jh-1)*2+1) * ssvmask(ji,jj) 306 out_v(ji,jj,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssvmask(ji,jj) 307 END DO 308 309 END DO 310 END DO 311 ! 312 ! clem: we could avoid this call if all the loops were from 1:jpi and 1:jpj 313 ! but I think this is the most efficient 314 CALL lbc_lnk_multi( 'dia_harm_end', out_eta, 'T', 1., out_u, 'U', -1. , out_v, 'V', -1. ) 363 315 ! 364 316 CALL dia_wri_harm ! Write results in files 317 ! 318 DEALLOCATE( out_eta, out_u, out_v ) 365 319 ! 366 320 END SUBROUTINE dia_harm_end … … 373 327 !! ** Purpose : Write tidal harmonic analysis results in a netcdf file 374 328 !!-------------------------------------------------------------------- 375 CHARACTER(LEN=lc) :: cltext376 CHARACTER(LEN=lc) :: &377 cdfile_name_T , & ! name of the file created (T-points)378 cdfile_name_U , & ! name of the file created (U-points)379 cdfile_name_V ! name of the file created (V-points)380 329 INTEGER :: jh 381 330 !!---------------------------------------------------------------------- … … 383 332 IF(lwp) WRITE(numout,*) ' ' 384 333 IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results' 385 IF(lwp) WRITE(numout,*) ' 334 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 386 335 387 336 ! A) Elevation 388 337 !///////////// 389 !390 338 DO jh = 1, nb_ana 391 339 CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) ) 392 CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:, nb_ana+jh) )340 CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,jh+nb_ana) ) 393 341 END DO 394 342 395 343 ! B) ubar 396 344 !///////// 397 !398 345 DO jh = 1, nb_ana 399 346 CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) ) 400 CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:, nb_ana+jh) )347 CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,jh+nb_ana) ) 401 348 END DO 402 349 403 350 ! C) vbar 404 351 !///////// 405 !406 352 DO jh = 1, nb_ana 407 353 CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh ) ) … … 412 358 413 359 414 SUBROUTINE SUR_DETERMINE(init) 360 SUBROUTINE SUR_DETERMINE_INIT 361 !!--------------------------------------------------------------------------------- 362 !! *** ROUTINE SUR_DETERMINE_INIT *** 363 !! 364 !!--------------------------------------------------------------------------------- 365 INTEGER :: ji_sd, jj_sd, ji1_sd, ji2_sd, jh1_sd, jh2_sd 366 INTEGER :: ipivot 367 REAL(wp) :: zval1, zval2, zcol1, zcol2 368 INTEGER , DIMENSION(jpincomax) :: ipos2 369 !!--------------------------------------------------------------------------------- 370 ! 371 ! 372 ztmp3(:,:) = 0._wp 373 ! 374 DO jh1_sd = 1, nsparse 375 DO jh2_sd = 1, nsparse 376 IF( nisparse(jh2_sd) == nisparse(jh1_sd) ) THEN 377 ztmp3(njsparse(jh1_sd),njsparse(jh2_sd)) = ztmp3(njsparse(jh1_sd),njsparse(jh2_sd)) & 378 & + valuesparse(jh1_sd)*valuesparse(jh2_sd) 379 ENDIF 380 END DO 381 END DO 382 ! 383 DO jj_sd = 1, ninco 384 ipos1(jj_sd) = jj_sd 385 ipos2(jj_sd) = jj_sd 386 END DO 387 ! 388 DO ji_sd = 1, ninco 389 ! 390 !find greatest non-zero pivot: 391 zval1 = ABS(ztmp3(ji_sd,ji_sd)) 392 ! 393 ipivot = ji_sd 394 DO jj_sd = ji_sd, ninco 395 zval2 = ABS(ztmp3(ji_sd,jj_sd)) 396 IF( zval2 >= zval1 )THEN 397 ipivot = jj_sd 398 zval1 = zval2 399 ENDIF 400 END DO 401 ! 402 DO ji1_sd = 1, ninco 403 zcol1 = ztmp3(ji1_sd,ji_sd) 404 zcol2 = ztmp3(ji1_sd,ipivot) 405 ztmp3(ji1_sd,ji_sd) = zcol2 406 ztmp3(ji1_sd,ipivot) = zcol1 407 END DO 408 ! 409 ipos2(ji_sd) = ipos1(ipivot) 410 ipos2(ipivot) = ipos1(ji_sd) 411 ipos1(ji_sd) = ipos2(ji_sd) 412 ipos1(ipivot) = ipos2(ipivot) 413 z1_pivot(ji_sd) = 1._wp / ztmp3(ji_sd,ji_sd) 414 DO jj_sd = 1, ninco 415 ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) * z1_pivot(ji_sd) 416 END DO 417 ! 418 DO ji2_sd = ji_sd+1, ninco 419 zpilier(ji2_sd,ji_sd) = ztmp3(ji2_sd,ji_sd) 420 DO jj_sd=1,ninco 421 ztmp3(ji2_sd,jj_sd) = ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd) 422 END DO 423 END DO 424 ! 425 END DO 426 ! 427 z1_tmp3 = 1._wp / ztmp3(ninco,ninco) 428 ! 429 END SUBROUTINE SUR_DETERMINE_INIT 430 431 432 SUBROUTINE SUR_DETERMINE 415 433 !!--------------------------------------------------------------------------------- 416 434 !! *** ROUTINE SUR_DETERMINE *** 417 435 !! 418 !! 419 !! 420 !!--------------------------------------------------------------------------------- 421 INTEGER, INTENT(in) :: init 422 ! 423 INTEGER :: ji_sd, jj_sd, ji1_sd, ji2_sd, jh1_sd, jh2_sd 424 REAL(wp) :: zval1, zval2, zx1 425 REAL(wp), DIMENSION(jpincomax) :: ztmpx, zcol1, zcol2 426 INTEGER , DIMENSION(jpincomax) :: ipos2, ipivot 427 !--------------------------------------------------------------------------------- 436 !!--------------------------------------------------------------------------------- 437 INTEGER :: ji_sd, jj_sd, ji1_sd, ji2_sd 438 REAL(wp) :: zx1 439 REAL(wp), DIMENSION(jpincomax) :: ztmpx 440 !!--------------------------------------------------------------------------------- 428 441 ! 429 IF( init == 1 ) THEN430 IF( nsparse > jpdimsparse ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')431 IF( ninco > jpincomax ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')432 !433 ztmp3(:,:) = 0._wp434 !435 DO jh1_sd = 1, nsparse436 DO jh2_sd = 1, nsparse437 nisparse(jh2_sd) = nisparse(jh2_sd)438 njsparse(jh2_sd) = njsparse(jh2_sd)439 IF( nisparse(jh2_sd) == nisparse(jh1_sd) ) THEN440 ztmp3(njsparse(jh1_sd),njsparse(jh2_sd)) = ztmp3(njsparse(jh1_sd),njsparse(jh2_sd)) &441 & + valuesparse(jh1_sd)*valuesparse(jh2_sd)442 ENDIF443 END DO444 END DO445 !446 DO jj_sd = 1 ,ninco447 ipos1(jj_sd) = jj_sd448 ipos2(jj_sd) = jj_sd449 END DO450 !451 DO ji_sd = 1 , ninco452 !453 !find greatest non-zero pivot:454 zval1 = ABS(ztmp3(ji_sd,ji_sd))455 !456 ipivot(ji_sd) = ji_sd457 DO jj_sd = ji_sd, ninco458 zval2 = ABS(ztmp3(ji_sd,jj_sd))459 IF( zval2 >= zval1 )THEN460 ipivot(ji_sd) = jj_sd461 zval1 = zval2462 ENDIF463 END DO464 !465 DO ji1_sd = 1, ninco466 zcol1(ji1_sd) = ztmp3(ji1_sd,ji_sd)467 zcol2(ji1_sd) = ztmp3(ji1_sd,ipivot(ji_sd))468 ztmp3(ji1_sd,ji_sd) = zcol2(ji1_sd)469 ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd)470 END DO471 !472 ipos2(ji_sd) = ipos1(ipivot(ji_sd))473 ipos2(ipivot(ji_sd)) = ipos1(ji_sd)474 ipos1(ji_sd) = ipos2(ji_sd)475 ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd))476 zpivot(ji_sd) = ztmp3(ji_sd,ji_sd)477 DO jj_sd = 1, ninco478 ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd)479 END DO480 !481 DO ji2_sd = ji_sd+1, ninco482 zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd)483 DO jj_sd=1,ninco484 ztmp3(ji2_sd,jj_sd)= ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)485 END DO486 END DO487 !488 END DO489 !490 ENDIF ! End init==1491 492 442 DO ji_sd = 1, ninco 493 ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd)443 ztmp4(ji_sd) = ztmp4(ji_sd) * z1_pivot(ji_sd) 494 444 DO ji2_sd = ji_sd+1, ninco 495 445 ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd) … … 498 448 499 449 !system solving: 500 ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco) 501 ji_sd = ninco 450 ztmpx(ninco) = ztmp4(ninco) * z1_tmp3 502 451 DO ji_sd = ninco-1, 1, -1 503 452 zx1 = 0._wp … … 505 454 zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd) 506 455 END DO 507 ztmpx(ji_sd) = ztmp4(ji_sd) -zx1508 END DO 509 510 DO jj_sd = 1, ninco511 ztmp7(ipos1(jj_sd)) =ztmpx(jj_sd)456 ztmpx(ji_sd) = ztmp4(ji_sd) - zx1 457 END DO 458 459 DO jj_sd = 1, ninco 460 ztmp7(ipos1(jj_sd)) = ztmpx(jj_sd) 512 461 END DO 513 462 ! 514 463 END SUBROUTINE SUR_DETERMINE 515 464 465 516 466 !!====================================================================== 517 467 END MODULE diaharm
Note: See TracChangeset
for help on using the changeset viewer.
