Changeset 4683
- Timestamp:
- 2014-06-21T13:23:16+02:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
r4624 r4683 18 18 USE daymod 19 19 USE tide_mod 20 ! 20 21 USE in_out_manager ! I/O units 21 22 USE iom ! I/0 library … … 34 35 INTEGER, PARAMETER :: jpdimsparse = jpincomax*300*24 35 36 36 ! !!!namelist variables37 ! !!** namelist variables ** 37 38 INTEGER :: nit000_han ! First time step used for harmonic analysis 38 39 INTEGER :: nitend_han ! Last time step used for harmonic analysis 39 40 INTEGER :: nstep_han ! Time step frequency for harmonic analysis 40 INTEGER :: nb_ana 41 INTEGER :: nb_ana ! Number of harmonics to analyse 41 42 42 43 INTEGER , ALLOCATABLE, DIMENSION(:) :: name … … 119 120 ENDIF 120 121 END DO 121 END DO122 END DO 122 123 ! 123 124 IF(lwp) THEN … … 158 159 ! ---------------------------- 159 160 ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) ) 160 ana_temp(:,:,:,:) = 0. e0161 ana_temp(:,:,:,:) = 0._wp 161 162 162 163 END SUBROUTINE dia_harm_init … … 179 180 IF( nn_timing == 1 ) CALL timing_start('dia_harm') 180 181 181 IF ( kt == nit000 ) CALL dia_harm_init 182 183 IF ( ((kt.GE.nit000_han).AND.(kt.LE.nitend_han)).AND. & 184 (MOD(kt,nstep_han).EQ.0) ) THEN 185 186 ztime = (kt-nit000+1)*rdt 182 IF( kt == nit000 ) CALL dia_harm_init 183 184 IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN 185 186 ztime = (kt-nit000+1) * rdt 187 187 188 nhc = 0189 DO jh = 1,nb_ana190 DO jc = 1,2191 nhc = nhc+1192 ztemp =( MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) &193 +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)))194 195 DO jj = 1,jpj196 DO ji = 1,jpi197 ! Elevation198 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj) *tmask(ji,jj,1)188 nhc = 0 189 DO jh = 1, nb_ana 190 DO jc = 1, 2 191 nhc = nhc+1 192 ztemp =( MOD(jc,2) * ft(jh) *COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) & 193 & +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh))) 194 195 DO jj = 1, jpj 196 DO ji = 1, jpi 197 ! Elevation 198 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj) *tmask(ji,jj,1) 199 199 #if defined key_dynspg_ts 200 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1)201 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1)202 #endif 203 END DO204 END DO205 206 END DO207 END DO208 200 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1) 201 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1) 202 #endif 203 END DO 204 END DO 205 ! 206 END DO 207 END DO 208 ! 209 209 END IF 210 210 … … 249 249 keq = keq + 1 250 250 kun = 0 251 DO jh = 1, nb_ana252 DO jc = 1, 2251 DO jh = 1, nb_ana 252 DO jc = 1, 2 253 253 kun = kun + 1 254 254 ksp = ksp + 1 … … 296 296 out_eta(ji,jj,jh ) = X1 * tmask(ji,jj,1) 297 297 out_eta(ji,jj,jh+nb_ana) = X2 * tmask(ji,jj,1) 298 END DO299 END DO300 END DO298 END DO 299 END DO 300 END DO 301 301 302 302 ! ubar: … … 309 309 kun = kun + 1 310 310 ztmp4(kun)=ana_temp(ji,jj,kun,2) 311 END DO312 END DO311 END DO 312 END DO 313 313 314 314 CALL SUR_DETERMINE(jj+1) … … 316 316 ! Fill output array 317 317 DO jh = 1, nb_ana 318 ana_amp(ji,jj,jh,1) =ztmp7((jh-1)*2+1)319 ana_amp(ji,jj,jh,2) =ztmp7((jh-1)*2+2)318 ana_amp(ji,jj,jh,1) = ztmp7((jh-1)*2+1) 319 ana_amp(ji,jj,jh,2) = ztmp7((jh-1)*2+2) 320 320 END DO 321 321 … … 326 326 DO ji = 1, jpi 327 327 DO jh = 1, nb_ana 328 X1 =ana_amp(ji,jj,jh,1)329 X2 =-ana_amp(ji,jj,jh,2)330 out_u(ji,jj,jh ) = X1 * umask(ji,jj,1)331 out_u 332 END DO333 END DO334 END DO328 X1 = ana_amp(ji,jj,jh,1) 329 X2 =-ana_amp(ji,jj,jh,2) 330 out_u(ji,jj,jh ) = X1 * umask(ji,jj,1) 331 out_u(ji,jj,nb_ana+jh) = X2 * umask(ji,jj,1) 332 END DO 333 END DO 334 END DO 335 335 336 336 ! vbar: … … 343 343 kun = kun + 1 344 344 ztmp4(kun)=ana_temp(ji,jj,kun,3) 345 END DO346 END DO345 END DO 346 END DO 347 347 348 348 CALL SUR_DETERMINE(jj+1) … … 364 364 out_v(ji,jj,jh)=X1 * vmask(ji,jj,1) 365 365 out_v(ji,jj,nb_ana+jh)=X2 * vmask(ji,jj,1) 366 END DO367 END DO368 END DO366 END DO 367 END DO 368 END DO 369 369 370 370 CALL dia_wri_harm ! Write results in files … … 437 437 #else 438 438 DO jh = 1, nb_ana 439 CALL iom_put( TRIM(tname(jh))//'x_v', out_ u(:,:,jh ) )440 CALL iom_put( TRIM(tname(jh))//'y_v', out_ u(:,:,jh+nb_ana) )441 END DO 442 #endif 443 439 CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh ) ) 440 CALL iom_put( TRIM(tname(jh))//'y_v', out_v(:,:,jh+nb_ana) ) 441 END DO 442 #endif 443 ! 444 444 END SUBROUTINE dia_wri_harm 445 445 446 446 447 447 SUBROUTINE SUR_DETERMINE(init) 448 !!---------------------------------------------------------------------------------449 !! *** ROUTINE SUR_DETERMINE ***450 !!451 !!452 !!453 !!---------------------------------------------------------------------------------454 INTEGER, INTENT(in) :: init455 !456 INTEGER :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd457 REAL(wp) :: zval1, zval2, zx1458 REAL(wp), POINTER, DIMENSION(:) :: ztmpx, zcol1, zcol2459 INTEGER , POINTER, DIMENSION(:) :: ipos2, ipivot460 !---------------------------------------------------------------------------------461 CALL wrk_alloc( jpincomax , ztmpx , zcol1 , zcol2 )462 CALL wrk_alloc( jpincomax , ipos2 , ipivot )448 !!--------------------------------------------------------------------------------- 449 !! *** ROUTINE SUR_DETERMINE *** 450 !! 451 !! 452 !! 453 !!--------------------------------------------------------------------------------- 454 INTEGER, INTENT(in) :: init 455 ! 456 INTEGER :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd 457 REAL(wp) :: zval1, zval2, zx1 458 REAL(wp), POINTER, DIMENSION(:) :: ztmpx, zcol1, zcol2 459 INTEGER , POINTER, DIMENSION(:) :: ipos2, ipivot 460 !--------------------------------------------------------------------------------- 461 CALL wrk_alloc( jpincomax , ztmpx , zcol1 , zcol2 ) 462 CALL wrk_alloc( jpincomax , ipos2 , ipivot ) 463 463 464 IF( init == 1 ) THEN 465 IF( nsparse > jpdimsparse ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse') 466 IF( ninco > jpincomax ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax') 467 ! 468 ztmp3(:,:) = 0._wp 469 ! 470 DO jk1_sd = 1, nsparse 471 DO jk2_sd = 1, nsparse 472 nisparse(jk2_sd) = nisparse(jk2_sd) 473 njsparse(jk2_sd) = njsparse(jk2_sd) 474 IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN 475 ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) & 476 + valuesparse(jk1_sd)*valuesparse(jk2_sd) 477 ENDIF 478 END DO 479 END DO 480 481 DO jj_sd = 1 ,ninco 482 ipos1(jj_sd) = jj_sd 483 ipos2(jj_sd) = jj_sd 484 ENDDO 485 486 DO ji_sd = 1 , ninco 487 488 !find greatest non-zero pivot: 489 zval1 = ABS(ztmp3(ji_sd,ji_sd)) 490 491 ipivot(ji_sd) = ji_sd 492 DO jj_sd = ji_sd, ninco 493 zval2 = ABS(ztmp3(ji_sd,jj_sd)) 494 IF( zval2.GE.zval1 )THEN 495 ipivot(ji_sd) = jj_sd 496 zval1 = zval2 497 ENDIF 498 ENDDO 499 500 DO ji1_sd = 1, ninco 501 zcol1(ji1_sd) = ztmp3(ji1_sd,ji_sd) 502 zcol2(ji1_sd) = ztmp3(ji1_sd,ipivot(ji_sd)) 503 ztmp3(ji1_sd,ji_sd) = zcol2(ji1_sd) 504 ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd) 505 ENDDO 506 507 ipos2(ji_sd) = ipos1(ipivot(ji_sd)) 508 ipos2(ipivot(ji_sd)) = ipos1(ji_sd) 509 ipos1(ji_sd) = ipos2(ji_sd) 510 ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd)) 511 zpivot(ji_sd) = ztmp3(ji_sd,ji_sd) 512 DO jj_sd = 1, ninco 513 ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd) 514 ENDDO 515 464 IF( init == 1 ) THEN 465 IF( nsparse > jpdimsparse ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse') 466 IF( ninco > jpincomax ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax') 467 ! 468 ztmp3(:,:) = 0._wp 469 ! 470 DO jk1_sd = 1, nsparse 471 DO jk2_sd = 1, nsparse 472 nisparse(jk2_sd) = nisparse(jk2_sd) 473 njsparse(jk2_sd) = njsparse(jk2_sd) 474 IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN 475 ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) & 476 & + valuesparse(jk1_sd)*valuesparse(jk2_sd) 477 ENDIF 478 END DO 479 END DO 480 ! 481 DO jj_sd = 1 ,ninco 482 ipos1(jj_sd) = jj_sd 483 ipos2(jj_sd) = jj_sd 484 END DO 485 ! 486 DO ji_sd = 1 , ninco 487 ! 488 !find greatest non-zero pivot: 489 zval1 = ABS(ztmp3(ji_sd,ji_sd)) 490 ! 491 ipivot(ji_sd) = ji_sd 492 DO jj_sd = ji_sd, ninco 493 zval2 = ABS(ztmp3(ji_sd,jj_sd)) 494 IF( zval2.GE.zval1 )THEN 495 ipivot(ji_sd) = jj_sd 496 zval1 = zval2 497 ENDIF 498 END DO 499 ! 500 DO ji1_sd = 1, ninco 501 zcol1(ji1_sd) = ztmp3(ji1_sd,ji_sd) 502 zcol2(ji1_sd) = ztmp3(ji1_sd,ipivot(ji_sd)) 503 ztmp3(ji1_sd,ji_sd) = zcol2(ji1_sd) 504 ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd) 505 END DO 506 ! 507 ipos2(ji_sd) = ipos1(ipivot(ji_sd)) 508 ipos2(ipivot(ji_sd)) = ipos1(ji_sd) 509 ipos1(ji_sd) = ipos2(ji_sd) 510 ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd)) 511 zpivot(ji_sd) = ztmp3(ji_sd,ji_sd) 512 DO jj_sd = 1, ninco 513 ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd) 514 END DO 515 ! 516 DO ji2_sd = ji_sd+1, ninco 517 zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd) 518 DO jj_sd=1,ninco 519 ztmp3(ji2_sd,jj_sd)= ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd) 520 END DO 521 END DO 522 ! 523 END DO 524 ! 525 ENDIF ! End init==1 526 527 DO ji_sd = 1, ninco 528 ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd) 516 529 DO ji2_sd = ji_sd+1, ninco 517 zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd) 518 DO jj_sd=1,ninco 519 ztmp3(ji2_sd,jj_sd)= ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd) 520 ENDDO 521 ENDDO 522 523 ENDDO 524 525 ENDIF ! End init==1 526 527 DO ji_sd = 1, ninco 528 ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd) 529 DO ji2_sd = ji_sd+1, ninco 530 ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd) 531 ENDDO 532 ENDDO 533 534 !system solving: 535 ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco) 536 ji_sd = ninco 537 DO ji_sd = ninco-1, 1, -1 538 zx1=0. 539 DO jj_sd = ji_sd+1, ninco 540 zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd) 541 ENDDO 542 ztmpx(ji_sd) = ztmp4(ji_sd)-zx1 543 ENDDO 544 545 DO jj_sd =1, ninco 546 ztmp7(ipos1(jj_sd))=ztmpx(jj_sd) 547 ENDDO 548 549 CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 ) 550 CALL wrk_dealloc( jpincomax , ipos2 , ipivot ) 551 552 END SUBROUTINE SUR_DETERMINE 530 ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd) 531 END DO 532 END DO 533 534 !system solving: 535 ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco) 536 ji_sd = ninco 537 DO ji_sd = ninco-1, 1, -1 538 zx1 = 0._wp 539 DO jj_sd = ji_sd+1, ninco 540 zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd) 541 END DO 542 ztmpx(ji_sd) = ztmp4(ji_sd)-zx1 543 END DO 544 545 DO jj_sd =1, ninco 546 ztmp7(ipos1(jj_sd))=ztmpx(jj_sd) 547 END DO 548 549 CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 ) 550 CALL wrk_dealloc( jpincomax , ipos2 , ipivot ) 551 ! 552 END SUBROUTINE SUR_DETERMINE 553 553 554 554 #else
Note: See TracChangeset
for help on using the changeset viewer.