Changeset 4292 for branches/2013/dev_MERGE_2013
- Timestamp:
- 2013-11-20T17:28:04+01:00 (7 years ago)
- Location:
- branches/2013/dev_MERGE_2013/NEMOGCM
- Files:
-
- 1 deleted
- 88 edited
- 2 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/SHARED/1_namelist_ref
r4230 r4292 507 507 !----------------------------------------------------------------------- 508 508 ln_tide_pot = .true. ! use tidal potential forcing 509 nb_harmo = 11 ! number of constituents used 509 510 clname(1) = 'M2' ! name of constituent 510 511 clname(2) = 'S2' … … 700 701 ln_dynadv_cen2= .false. ! flux form - 2nd order centered scheme 701 702 ln_dynadv_ubs = .false. ! flux form - 3rd order UBS scheme 703 / 704 !----------------------------------------------------------------------- 705 &nam_vvl ! vertical coordinate options 706 !----------------------------------------------------------------------- 707 ln_vvl_zstar = .true. ! zstar vertical coordinate 708 ln_vvl_ztilde = .false. ! ztilde vertical coordinate: only high frequency variations 709 ln_vvl_layer = .false. ! full layer vertical coordinate 710 ln_vvl_ztilde_as_zstar = .false. ! ztilde vertical coordinate emulating zstar 711 rn_ahe3 = 0.0e0 ! thickness diffusion coefficient 712 rn_rst_e3t = 30.e0 ! ztilde to zstar restoration timescale [days] 713 rn_lf_cutoff = 5.0e0 ! cutoff frequency for low-pass filter [days] 714 rn_zdef_max = 0.9e0 ! maximum fractional e3t deformation 715 ln_vvl_dbg = .true. ! debug prints (T/F) 702 716 / 703 717 !----------------------------------------------------------------------- -
branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/SHARED/namelist_ref
r4245 r4292 513 513 !----------------------------------------------------------------------- 514 514 ln_tide_pot = .true. ! use tidal potential forcing 515 nb_harmo = 11 ! number of constituents used 515 516 clname(1) = 'M2' ! name of constituent 516 517 clname(2) = 'S2' … … 712 713 ln_dynadv_cen2= .false. ! flux form - 2nd order centered scheme 713 714 ln_dynadv_ubs = .false. ! flux form - 3rd order UBS scheme 715 / 716 !----------------------------------------------------------------------- 717 &nam_vvl ! vertical coordinate options 718 !----------------------------------------------------------------------- 719 ln_vvl_zstar = .true. ! zstar vertical coordinate 720 ln_vvl_ztilde = .false. ! ztilde vertical coordinate: only high frequency variations 721 ln_vvl_layer = .false. ! full layer vertical coordinate 722 ln_vvl_ztilde_as_zstar = .false. ! ztilde vertical coordinate emulating zstar 723 rn_ahe3 = 0.0e0 ! thickness diffusion coefficient 724 rn_rst_e3t = 30.e0 ! ztilde to zstar restoration timescale [days] 725 rn_lf_cutoff = 5.0e0 ! cutoff frequency for low-pass filter [days] 726 rn_zdef_max = 0.9e0 ! maximum fractional e3t deformation 727 ln_vvl_dbg = .true. ! debug prints (T/F) 714 728 / 715 729 !----------------------------------------------------------------------- -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90
r4148 r4292 24 24 USE phycst ! physical constants 25 25 USE dom_oce ! ocean domain 26 USE domvvl ! ocean vertical scale factors 26 27 USE dom_ice_2 ! LIM-2: ice domain 27 28 USE ice_2 ! LIM-2: ice variables … … 59 60 !! * Substitutions 60 61 # include "vectopt_loop_substitute.h90" 62 # include "domzgr_substitute.h90" 61 63 !!---------------------------------------------------------------------- 62 64 !! NEMO/LIM2 4.0 , UCL - NEMO Consortium (2011) … … 446 448 !!------------------------------------------------------------------- 447 449 ! 450 INTEGER :: jk ! local integer 451 ! 448 452 IF(lwp) WRITE(numout,*) 449 453 IF(lwp) WRITE(numout,*) 'lim_sbc_init_2 : LIM-2 sea-ice - surface boundary condition' … … 472 476 snwice_mass (:,:) = 0.e0 ! no mass exchanges 473 477 snwice_mass_b(:,:) = 0.e0 ! no mass exchanges 478 snwice_fmass (:,:) = 0.e0 ! no mass exchanges 474 479 ENDIF 475 480 IF( nn_ice_embd == 2 .AND. & ! full embedment (case 2) & no restart : … … 477 482 sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 478 483 sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 484 do jk = 1,jpkm1 ! adjust initial vertical scale factors 485 fse3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 486 fse3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 487 end do 488 fse3t_a(:,:,:) = fse3t_b(:,:,:) 489 ! Reconstruction of all vertical scale factors at now and before time steps 490 ! ============================================================================= 491 ! Horizontal scale factor interpolations 492 ! -------------------------------------- 493 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 494 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 495 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 496 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 497 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' ) 498 ! Vertical scale factor interpolations 499 ! ------------------------------------ 500 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 501 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 502 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 503 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 504 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 505 ! t- and w- points depth 506 ! ---------------------- 507 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 508 fsdepw_n(:,:,1) = 0.0_wp 509 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 510 DO jk = 2, jpk 511 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 512 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 513 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 514 END DO 479 515 ENDIF 480 516 ! -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90
r4147 r4292 237 237 238 238 ! energy needed to bring ocean surface layer until its freezing 239 qcmif (ji,jj) = rau0 * rcp * fse3t_m(ji,jj ,1) * ( tfu(ji,jj) - sst_m(ji,jj) - rt0 ) * ( 1 - zinda )239 qcmif (ji,jj) = rau0 * rcp * fse3t_m(ji,jj) * ( tfu(ji,jj) - sst_m(ji,jj) - rt0 ) * ( 1 - zinda ) 240 240 241 241 ! calculate oceanic heat flux. -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90
r4147 r4292 40 40 INTEGER :: e1u_id, e2v_id, sshn_id, gcb_id 41 41 INTEGER :: trn_id, trb_id, tra_id 42 INTEGER :: unb_id, vnb_id 42 43 43 44 !!---------------------------------------------------------------------- -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
r3294 r4292 27 27 USE agrif_opa_sponge 28 28 USE lib_mpp 29 USE wrk_nemo 29 USE wrk_nemo 30 USE dynspg_oce 30 31 31 32 IMPLICIT NONE 32 33 PRIVATE 34 35 ! Barotropic arrays used to store open boundary data during 36 ! time-splitting loop: 37 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: ubdy_w, vbdy_w, hbdy_w 38 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: ubdy_e, vbdy_e, hbdy_e 39 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: ubdy_n, vbdy_n, hbdy_n 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: ubdy_s, vbdy_s, hbdy_s 33 41 34 PUBLIC Agrif_tra, Agrif_dyn, Agrif_ssh, interpu, interpv 42 PUBLIC Agrif_tra, Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ssh_ts 43 PUBLIC interpu, interpv, interpunb, interpvnb, interpsshn 35 44 36 45 # include "domzgr_substitute.h90" … … 169 178 REAL(wp) :: timeref 170 179 REAL(wp) :: z2dt, znugdt 171 REAL(wp) :: zrhox, rhoy180 REAL(wp) :: zrhox, zrhoy 172 181 REAL(wp), POINTER, DIMENSION(:,:,:) :: zua, zva 173 182 REAL(wp), POINTER, DIMENSION(:,:) :: spgv1, spgu1, zua2d, zva2d … … 180 189 181 190 zrhox = Agrif_Rhox() 182 rhoy = Agrif_Rhoy()191 zrhoy = Agrif_Rhoy() 183 192 184 193 timeref = 1. … … 201 210 zva2d = 0. 202 211 212 #if defined key_dynspg_flt 203 213 Agrif_SpecialValue=0. 204 214 Agrif_UseSpecialValue = ln_spc_dyn 205 215 CALL Agrif_Bc_variable(zua2d,e1u_id,calledweight=1.,procname=interpu2d) 206 216 CALL Agrif_Bc_variable(zva2d,e2v_id,calledweight=1.,procname=interpv2d) 217 #endif 207 218 Agrif_UseSpecialValue = .FALSE. 208 219 … … 210 221 IF((nbondi == -1).OR.(nbondi == 2)) THEN 211 222 223 #if defined key_dynspg_flt 212 224 DO jj=1,jpj 213 laplacu(2,jj) = timeref * (zua2d(2,jj)/(rhoy*e2u(2,jj)))*umask(2,jj,1) 214 END DO 215 216 DO jk=1,jpkm1 217 DO jj=1,jpj 218 ua(1:2,jj,jk) = (zua(1:2,jj,jk)/(rhoy*e2u(1:2,jj))) 225 laplacu(2,jj) = timeref * (zua2d(2,jj)/(zrhoy*e2u(2,jj)))*umask(2,jj,1) 226 END DO 227 #endif 228 229 DO jk=1,jpkm1 230 DO jj=1,jpj 231 ua(1:2,jj,jk) = (zua(1:2,jj,jk)/(zrhoy*e2u(1:2,jj))) 219 232 ua(1:2,jj,jk) = ua(1:2,jj,jk) / fse3u(1:2,jj,jk) 220 233 END DO 221 234 END DO 222 235 236 #if defined key_dynspg_flt 223 237 DO jk=1,jpkm1 224 238 DO jj=1,jpj … … 240 254 ENDIF 241 255 END DO 256 #else 257 spgu(2,:) = ua_b(2,:) 258 #endif 242 259 243 260 DO jk=1,jpkm1 … … 278 295 279 296 IF((nbondi == 1).OR.(nbondi == 2)) THEN 280 297 #if defined key_dynspg_flt 281 298 DO jj=1,jpj 282 laplacu(nlci-2,jj) = timeref * (zua2d(nlci-2,jj)/(rhoy*e2u(nlci-2,jj))) 283 END DO 284 285 DO jk=1,jpkm1 286 DO jj=1,jpj 287 ua(nlci-2:nlci-1,jj,jk) = (zua(nlci-2:nlci-1,jj,jk)/(rhoy*e2u(nlci-2:nlci-1,jj))) 299 laplacu(nlci-2,jj) = timeref * (zua2d(nlci-2,jj)/(zrhoy*e2u(nlci-2,jj))) 300 END DO 301 #endif 302 303 DO jk=1,jpkm1 304 DO jj=1,jpj 305 ua(nlci-2:nlci-1,jj,jk) = (zua(nlci-2:nlci-1,jj,jk)/(zrhoy*e2u(nlci-2:nlci-1,jj))) 288 306 289 307 ua(nlci-2:nlci-1,jj,jk) = ua(nlci-2:nlci-1,jj,jk) / fse3u(nlci-2:nlci-1,jj,jk) … … 292 310 END DO 293 311 312 #if defined key_dynspg_flt 294 313 DO jk=1,jpkm1 295 314 DO jj=1,jpj … … 312 331 ENDIF 313 332 END DO 333 #else 334 spgu(nlci-2,:) = ua_b(nlci-2,:) 335 #endif 314 336 315 337 DO jk=1,jpkm1 … … 353 375 IF((nbondj == -1).OR.(nbondj == 2)) THEN 354 376 377 #if defined key_dynspg_flt 355 378 DO ji=1,jpi 356 379 laplacv(ji,2) = timeref * (zva2d(ji,2)/(zrhox*e1v(ji,2))) 357 380 END DO 381 #endif 358 382 359 383 DO jk=1,jpkm1 … … 364 388 END DO 365 389 390 #if defined key_dynspg_flt 366 391 DO jk=1,jpkm1 367 392 DO ji=1,jpi … … 383 408 ENDIF 384 409 END DO 410 #else 411 spgv(:,2)=va_b(:,2) 412 #endif 385 413 386 414 DO jk=1,jpkm1 … … 413 441 DO jk=1,jpkm1 414 442 DO ji=1,jpi 415 ua(ji,2,jk) = (zua(ji,2,jk)/( rhoy*e2u(ji,2)))*umask(ji,2,jk)443 ua(ji,2,jk) = (zua(ji,2,jk)/(zrhoy*e2u(ji,2)))*umask(ji,2,jk) 416 444 ua(ji,2,jk) = ua(ji,2,jk) / fse3u(ji,2,jk) 417 445 END DO … … 422 450 IF((nbondj == 1).OR.(nbondj == 2)) THEN 423 451 452 #if defined key_dynspg_flt 424 453 DO ji=1,jpi 425 454 laplacv(ji,nlcj-2) = timeref * (zva2d(ji,nlcj-2)/(zrhox*e1v(ji,nlcj-2))) 426 455 END DO 456 #endif 427 457 428 458 DO jk=1,jpkm1 … … 433 463 END DO 434 464 465 #if defined key_dynspg_flt 435 466 DO jk=1,jpkm1 436 467 DO ji=1,jpi … … 438 469 END DO 439 470 END DO 440 441 471 442 472 spgv(:,nlcj-2)=0. … … 453 483 ENDIF 454 484 END DO 485 #else 486 spgv(:,nlcj-2)=va_b(:,nlcj-2) 487 #endif 455 488 456 489 DO jk=1,jpkm1 … … 483 516 DO jk=1,jpkm1 484 517 DO ji=1,jpi 485 ua(ji,nlcj-1,jk) = (zua(ji,nlcj-1,jk)/( rhoy*e2u(ji,nlcj-1)))*umask(ji,nlcj-1,jk)518 ua(ji,nlcj-1,jk) = (zua(ji,nlcj-1,jk)/(zrhoy*e2u(ji,nlcj-1)))*umask(ji,nlcj-1,jk) 486 519 ua(ji,nlcj-1,jk) = ua(ji,nlcj-1,jk) / fse3u(ji,nlcj-1,jk) 487 520 END DO … … 495 528 END SUBROUTINE Agrif_dyn 496 529 530 SUBROUTINE Agrif_dyn_ts( kt, jn ) 531 !!---------------------------------------------------------------------- 532 !! *** ROUTINE Agrif_dyn_ts *** 533 !!---------------------------------------------------------------------- 534 !! 535 INTEGER, INTENT(in) :: kt, jn 536 !! 537 INTEGER :: ji, jj 538 REAL(wp) :: zrhox, zrhoy 539 REAL(wp), POINTER, DIMENSION(:,:) :: spgv1, spgu1 540 REAL(wp), POINTER, DIMENSION(:,:) :: zunb, zvnb, zsshn 541 !!---------------------------------------------------------------------- 542 543 IF( Agrif_Root() ) RETURN 544 545 IF ((kt==nit000).AND.(jn==1)) THEN 546 ALLOCATE( ubdy_w(jpj), vbdy_w(jpj), hbdy_w(jpj)) 547 ALLOCATE( ubdy_e(jpj), vbdy_e(jpj), hbdy_e(jpj)) 548 ALLOCATE( ubdy_n(jpi), vbdy_n(jpi), hbdy_n(jpi)) 549 ALLOCATE( ubdy_s(jpi), vbdy_s(jpi), hbdy_s(jpi)) 550 ENDIF 551 552 IF (jn==1) THEN 553 ! Fill boundary arrays at each baroclinic step 554 ! with Parent grid barotropic fluxes and sea level 555 ! 556 CALL wrk_alloc( jpi, jpj, zunb, zvnb, zsshn ) 557 558 zrhox = Agrif_Rhox() 559 zrhoy = Agrif_Rhoy() 560 561 !alt Agrif_SpecialValue = 0.e0 562 !alt Agrif_UseSpecialValue = .TRUE. 563 !alt CALL Agrif_Bc_variable(zsshn, sshn_id, procname=interpsshn ) 564 !alt Agrif_UseSpecialValue = .FALSE. 565 566 Agrif_SpecialValue=0. 567 Agrif_UseSpecialValue = ln_spc_dyn 568 zunb(:,:) = 0._wp ; zvnb(:,:) = 0._wp 569 CALL Agrif_Bc_variable(zunb,unb_id,procname=interpunb) 570 CALL Agrif_Bc_variable(zvnb,vnb_id,procname=interpvnb) 571 Agrif_UseSpecialValue = .FALSE. 572 573 IF((nbondi == -1).OR.(nbondi == 2)) THEN 574 DO jj=1,jpj 575 ubdy_w(jj) = (zunb(2,jj)/(zrhoy*e2u(2,jj))) 576 vbdy_w(jj) = (zvnb(2,jj)/(zrhox*e1v(2,jj))) 577 hbdy_w(jj) = zsshn(2,jj) 578 END DO 579 ENDIF 580 581 IF((nbondi == 1).OR.(nbondi == 2)) THEN 582 DO jj=1,jpj 583 ubdy_e(jj) = zunb(nlci-2,jj)/(zrhoy*e2u(nlci-2,jj)) 584 vbdy_e(jj) = zvnb(nlci-1,jj)/(zrhox*e1v(nlci-1,jj)) 585 hbdy_e(jj) = zsshn(nlci-1,jj) 586 END DO 587 ENDIF 588 589 IF((nbondj == -1).OR.(nbondj == 2)) THEN 590 DO ji=1,jpi 591 ubdy_s(ji) = zunb(ji,2)/(zrhoy*e2u(ji,2)) 592 vbdy_s(ji) = zvnb(ji,2)/(zrhox*e1v(ji,2)) 593 hbdy_s(ji) = zsshn(ji,2) 594 END DO 595 ENDIF 596 597 IF((nbondj == 1).OR.(nbondj == 2)) THEN 598 DO ji=1,jpi 599 ubdy_n(ji) = zunb(ji,nlcj-1)/(zrhoy*e2u(ji,nlcj-1)) 600 vbdy_n(ji) = zvnb(ji,nlcj-2)/(zrhox*e1v(ji,nlcj-2)) 601 hbdy_n(ji) = zsshn(ji,nlcj-1) 602 END DO 603 ENDIF 604 605 CALL wrk_dealloc( jpi, jpj, zunb, zvnb, zsshn ) 606 ENDIF ! jn==1 607 608 ! Then update velocities at each barotropic time step 609 IF((nbondi == -1).OR.(nbondi == 2)) THEN 610 DO jj=1,jpj 611 va_e(2,jj) = vbdy_w(jj) * hvr_e(2,jj) 612 ! Specified fluxes: 613 ua_e(2,jj) = ubdy_w(jj) * hur_e(2,jj) 614 ! Characteristics method: 615 !alt ua_e(2,jj) = 0.5_wp * ( ubdy_w(jj) * hur_e(2,jj) + ua_e(3,jj) & 616 !alt & - sqrt(grav * hur_e(2,jj)) * (sshn_e(3,jj) - hbdy_w(jj)) ) 617 END DO 618 ENDIF 619 620 IF((nbondi == 1).OR.(nbondi == 2)) THEN 621 DO jj=1,jpj 622 va_e(nlci-1,jj) = vbdy_e(jj) * hvr_e(nlci-1,jj) 623 ! Specified fluxes: 624 ua_e(nlci-2,jj) = ubdy_e(jj) * hur_e(nlci-2,jj) 625 ! Characteristics method: 626 !alt ua_e(nlci-2,jj) = 0.5_wp * ( ubdy_e(jj) * hur_e(nlci-2,jj) + ua_e(nlci-3,jj) & 627 !alt & + sqrt(grav * hur_e(nlci-2,jj)) * (sshn_e(nlci-2,jj) - hbdy_e(jj)) ) 628 END DO 629 ENDIF 630 631 IF((nbondj == -1).OR.(nbondj == 2)) THEN 632 DO ji=1,jpi 633 ua_e(ji,2) = ubdy_s(ji) * hur_e(ji,2) 634 ! Specified fluxes: 635 va_e(ji,2) = vbdy_s(ji) * hvr_e(ji,2) 636 ! Characteristics method: 637 !alt va_e(ji,2) = 0.5_wp * ( vbdy_s(ji) * hvr_e(ji,2) + va_e(ji,3) & 638 !alt & - sqrt(grav * hvr_e(ji,2)) * (sshn_e(ji,3) - hbdy_s(ji)) ) 639 END DO 640 ENDIF 641 642 IF((nbondj == 1).OR.(nbondj == 2)) THEN 643 DO ji=1,jpi 644 ua_e(ji,nlcj-1) = ubdy_n(ji) * hur_e(ji,nlcj-1) 645 ! Specified fluxes: 646 va_e(ji,nlcj-2) = vbdy_n(ji) * hvr_e(ji,nlcj-2) 647 ! Characteristics method: 648 !alt va_e(ji,nlcj-2) = 0.5_wp * ( vbdy_n(ji) * hvr_e(ji,nlcj-2) + va_e(ji,nlcj-3) & 649 !alt & + sqrt(grav * hvr_e(ji,nlcj-2)) * (sshn_e(ji,nlcj-2) - hbdy_n(ji)) ) 650 END DO 651 ENDIF 652 ! 653 END SUBROUTINE Agrif_dyn_ts 497 654 498 655 SUBROUTINE Agrif_ssh( kt ) … … 518 675 519 676 IF((nbondj == -1).OR.(nbondj == 2)) THEN 520 ssha(:,2)=ssh n(:,3)521 sshn(:,2)=ssh b(:,3)677 ssha(:,2)=ssha(:,3) 678 sshn(:,2)=sshn(:,3) 522 679 ENDIF 523 680 524 681 IF((nbondj == 1).OR.(nbondj == 2)) THEN 525 682 ssha(:,nlcj-1)=ssha(:,nlcj-2) 526 ssh a(:,nlcj-1)=sshn(:,nlcj-2)683 sshn(:,nlcj-1)=sshn(:,nlcj-2) 527 684 ENDIF 528 685 529 686 END SUBROUTINE Agrif_ssh 530 687 688 SUBROUTINE Agrif_ssh_ts( kt ) 689 !!---------------------------------------------------------------------- 690 !! *** ROUTINE Agrif_ssh_ts *** 691 !!---------------------------------------------------------------------- 692 INTEGER, INTENT(in) :: kt 693 !! 694 !!---------------------------------------------------------------------- 695 696 IF((nbondi == -1).OR.(nbondi == 2)) THEN 697 ssha_e(2,:) = ssha_e(3,:) 698 ENDIF 699 700 IF((nbondi == 1).OR.(nbondi == 2)) THEN 701 ssha_e(nlci-1,:) = ssha_e(nlci-2,:) 702 ENDIF 703 704 IF((nbondj == -1).OR.(nbondj == 2)) THEN 705 ssha_e(:,2) = ssha_e(:,3) 706 ENDIF 707 708 IF((nbondj == 1).OR.(nbondj == 2)) THEN 709 ssha_e(:,nlcj-1) = ssha_e(:,nlcj-2) 710 ENDIF 711 712 END SUBROUTINE Agrif_ssh_ts 713 714 SUBROUTINE interpsshn(tabres,i1,i2,j1,j2) 715 !!---------------------------------------------------------------------- 716 !! *** ROUTINE interpsshn *** 717 !!---------------------------------------------------------------------- 718 INTEGER, INTENT(in) :: i1,i2,j1,j2 719 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 720 !! 721 INTEGER :: ji,jj 722 !!---------------------------------------------------------------------- 723 724 tabres(i1:i2,j1:j2) = sshn(i1:i2,j1:j2) 725 726 END SUBROUTINE interpsshn 531 727 532 728 SUBROUTINE interpu(tabres,i1,i2,j1,j2,k1,k2) … … 611 807 612 808 END SUBROUTINE interpv2d 809 810 SUBROUTINE interpunb(tabres,i1,i2,j1,j2) 811 !!---------------------------------------------------------------------- 812 !! *** ROUTINE interpunb *** 813 !!---------------------------------------------------------------------- 814 INTEGER, INTENT(in) :: i1,i2,j1,j2 815 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 816 !! 817 INTEGER :: ji,jj,jk 818 !!---------------------------------------------------------------------- 819 820 tabres(:,:) = 0.e0 821 DO jk=1,jpkm1 822 DO jj=j1,j2 823 DO ji=i1,i2 824 tabres(ji,jj) = tabres(ji,jj) + e2u(ji,jj) * un(ji,jj,jk) & 825 * umask(ji,jj,jk) * fse3u(ji,jj,jk) 826 END DO 827 END DO 828 END DO 829 830 END SUBROUTINE interpunb 831 832 SUBROUTINE interpvnb(tabres,i1,i2,j1,j2) 833 !!---------------------------------------------------------------------- 834 !! *** ROUTINE interpvnb *** 835 !!---------------------------------------------------------------------- 836 INTEGER, INTENT(in) :: i1,i2,j1,j2 837 REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 838 !! 839 INTEGER :: ji,jj,jk 840 !!---------------------------------------------------------------------- 841 842 tabres(:,:) = 0.e0 843 DO jk=1,jpkm1 844 DO jj=j1,j2 845 DO ji=i1,i2 846 tabres(ji,jj) = tabres(ji,jj) + e1v(ji,jj) * vn(ji,jj,jk) & 847 * vmask(ji,jj,jk) * fse3v(ji,jj,jk) 848 END DO 849 END DO 850 END DO 851 852 END SUBROUTINE interpvnb 613 853 614 854 #else -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OFF_SRC/domain.F90
r4248 r4292 295 295 !! vertical scale factors. 296 296 !! 297 !! ** Method : - reference 1D vertical coordinate (gdep._ 0, e3._0)297 !! ** Method : - reference 1D vertical coordinate (gdep._1d, e3._1d) 298 298 !! - read/set ocean depth and ocean levels (bathy, mbathy) 299 299 !! - vertical coordinate (gdep., e3.) depending on the -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OFF_SRC/domrea.F90
r3680 r4292 25 25 26 26 PUBLIC dom_rea ! routine called by inidom.F90 27 !! * Substitutions 28 # include "domzgr_substitute.h90" 27 29 !!---------------------------------------------------------------------- 28 30 !! NEMO/OFF 3.3 , NEMO Consortium (2010) … … 173 175 CALL iom_get( inum4, jpdom_unknown, 'esigw', esigw ) 174 176 175 CALL iom_get( inum4, jpdom_data, 'e3t', e3t) ! scale factors176 CALL iom_get( inum4, jpdom_data, 'e3u', e3u)177 CALL iom_get( inum4, jpdom_data, 'e3v', e3v)178 CALL iom_get( inum4, jpdom_data, 'e3w', e3w)179 180 CALL iom_get( inum4, jpdom_unknown, 'gdept_ 0', gdept_0) ! depth181 CALL iom_get( inum4, jpdom_unknown, 'gdepw_ 0', gdepw_0)177 CALL iom_get( inum4, jpdom_data, 'e3t', fse3t_n(:,:,:) ) ! scale factors 178 CALL iom_get( inum4, jpdom_data, 'e3u', fse3u_n(:,:,:) ) 179 CALL iom_get( inum4, jpdom_data, 'e3v', fse3v_n(:,:,:) ) 180 CALL iom_get( inum4, jpdom_data, 'e3w', fse3w_n(:,:,:) ) 181 182 CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth 183 CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d ) 182 184 ENDIF 183 185 184 186 185 187 IF( ln_zps ) THEN ! z-coordinate - partial steps 186 CALL iom_get( inum4, jpdom_unknown, 'gdept_ 0', gdept_0 )! reference depth187 CALL iom_get( inum4, jpdom_unknown, 'gdepw_ 0', gdepw_0)188 CALL iom_get( inum4, jpdom_unknown, 'e3t_ 0' , e3t_0) ! reference scale factors189 CALL iom_get( inum4, jpdom_unknown, 'e3w_ 0' , e3w_0)188 CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! reference depth 189 CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d ) 190 CALL iom_get( inum4, jpdom_unknown, 'e3t_1d' , e3t_1d ) ! reference scale factors 191 CALL iom_get( inum4, jpdom_unknown, 'e3w_1d' , e3w_1d ) 190 192 ! 191 193 IF( nmsh <= 6 ) THEN ! 3D vertical scale factors 192 CALL iom_get( inum4, jpdom_data, 'e3t', e3t)193 CALL iom_get( inum4, jpdom_data, 'e3u', e3u)194 CALL iom_get( inum4, jpdom_data, 'e3v', e3v)195 CALL iom_get( inum4, jpdom_data, 'e3w', e3w)194 CALL iom_get( inum4, jpdom_data, 'e3t', fse3t_n(:,:,:) ) 195 CALL iom_get( inum4, jpdom_data, 'e3u', fse3u_n(:,:,:) ) 196 CALL iom_get( inum4, jpdom_data, 'e3v', fse3v_n(:,:,:) ) 197 CALL iom_get( inum4, jpdom_data, 'e3w', fse3w_n(:,:,:) ) 196 198 ELSE ! 2D bottom scale factors 197 199 CALL iom_get( inum4, jpdom_data, 'e3t_ps', e3tp ) … … 199 201 ! ! deduces the 3D scale factors 200 202 DO jk = 1, jpk 201 e3t(:,:,jk) = e3t_0(jk)! set to the ref. factors202 e3u(:,:,jk) = e3t_0(jk)203 e3v(:,:,jk) = e3t_0(jk)204 e3w(:,:,jk) = e3w_0(jk)203 fse3t_n(:,:,jk) = e3t_1d(jk) ! set to the ref. factors 204 fse3u_n(:,:,jk) = e3t_1d(jk) 205 fse3v_n(:,:,jk) = e3t_1d(jk) 206 fse3w_n(:,:,jk) = e3w_1d(jk) 205 207 END DO 206 208 DO jj = 1,jpj ! adjust the deepest values 207 209 DO ji = 1,jpi 208 210 ik = mbkt(ji,jj) 209 e3t(ji,jj,ik) = e3tp(ji,jj) * tmask(ji,jj,1) + e3t_0(1) * ( 1._wp - tmask(ji,jj,1) )210 e3w(ji,jj,ik) = e3wp(ji,jj) * tmask(ji,jj,1) + e3w_0(1) * ( 1._wp - tmask(ji,jj,1) )211 fse3t_n(ji,jj,ik) = e3tp(ji,jj) * tmask(ji,jj,1) + e3t_1d(1) * ( 1._wp - tmask(ji,jj,1) ) 212 fse3w_n(ji,jj,ik) = e3wp(ji,jj) * tmask(ji,jj,1) + e3w_1d(1) * ( 1._wp - tmask(ji,jj,1) ) 211 213 END DO 212 214 END DO … … 214 216 DO jj = 1, jpjm1 215 217 DO ji = 1, jpim1 216 e3u(ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) )217 e3v(ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) )218 fse3u_n(ji,jj,jk) = MIN( fse3t_n(ji,jj,jk), fse3t_n(ji+1,jj,jk) ) 219 fse3v_n(ji,jj,jk) = MIN( fse3t_n(ji,jj,jk), fse3t_n(ji,jj+1,jk) ) 218 220 END DO 219 221 END DO 220 222 END DO 221 CALL lbc_lnk( e3u , 'U', 1._wp ) ; CALL lbc_lnk( e3uw, 'U', 1._wp ) ! lateral boundary conditions222 CALL lbc_lnk( e3v , 'V', 1._wp ) ; CALL lbc_lnk( e3vw, 'V', 1._wp )223 CALL lbc_lnk( fse3u_n(:,:,:) , 'U', 1._wp ) ; CALL lbc_lnk( fse3uw_n(:,:,:), 'U', 1._wp ) ! lateral boundary conditions 224 CALL lbc_lnk( fse3v_n(:,:,:) , 'V', 1._wp ) ; CALL lbc_lnk( fse3vw_n(:,:,:), 'V', 1._wp ) 223 225 ! 224 226 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 225 WHERE( e3u(:,:,jk) == 0._wp ) e3u(:,:,jk) = e3t_0(jk)226 WHERE( e3v(:,:,jk) == 0._wp ) e3v(:,:,jk) = e3t_0(jk)227 WHERE( fse3u_n(:,:,jk) == 0._wp ) fse3u_n(:,:,jk) = e3t_1d(jk) 228 WHERE( fse3v_n(:,:,jk) == 0._wp ) fse3v_n(:,:,jk) = e3t_1d(jk) 227 229 END DO 228 230 END IF 229 231 230 232 IF( iom_varid( inum4, 'gdept', ldstop = .FALSE. ) > 0 ) THEN ! 3D depth of t- and w-level 231 CALL iom_get( inum4, jpdom_data, 'gdept', gdept)232 CALL iom_get( inum4, jpdom_data, 'gdepw', gdepw)233 CALL iom_get( inum4, jpdom_data, 'gdept', fsdept_n(:,:,:) ) 234 CALL iom_get( inum4, jpdom_data, 'gdepw', fsdepw_n(:,:,:) ) 233 235 ELSE ! 2D bottom depth 234 236 CALL iom_get( inum4, jpdom_data, 'hdept', zprt ) … … 236 238 ! 237 239 DO jk = 1, jpk ! deduces the 3D depth 238 gdept(:,:,jk) = gdept_0(jk)239 gdepw(:,:,jk) = gdepw_0(jk)240 fsdept_n(:,:,jk) = gdept_1d(jk) 241 fsdepw_n(:,:,jk) = gdepw_1d(jk) 240 242 END DO 241 243 DO jj = 1, jpj … … 243 245 ik = mbkt(ji,jj) 244 246 IF( ik > 0 ) THEN 245 gdepw(ji,jj,ik+1) = zprw(ji,jj)246 gdept(ji,jj,ik ) = zprt(ji,jj)247 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)247 fsdepw_n(ji,jj,ik+1) = zprw(ji,jj) 248 fsdept_n(ji,jj,ik ) = zprt(ji,jj) 249 fsdept_n(ji,jj,ik+1) = fsdept_n(ji,jj,ik) + fse3t_n(ji,jj,ik) 248 250 ENDIF 249 251 END DO … … 254 256 255 257 IF( ln_zco ) THEN ! Vertical coordinates and scales factors 256 CALL iom_get( inum4, jpdom_unknown, 'gdept_ 0', gdept_0) ! depth257 CALL iom_get( inum4, jpdom_unknown, 'gdepw_ 0', gdepw_0)258 CALL iom_get( inum4, jpdom_unknown, 'e3t_ 0' , e3t_0)259 CALL iom_get( inum4, jpdom_unknown, 'e3w_ 0' , e3w_0)258 CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth 259 CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d ) 260 CALL iom_get( inum4, jpdom_unknown, 'e3t_1d' , e3t_1d ) 261 CALL iom_get( inum4, jpdom_unknown, 'e3w_1d' , e3w_1d ) 260 262 DO jk = 1, jpk 261 e3t (:,:,jk) = e3t_0(jk)! set to the ref. factors262 e3u (:,:,jk) = e3t_0(jk)263 e3v (:,:,jk) = e3t_0(jk)264 e3w (:,:,jk) = e3w_0(jk)265 gdept(:,:,jk) = gdept_0(jk)266 gdepw(:,:,jk) = gdepw_0(jk)263 fse3t_n(:,:,jk) = e3t_1d(jk) ! set to the ref. factors 264 fse3u_n(:,:,jk) = e3t_1d(jk) 265 fse3v_n(:,:,jk) = e3t_1d(jk) 266 fse3w_n(:,:,jk) = e3w_1d(jk) 267 fsdept_n(:,:,jk) = gdept_1d(jk) 268 fsdepw_n(:,:,jk) = gdepw_1d(jk) 267 269 END DO 268 270 ENDIF … … 270 272 !!gm BUG in s-coordinate this does not work! 271 273 ! deepest/shallowest W level Above/Below ~10m 272 zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_ 0) )! ref. depth with tolerance (10% of minimum layer thickness)273 nlb10 = MINLOC( gdepw_ 0, mask = gdepw_0 > zrefdep, dim = 1 )! shallowest W level Below ~10m274 zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_1d) ) ! ref. depth with tolerance (10% of minimum layer thickness) 275 nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 274 276 nla10 = nlb10 - 1 ! deepest W level Above ~10m 275 277 !!gm end bug … … 312 314 WRITE(numout,*) ' Reference z-coordinate depth and scale factors:' 313 315 WRITE(numout, "(9x,' level gdept gdepw e3t e3w ')" ) 314 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_ 0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk )316 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk ) 315 317 ENDIF 316 318 317 319 DO jk = 1, jpk 318 IF( e3w_ 0 (jk) <= 0._wp .OR. e3t_0 (jk) <= 0._wp ) CALL ctl_stop( ' e3w_0 or e3t_0=< 0 ' )319 IF( gdepw_ 0(jk) < 0._wp .OR. gdept_0(jk) < 0._wp ) CALL ctl_stop( ' gdepw_0 or gdept_0< 0 ' )320 IF( e3w_1d (jk) <= 0._wp .OR. e3t_1d (jk) <= 0._wp ) CALL ctl_stop( ' e3w_1d or e3t_1d =< 0 ' ) 321 IF( gdepw_1d(jk) < 0._wp .OR. gdept_1d(jk) < 0._wp ) CALL ctl_stop( ' gdepw_1d or gdept_1d < 0 ' ) 320 322 END DO 321 323 ! ! ============================ -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
r3651 r4292 8 8 !! 3.3 ! 2010-09 (D. Storkey) add ice boundary conditions 9 9 !! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge 10 !! 3.6 ! 2012-01 (C. Rousset) add ice boundary conditions for lim3 10 11 !!---------------------------------------------------------------------- 11 12 #if defined key_bdy … … 27 28 INTEGER, POINTER, DIMENSION(:,:) :: nbr 28 29 INTEGER, POINTER, DIMENSION(:,:) :: nbmap 29 REAL , POINTER, DIMENSION(:,:) :: nbw 30 REAL , POINTER, DIMENSION(:,:) :: nbd 31 REAL , POINTER, DIMENSION(:) :: flagu 32 REAL , POINTER, DIMENSION(:) :: flagv 30 REAL(wp) , POINTER, DIMENSION(:,:) :: nbw 31 REAL(wp) , POINTER, DIMENSION(:,:) :: nbd 32 REAL(wp) , POINTER, DIMENSION(:,:) :: nbdout 33 REAL(wp) , POINTER, DIMENSION(:,:) :: flagu 34 REAL(wp) , POINTER, DIMENSION(:,:) :: flagv 33 35 END TYPE OBC_INDEX 34 36 37 !! Logicals in OBC_DATA structure are true if the chosen algorithm requires this 38 !! field as external data. If true the data can come from external files 39 !! or model initial conditions. If false then no "external" data array 40 !! is required for this field. 41 35 42 TYPE, PUBLIC :: OBC_DATA !: Storage for external data 36 REAL, POINTER, DIMENSION(:) :: ssh 37 REAL, POINTER, DIMENSION(:) :: u2d 38 REAL, POINTER, DIMENSION(:) :: v2d 39 REAL, POINTER, DIMENSION(:,:) :: u3d 40 REAL, POINTER, DIMENSION(:,:) :: v3d 41 REAL, POINTER, DIMENSION(:,:) :: tem 42 REAL, POINTER, DIMENSION(:,:) :: sal 43 INTEGER, DIMENSION(2) :: nread 44 LOGICAL :: ll_ssh 45 LOGICAL :: ll_u2d 46 LOGICAL :: ll_v2d 47 LOGICAL :: ll_u3d 48 LOGICAL :: ll_v3d 49 LOGICAL :: ll_tem 50 LOGICAL :: ll_sal 51 REAL(wp), POINTER, DIMENSION(:) :: ssh 52 REAL(wp), POINTER, DIMENSION(:) :: u2d 53 REAL(wp), POINTER, DIMENSION(:) :: v2d 54 REAL(wp), POINTER, DIMENSION(:,:) :: u3d 55 REAL(wp), POINTER, DIMENSION(:,:) :: v3d 56 REAL(wp), POINTER, DIMENSION(:,:) :: tem 57 REAL(wp), POINTER, DIMENSION(:,:) :: sal 43 58 #if defined key_lim2 44 REAL, POINTER, DIMENSION(:) :: frld 45 REAL, POINTER, DIMENSION(:) :: hicif 46 REAL, POINTER, DIMENSION(:) :: hsnif 59 LOGICAL :: ll_frld 60 LOGICAL :: ll_hicif 61 LOGICAL :: ll_hsnif 62 REAL(wp), POINTER, DIMENSION(:) :: frld 63 REAL(wp), POINTER, DIMENSION(:) :: hicif 64 REAL(wp), POINTER, DIMENSION(:) :: hsnif 65 #elif defined key_lim3 66 LOGICAL :: ll_a_i 67 LOGICAL :: ll_ht_i 68 LOGICAL :: ll_ht_s 69 REAL, POINTER, DIMENSION(:,:) :: a_i !: now ice leads fraction climatology 70 REAL, POINTER, DIMENSION(:,:) :: ht_i !: Now ice thickness climatology 71 REAL, POINTER, DIMENSION(:,:) :: ht_s !: now snow thickness 47 72 #endif 48 73 END TYPE OBC_DATA … … 63 88 INTEGER :: nn_volctl !: = 0 the total volume will have the variability of the surface Flux E-P 64 89 ! ! = 1 the volume will be constant during all the integration. 65 INTEGER, DIMENSION(jp_bdy) :: nn_dyn2d! Choice of boundary condition for barotropic variables (U,V,SSH)66 INTEGER, DIMENSION(jp_bdy) :: nn_dyn2d_dta!: = 0 use the initial state as bdy dta ;90 CHARACTER(len=20), DIMENSION(jp_bdy) :: cn_dyn2d ! Choice of boundary condition for barotropic variables (U,V,SSH) 91 INTEGER, DIMENSION(jp_bdy) :: nn_dyn2d_dta !: = 0 use the initial state as bdy dta ; 67 92 !: = 1 read it in a NetCDF file 68 93 !: = 2 read tidal harmonic forcing from a NetCDF file 69 94 !: = 3 read external data AND tidal harmonic forcing from NetCDF files 70 INTEGER, DIMENSION(jp_bdy) :: nn_dyn3d! Choice of boundary condition for baroclinic velocities71 INTEGER, DIMENSION(jp_bdy) :: nn_dyn3d_dta!: = 0 use the initial state as bdy dta ;95 CHARACTER(len=20), DIMENSION(jp_bdy) :: cn_dyn3d ! Choice of boundary condition for baroclinic velocities 96 INTEGER, DIMENSION(jp_bdy) :: nn_dyn3d_dta !: = 0 use the initial state as bdy dta ; 72 97 !: = 1 read it in a NetCDF file 73 INTEGER, DIMENSION(jp_bdy) :: nn_tra! Choice of boundary condition for active tracers (T and S)74 INTEGER, DIMENSION(jp_bdy) :: nn_tra_dta!: = 0 use the initial state as bdy dta ;98 CHARACTER(len=20), DIMENSION(jp_bdy) :: cn_tra ! Choice of boundary condition for active tracers (T and S) 99 INTEGER, DIMENSION(jp_bdy) :: nn_tra_dta !: = 0 use the initial state as bdy dta ; 75 100 !: = 1 read it in a NetCDF file 76 101 LOGICAL, DIMENSION(jp_bdy) :: ln_tra_dmp !: =T Tracer damping 77 102 LOGICAL, DIMENSION(jp_bdy) :: ln_dyn3d_dmp !: =T Baroclinic velocity damping 78 REAL, DIMENSION(jp_bdy) :: rn_time_dmp !: Damping time scale in days 103 REAL(wp), DIMENSION(jp_bdy) :: rn_time_dmp !: Damping time scale in days 104 REAL(wp), DIMENSION(jp_bdy) :: rn_time_dmp_out !: Damping time scale in days at radiation outflow points 79 105 80 106 #if defined key_lim2 81 INTEGER, DIMENSION(jp_bdy) :: nn_ice_lim2! Choice of boundary condition for sea ice variables82 INTEGER, DIMENSION(jp_bdy) :: nn_ice_lim2_dta!: = 0 use the initial state as bdy dta ;83 !: = 1 read it in a NetCDF file107 CHARACTER(len=20), DIMENSION(jp_bdy) :: nn_ice_lim2 ! Choice of boundary condition for sea ice variables 108 INTEGER, DIMENSION(jp_bdy) :: nn_ice_lim2_dta !: = 0 use the initial state as bdy dta ; 109 !: = 1 read it in a NetCDF file 84 110 #endif 85 111 ! … … 88 114 !! Global variables 89 115 !!---------------------------------------------------------------------- 90 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: bdytmask !: Mask defining computational domain at T-points91 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: bdyumask !: Mask defining computational domain at U-points92 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: bdyvmask !: Mask defining computational domain at V-points116 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET :: bdytmask !: Mask defining computational domain at T-points 117 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET :: bdyumask !: Mask defining computational domain at U-points 118 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET :: bdyvmask !: Mask defining computational domain at V-points 93 119 94 120 REAL(wp) :: bdysurftot !: Lateral surface of unstructured open boundary 95 121 96 REAL(wp), POINTER, DIMENSION(:,:) :: pssh !:97 REAL(wp), POINTER, DIMENSION(:,:) :: phur !:98 REAL(wp), POINTER, DIMENSION(:,:) :: phvr !: Pointers for barotropic fields99 REAL(wp), POINTER, DIMENSION(:,:) :: pu 2d!:100 REAL(wp), POINTER, DIMENSION(:,:) :: pv 2d!:122 REAL(wp), POINTER, DIMENSION(:,:) :: pssh !: 123 REAL(wp), POINTER, DIMENSION(:,:) :: phur !: 124 REAL(wp), POINTER, DIMENSION(:,:) :: phvr !: Pointers for barotropic fields 125 REAL(wp), POINTER, DIMENSION(:,:) :: pub2d, pun2d, pua2d !: 126 REAL(wp), POINTER, DIMENSION(:,:) :: pvb2d, pvn2d, pva2d !: 101 127 102 128 !!---------------------------------------------------------------------- … … 109 135 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET :: dta_global2 !: workspace for reading in global data arrays (struct. bdy) 110 136 TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET :: idx_bdy !: bdy indices (local process) 111 TYPE(OBC_DATA) , DIMENSION(jp_bdy) 137 TYPE(OBC_DATA) , DIMENSION(jp_bdy), TARGET :: dta_bdy !: bdy external data (local process) 112 138 113 139 !!---------------------------------------------------------------------- … … 125 151 !!---------------------------------------------------------------------- 126 152 ! 127 ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj), 153 ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj), & 128 154 & STAT=bdy_oce_alloc ) 129 155 ! -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_par.F90
r3294 r4292 23 23 # endif 24 24 INTEGER, PUBLIC, PARAMETER :: jp_bdy = 10 !: Maximum number of bdy sets 25 INTEGER, PUBLIC, PARAMETER :: jpbtime = 1000 !: Max number of time dumps per file26 25 INTEGER, PUBLIC, PARAMETER :: jpbgrd = 3 !: Number of horizontal grid types used (T, U, V) 27 26 28 !! Flags for choice of schemes29 INTEGER, PUBLIC, PARAMETER :: jp_none = 0 !: Flag for no open boundary condition30 INTEGER, PUBLIC, PARAMETER :: jp_frs = 1 !: Flag for Flow Relaxation Scheme31 INTEGER, PUBLIC, PARAMETER :: jp_flather = 2 !: Flag for Flather32 27 #else 33 28 !!---------------------------------------------------------------------- -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
r4230 r4292 11 11 !! 3.3 ! 2010-09 (D.Storkey) add ice boundary conditions 12 12 !! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge 13 !! 3.6 ! 2012-01 (C. Rousset) add ice boundary conditions for lim3 13 14 !!---------------------------------------------------------------------- 14 15 #if defined key_bdy … … 29 30 USE iom ! IOM library 30 31 USE in_out_manager ! I/O logical units 32 USE dynspg_oce, ONLY: lk_dynspg_ts ! Split-explicit free surface flag 31 33 #if defined key_lim2 32 34 USE ice_2 35 #elif defined key_lim3 36 USE par_ice 37 USE ice 38 USE limcat_1D ! redistribute ice input into categories 33 39 #endif 34 40 USE sbcapr … … 49 55 50 56 TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:) :: nbmap_ptr ! array of pointers to nbmap 57 58 #if defined key_lim3 59 LOGICAL :: ll_bdylim3 ! determine whether ice input is lim2 (F) or lim3 (T) type 60 INTEGER :: jfld_hti, jfld_hts, jfld_ai ! indices of ice thickness, snow thickness and concentration in bf structure 61 #endif 51 62 52 63 # include "domzgr_substitute.h90" … … 77 88 ! etc. 78 89 !! 79 INTEGER :: ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd ! local indices90 INTEGER :: ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd, jl ! local indices 80 91 INTEGER, DIMENSION(jpbgrd) :: ilen1 81 92 INTEGER, POINTER, DIMENSION(:) :: nblen, nblenrim ! short cuts 93 TYPE(OBC_DATA), POINTER :: dta ! short cut 82 94 !! 83 95 !!--------------------------------------------------------------------------- … … 91 103 ! Calculate depth-mean currents 92 104 !----------------------------- 93 CALL wrk_alloc(jpi,jpj,pu2d,pv2d) 94 95 pu2d(:,:) = 0.e0 96 pv2d(:,:) = 0.e0 97 105 CALL wrk_alloc(jpi,jpj,pun2d,pvn2d) 106 107 pun2d(:,:) = 0.e0 108 pvn2d(:,:) = 0.e0 98 109 DO ik = 1, jpkm1 !! Vertically integrated momentum trends 99 pu 2d(:,:) = pu2d(:,:) + fse3u(:,:,ik) * umask(:,:,ik) * un(:,:,ik)100 pv 2d(:,:) = pv2d(:,:) + fse3v(:,:,ik) * vmask(:,:,ik) * vn(:,:,ik)110 pun2d(:,:) = pun2d(:,:) + fse3u(:,:,ik) * umask(:,:,ik) * un(:,:,ik) 111 pvn2d(:,:) = pvn2d(:,:) + fse3v(:,:,ik) * vmask(:,:,ik) * vn(:,:,ik) 101 112 END DO 102 pu 2d(:,:) = pu2d(:,:) * hur(:,:)103 pv 2d(:,:) = pv2d(:,:) * hvr(:,:)113 pun2d(:,:) = pun2d(:,:) * hur(:,:) 114 pvn2d(:,:) = pvn2d(:,:) * hvr(:,:) 104 115 105 116 DO ib_bdy = 1, nb_bdy … … 107 118 nblen => idx_bdy(ib_bdy)%nblen 108 119 nblenrim => idx_bdy(ib_bdy)%nblenrim 109 110 IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 0 ) THEN 120 dta => dta_bdy(ib_bdy) 121 122 IF( nn_dyn2d_dta(ib_bdy) .eq. 0 ) THEN 111 123 ilen1(:) = nblen(:) 112 igrd = 1 113 DO ib = 1, ilen1(igrd) 114 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 115 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 116 dta_bdy(ib_bdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1) 117 END DO 118 igrd = 2 119 DO ib = 1, ilen1(igrd) 120 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 121 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 122 dta_bdy(ib_bdy)%u2d(ib) = pu2d(ii,ij) * umask(ii,ij,1) 123 END DO 124 igrd = 3 125 DO ib = 1, ilen1(igrd) 126 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 127 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 128 dta_bdy(ib_bdy)%v2d(ib) = pv2d(ii,ij) * vmask(ii,ij,1) 129 END DO 130 ENDIF 131 132 IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 133 ilen1(:) = nblen(:) 134 igrd = 2 135 DO ib = 1, ilen1(igrd) 136 DO ik = 1, jpkm1 124 IF( dta%ll_ssh ) THEN 125 igrd = 1 126 DO ib = 1, ilen1(igrd) 137 127 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 138 128 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 139 dta_bdy(ib_bdy)% u3d(ib,ik) = ( un(ii,ij,ik) - pu2d(ii,ij) ) * umask(ii,ij,ik)140 END DO 141 END DO142 igrd = 3143 DO ib = 1, ilen1(igrd)144 DO i k = 1, jpkm1129 dta_bdy(ib_bdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1) 130 END DO 131 END IF 132 IF( dta%ll_u2d ) THEN 133 igrd = 2 134 DO ib = 1, ilen1(igrd) 145 135 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 146 136 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 147 dta_bdy(ib_bdy)%v3d(ib,ik) = ( vn(ii,ij,ik) - pv2d(ii,ij) ) * vmask(ii,ij,ik) 148 END DO 149 END DO 150 ENDIF 151 152 IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 0 ) THEN 153 ilen1(:) = nblen(:) 154 igrd = 1 ! Everything is at T-points here 155 DO ib = 1, ilen1(igrd) 156 DO ik = 1, jpkm1 137 dta_bdy(ib_bdy)%u2d(ib) = pun2d(ii,ij) * umask(ii,ij,1) 138 END DO 139 END IF 140 IF( dta%ll_v2d ) THEN 141 igrd = 3 142 DO ib = 1, ilen1(igrd) 157 143 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 158 144 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 159 dta_bdy(ib_bdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik) 160 dta_bdy(ib_bdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik) 145 dta_bdy(ib_bdy)%v2d(ib) = pvn2d(ii,ij) * vmask(ii,ij,1) 146 END DO 147 END IF 148 ENDIF 149 150 IF( nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 151 ilen1(:) = nblen(:) 152 IF( dta%ll_u3d ) THEN 153 igrd = 2 154 DO ib = 1, ilen1(igrd) 155 DO ik = 1, jpkm1 156 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 157 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 158 dta_bdy(ib_bdy)%u3d(ib,ik) = ( un(ii,ij,ik) - pun2d(ii,ij) ) * umask(ii,ij,ik) 159 END DO 160 END DO 161 END IF 162 IF( dta%ll_v3d ) THEN 163 igrd = 3 164 DO ib = 1, ilen1(igrd) 165 DO ik = 1, jpkm1 166 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 167 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 168 dta_bdy(ib_bdy)%v3d(ib,ik) = ( vn(ii,ij,ik) - pvn2d(ii,ij) ) * vmask(ii,ij,ik) 169 END DO 170 END DO 171 END IF 172 ENDIF 173 174 IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN 175 ilen1(:) = nblen(:) 176 IF( dta%ll_tem ) THEN 177 igrd = 1 178 DO ib = 1, ilen1(igrd) 179 DO ik = 1, jpkm1 180 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 181 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 182 dta_bdy(ib_bdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik) 183 END DO 184 END DO 185 END IF 186 IF( dta%ll_sal ) THEN 187 igrd = 1 188 DO ib = 1, ilen1(igrd) 189 DO ik = 1, jpkm1 190 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 191 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 192 dta_bdy(ib_bdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik) 193 END DO 194 END DO 195 END IF 196 ENDIF 197 198 #if defined key_lim2 199 IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 200 ilen1(:) = nblen(:) 201 IF( dta%ll_frld ) THEN 202 igrd = 1 203 DO ib = 1, ilen1(igrd) 204 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 205 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 206 dta_bdy(ib_bdy)%frld(ib) = frld(ii,ij) * tmask(ii,ij,1) 207 END DO 208 END IF 209 IF( dta%ll_hicif ) THEN 210 igrd = 1 211 DO ib = 1, ilen1(igrd) 212 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 213 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 214 dta_bdy(ib_bdy)%hicif(ib) = hicif(ii,ij) * tmask(ii,ij,1) 215 END DO 216 END IF 217 IF( dta%ll_hsnif ) THEN 218 igrd = 1 219 DO ib = 1, ilen1(igrd) 220 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 221 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 222 dta_bdy(ib_bdy)%hsnif(ib) = hsnif(ii,ij) * tmask(ii,ij,1) 223 END DO 224 END IF 225 ENDIF 226 #elif defined key_lim3 227 IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN 228 ilen1(:) = nblen(:) 229 IF( dta%ll_a_i ) THEN 230 igrd = 1 231 DO jl = 1, jpl 232 DO ib = 1, ilen1(igrd) 233 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 234 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 235 dta_bdy(ib_bdy)%a_i (ib,jl) = a_i(ii,ij,jl) * tmask(ii,ij,1) 236 END DO 161 237 END DO 162 END DO 163 ENDIF 164 165 #if defined key_lim2 166 IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 167 ilen1(:) = nblen(:) 168 igrd = 1 ! Everything is at T-points here 169 DO ib = 1, ilen1(igrd) 170 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 171 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 172 dta_bdy(ib_bdy)%frld(ib) = frld(ii,ij) * tmask(ii,ij,1) 173 dta_bdy(ib_bdy)%hicif(ib) = hicif(ii,ij) * tmask(ii,ij,1) 174 dta_bdy(ib_bdy)%hsnif(ib) = hsnif(ii,ij) * tmask(ii,ij,1) 175 END DO 238 ENDIF 239 IF( dta%ll_ht_i ) THEN 240 igrd = 1 241 DO jl = 1, jpl 242 DO ib = 1, ilen1(igrd) 243 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 244 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 245 dta_bdy(ib_bdy)%ht_i (ib,jl) = ht_i(ii,ij,jl) * tmask(ii,ij,1) 246 END DO 247 END DO 248 ENDIF 249 IF( dta%ll_ht_s ) THEN 250 igrd = 1 251 DO jl = 1, jpl 252 DO ib = 1, ilen1(igrd) 253 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 254 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 255 dta_bdy(ib_bdy)%ht_s (ib,jl) = ht_s(ii,ij,jl) * tmask(ii,ij,1) 256 END DO 257 END DO 258 ENDIF 176 259 ENDIF 177 260 #endif … … 179 262 ENDDO ! ib_bdy 180 263 181 CALL wrk_dealloc(jpi,jpj,pu 2d,pv2d)264 CALL wrk_dealloc(jpi,jpj,pun2d,pvn2d) 182 265 183 266 ENDIF ! kt .eq. nit000 … … 188 271 jstart = 1 189 272 DO ib_bdy = 1, nb_bdy 273 dta => dta_bdy(ib_bdy) 190 274 IF( nn_dta(ib_bdy) .eq. 1 ) THEN ! skip this bit if no external data required 191 275 … … 193 277 ! Update barotropic boundary conditions only 194 278 ! jit is optional argument for fld_read and bdytide_update 195 IF( nn_dyn2d(ib_bdy) .gt. 0) THEN279 IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN 196 280 IF( nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 197 dta_bdy(ib_bdy)%ssh(:) = 0.0198 dta_bdy(ib_bdy)%u2d(:) = 0.0199 dta_bdy(ib_bdy)%v2d(:) = 0.0281 IF( dta%ll_ssh ) dta%ssh(:) = 0.0 282 IF( dta%ll_u2d ) dta%u2d(:) = 0.0 283 IF( dta%ll_u3d ) dta%v2d(:) = 0.0 200 284 ENDIF 201 IF (nn_tra(ib_bdy).ne.4) THEN 202 IF( nn_dyn2d_dta(ib_bdy) .EQ. 1 .OR. nn_dyn2d_dta(ib_bdy) .EQ. 3 .OR. & 203 & (ln_full_vel_array(ib_bdy) .AND. nn_dyn3d_dta(ib_bdy).eq.1) )THEN 204 205 ! For the runoff case, no need to update the forcing (already done in the baroclinic part) 206 jend = nb_bdy_fld(ib_bdy) 207 IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend - 2 285 IF (cn_tra(ib_bdy) /= 'runoff') THEN 286 IF( nn_dyn2d_dta(ib_bdy) .EQ. 1 .OR. nn_dyn2d_dta(ib_bdy) .EQ. 3 ) THEN 287 288 jend = jstart + dta%nread(2) - 1 208 289 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), & 209 290 & kit=jit, kt_offset=time_offset ) 210 IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend + 2 211 212 ! If full velocities in boundary data then split into barotropic and baroclinic data 291 292 ! If full velocities in boundary data then extract barotropic velocities from 3D fields 213 293 IF( ln_full_vel_array(ib_bdy) .AND. & 214 294 & ( nn_dyn2d_dta(ib_bdy) .EQ. 1 .OR. nn_dyn2d_dta(ib_bdy) .EQ. 3 .OR. & … … 216 296 217 297 igrd = 2 ! zonal velocity 218 dta _bdy(ib_bdy)%u2d(:) = 0.0298 dta%u2d(:) = 0.0 219 299 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 220 300 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 221 301 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 222 302 DO ik = 1, jpkm1 223 dta _bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) &224 & + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta _bdy(ib_bdy)%u3d(ib,ik)303 dta%u2d(ib) = dta%u2d(ib) & 304 & + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik) 225 305 END DO 226 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 227 DO ik = 1, jpkm1 228 dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib) 229 END DO 306 dta%u2d(ib) = dta%u2d(ib) * hur(ii,ij) 230 307 END DO 231 308 igrd = 3 ! meridional velocity 232 dta _bdy(ib_bdy)%v2d(:) = 0.0309 dta%v2d(:) = 0.0 233 310 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 234 311 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 235 312 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 236 313 DO ik = 1, jpkm1 237 dta _bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) &238 & + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta _bdy(ib_bdy)%v3d(ib,ik)314 dta%v2d(ib) = dta%v2d(ib) & 315 & + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik) 239 316 END DO 240 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 241 DO ik = 1, jpkm1 242 dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib) 243 END DO 317 dta%v2d(ib) = dta%v2d(ib) * hvr(ii,ij) 244 318 END DO 245 319 ENDIF 246 320 ENDIF 247 321 IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 248 CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta _bdy(ib_bdy), td=tides(ib_bdy), &322 CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta, td=tides(ib_bdy), & 249 323 & jit=jit, time_offset=time_offset ) 250 324 ENDIF … … 252 326 ENDIF 253 327 ELSE 254 IF ( nn_tra(ib_bdy).eq.4) then ! runoff condition328 IF (cn_tra(ib_bdy) == 'runoff') then ! runoff condition 255 329 jend = nb_bdy_fld(ib_bdy) 256 330 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), & … … 261 335 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 262 336 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 263 dta _bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )337 dta%u2d(ib) = dta%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) ) 264 338 END DO 265 339 ! … … 268 342 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 269 343 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 270 dta _bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )344 dta%v2d(ib) = dta%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) ) 271 345 END DO 272 346 ELSE 273 IF( nn_dyn2d (ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays274 dta_bdy(ib_bdy)%ssh(:) = 0.0275 dta_bdy(ib_bdy)%u2d(:) = 0.0276 dta_bdy(ib_bdy)%v2d(:) = 0.0347 IF( nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 348 IF( dta%ll_ssh ) dta%ssh(:) = 0.0 349 IF( dta%ll_u2d ) dta%u2d(:) = 0.0 350 IF( dta%ll_v2d ) dta%v2d(:) = 0.0 277 351 ENDIF 278 IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data279 jend = nb_bdy_fld(ib_bdy)352 IF( dta%nread(1) .gt. 0 ) THEN ! update external data 353 jend = jstart + dta%nread(1) - 1 280 354 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), & 281 355 & map=nbmap_ptr(jstart:jend), kt_offset=time_offset ) … … 286 360 & nn_dyn3d_dta(ib_bdy) .EQ. 1 ) ) THEN 287 361 igrd = 2 ! zonal velocity 288 dta _bdy(ib_bdy)%u2d(:) = 0.0362 dta%u2d(:) = 0.0 289 363 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 290 364 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 291 365 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 292 366 DO ik = 1, jpkm1 293 dta _bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) &294 & + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta _bdy(ib_bdy)%u3d(ib,ik)367 dta%u2d(ib) = dta%u2d(ib) & 368 & + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik) 295 369 END DO 296 dta _bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij)370 dta%u2d(ib) = dta%u2d(ib) * hur(ii,ij) 297 371 DO ik = 1, jpkm1 298 dta _bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib)372 dta%u3d(ib,ik) = dta%u3d(ib,ik) - dta%u2d(ib) 299 373 END DO 300 374 END DO 301 375 igrd = 3 ! meridional velocity 302 dta _bdy(ib_bdy)%v2d(:) = 0.0376 dta%v2d(:) = 0.0 303 377 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 304 378 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 305 379 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 306 380 DO ik = 1, jpkm1 307 dta _bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) &308 & + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta _bdy(ib_bdy)%v3d(ib,ik)381 dta%v2d(ib) = dta%v2d(ib) & 382 & + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik) 309 383 END DO 310 dta _bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij)384 dta%v2d(ib) = dta%v2d(ib) * hvr(ii,ij) 311 385 DO ik = 1, jpkm1 312 dta _bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib)386 dta%v3d(ib,ik) = dta%v3d(ib,ik) - dta%v2d(ib) 313 387 END DO 314 388 END DO 315 389 ENDIF 316 IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 317 CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), & 318 & td=tides(ib_bdy), time_offset=time_offset ) 319 ENDIF 320 ENDIF 321 ENDIF 322 jstart = jend+1 390 391 ENDIF 392 #if defined key_lim3 393 IF( .NOT. ll_bdylim3 .AND. nn_ice_lim(ib_bdy) > 0 .AND. nn_ice_lim_dta(ib_bdy) == 1 ) THEN ! bdy ice input (case input is lim2 type) 394 CALL lim_cat_1D ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), & 395 & dta_bdy(ib_bdy)%ht_i, dta_bdy(ib_bdy)%ht_s, dta_bdy(ib_bdy)%a_i ) 396 ENDIF 397 #endif 398 ENDIF 399 jstart = jstart + dta%nread(1) 323 400 END IF ! nn_dta(ib_bdy) = 1 324 401 END DO ! ib_bdy 325 402 403 ! bg jchanut tschanges 404 #if defined key_tide 405 ! Add tides if not split-explicit free surface else this is done in ts loop 406 IF (.NOT.lk_dynspg_ts) CALL bdy_dta_tides( kt=kt, time_offset=time_offset ) 407 #endif 408 ! end jchanut tschanges 409 326 410 IF ( ln_apr_obc ) THEN 327 411 DO ib_bdy = 1, nb_bdy 328 IF ( nn_tra(ib_bdy).NE.4)THEN412 IF (cn_tra(ib_bdy) /= 'runoff')THEN 329 413 igrd = 1 ! meridional velocity 330 414 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) … … 349 433 !! for open boundary conditions 350 434 !! 351 !! ** Method : Use fldread.F90435 !! ** Method : 352 436 !! 353 437 !!---------------------------------------------------------------------- … … 362 446 ! =F => baroclinic velocities in 3D boundary data 363 447 INTEGER :: ilen_global ! Max length required for global bdy dta arrays 364 INTEGER, DIMENSION(jpbgrd) :: ilen0 ! size of local arrays365 448 INTEGER, ALLOCATABLE, DIMENSION(:) :: ilen1, ilen3 ! size of 1st and 3rd dimensions of local arrays 366 449 INTEGER, ALLOCATABLE, DIMENSION(:) :: ibdy ! bdy set for a particular jfld 367 450 INTEGER, ALLOCATABLE, DIMENSION(:) :: igrid ! index for grid type (1,2,3 = T,U,V) 368 451 INTEGER, POINTER, DIMENSION(:) :: nblen, nblenrim ! short cuts 452 TYPE(OBC_DATA), POINTER :: dta ! short cut 453 #if defined key_lim3 454 INTEGER, DIMENSION(3) :: zdimsz ! number of elements in each of the 4 dimensions (i.e. i,j,t,ice-cat) for an array 455 INTEGER :: zndims ! number of dimensions in an array (i.e. 3 = wo ice cat; 4 = w ice cat) 456 INTEGER :: inum,id1 ! local integer 457 #endif 369 458 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: blf_i ! array of namelist information structures 370 459 TYPE(FLD_N) :: bn_tem, bn_sal, bn_u3d, bn_v3d ! … … 372 461 #if defined key_lim2 373 462 TYPE(FLD_N) :: bn_frld, bn_hicif, bn_hsnif ! 463 #elif defined key_lim3 464 TYPE(FLD_N) :: bn_a_i, bn_ht_i, bn_ht_s 374 465 #endif 375 466 NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d 376 467 #if defined key_lim2 377 468 NAMELIST/nambdy_dta/ bn_frld, bn_hicif, bn_hsnif 469 #elif defined key_lim3 470 NAMELIST/nambdy_dta/ bn_a_i, bn_ht_i, bn_ht_s 378 471 #endif 379 472 NAMELIST/nambdy_dta/ ln_full_vel … … 392 485 ,nn_dyn3d_dta(ib_bdy) & 393 486 ,nn_tra_dta(ib_bdy) & 394 #if defined key_lim2395 ,nn_ice_lim2_dta(ib_bdy) &487 #if ( defined key_lim2 || defined key_lim3 ) 488 ,nn_ice_lim_dta(ib_bdy) & 396 489 #endif 397 490 ) … … 404 497 nb_bdy_fld(:) = 0 405 498 DO ib_bdy = 1, nb_bdy 406 IF( nn_dyn2d(ib_bdy) .gt. 0.and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN499 IF( cn_dyn2d(ib_bdy) /= 'none' .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN 407 500 nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3 408 501 ENDIF 409 IF( nn_dyn3d(ib_bdy) .gt. 0.and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) THEN502 IF( cn_dyn3d(ib_bdy) /= 'none' .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) THEN 410 503 nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2 411 504 ENDIF 412 IF( nn_tra(ib_bdy) .gt. 0.and. nn_tra_dta(ib_bdy) .eq. 1 ) THEN505 IF( cn_tra(ib_bdy) /= 'none' .and. nn_tra_dta(ib_bdy) .eq. 1 ) THEN 413 506 nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2 414 507 ENDIF 415 #if defined key_lim2416 IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1 ) THEN508 #if ( defined key_lim2 || defined key_lim3 ) 509 IF( cn_ice_lim(ib_bdy) /= 'none' .and. nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN 417 510 nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3 418 511 ENDIF … … 458 551 nblen => idx_bdy(ib_bdy)%nblen 459 552 nblenrim => idx_bdy(ib_bdy)%nblenrim 553 dta => dta_bdy(ib_bdy) 554 dta%nread(2) = 0 460 555 461 556 ! Only read in necessary fields for this set. 462 557 ! Important that barotropic variables come first. 463 IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN 464 465 IF( nn_dyn2d(ib_bdy) .ne. jp_frs .and. nn_tra(ib_bdy) .ne. 4 ) THEN ! runoff condition : no ssh reading 558 IF( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN 559 560 IF( dta%ll_ssh ) THEN 561 if(lwp) write(numout,*) '++++++ reading in ssh field' 466 562 jfld = jfld + 1 467 563 blf_i(jfld) = bn_ssh … … 470 566 ilen1(jfld) = nblen(igrid(jfld)) 471 567 ilen3(jfld) = 1 472 ENDIF 473 474 IF( .not. ln_full_vel_array(ib_bdy) ) THEN 568 dta%nread(2) = dta%nread(2) + 1 569 ENDIF 570 571 IF( dta%ll_u2d .and. .not. ln_full_vel_array(ib_bdy) ) THEN 572 if(lwp) write(numout,*) '++++++ reading in u2d field' 475 573 jfld = jfld + 1 476 574 blf_i(jfld) = bn_u2d … … 479 577 ilen1(jfld) = nblen(igrid(jfld)) 480 578 ilen3(jfld) = 1 481 579 dta%nread(2) = dta%nread(2) + 1 580 ENDIF 581 582 IF( dta%ll_v2d .and. .not. ln_full_vel_array(ib_bdy) ) THEN 583 if(lwp) write(numout,*) '++++++ reading in v2d field' 482 584 jfld = jfld + 1 483 585 blf_i(jfld) = bn_v2d … … 486 588 ilen1(jfld) = nblen(igrid(jfld)) 487 589 ilen3(jfld) = 1 488 ENDIF 489 490 ENDIF 491 492 ! baroclinic velocities 493 IF( ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) .or. & 494 & ( ln_full_vel_array(ib_bdy) .and. nn_dyn2d(ib_bdy) .gt. 0 .and. & 495 & ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN 496 497 jfld = jfld + 1 498 blf_i(jfld) = bn_u3d 499 ibdy(jfld) = ib_bdy 500 igrid(jfld) = 2 501 ilen1(jfld) = nblen(igrid(jfld)) 502 ilen3(jfld) = jpk 503 504 jfld = jfld + 1 505 blf_i(jfld) = bn_v3d 506 ibdy(jfld) = ib_bdy 507 igrid(jfld) = 3 508 ilen1(jfld) = nblen(igrid(jfld)) 509 ilen3(jfld) = jpk 590 dta%nread(2) = dta%nread(2) + 1 591 ENDIF 592 593 ENDIF 594 595 ! read 3D velocities if baroclinic velocities require OR if 596 ! barotropic velocities required and ln_full_vel set to .true. 597 IF( nn_dyn3d_dta(ib_bdy) .eq. 1 .or. & 598 & ( ln_full_vel_array(ib_bdy) .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN 599 600 IF( dta%ll_u3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) ) THEN 601 if(lwp) write(numout,*) '++++++ reading in u3d field' 602 jfld = jfld + 1 603 blf_i(jfld) = bn_u3d 604 ibdy(jfld) = ib_bdy 605 igrid(jfld) = 2 606 ilen1(jfld) = nblen(igrid(jfld)) 607 ilen3(jfld) = jpk 608 IF( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) dta%nread(2) = dta%nread(2) + 1 609 ENDIF 610 611 IF( dta%ll_v3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) ) THEN 612 if(lwp) write(numout,*) '++++++ reading in v3d field' 613 jfld = jfld + 1 614 blf_i(jfld) = bn_v3d 615 ibdy(jfld) = ib_bdy 616 igrid(jfld) = 3 617 ilen1(jfld) = nblen(igrid(jfld)) 618 ilen3(jfld) = jpk 619 IF( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) dta%nread(2) = dta%nread(2) + 1 620 ENDIF 510 621 511 622 ENDIF 512 623 513 624 ! temperature and salinity 514 IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 1 ) THEN 515 516 jfld = jfld + 1 517 blf_i(jfld) = bn_tem 518 ibdy(jfld) = ib_bdy 519 igrid(jfld) = 1 520 ilen1(jfld) = nblen(igrid(jfld)) 521 ilen3(jfld) = jpk 522 523 jfld = jfld + 1 524 blf_i(jfld) = bn_sal 525 ibdy(jfld) = ib_bdy 526 igrid(jfld) = 1 527 ilen1(jfld) = nblen(igrid(jfld)) 528 ilen3(jfld) = jpk 625 IF( nn_tra_dta(ib_bdy) .eq. 1 ) THEN 626 627 IF( dta%ll_tem ) THEN 628 if(lwp) write(numout,*) '++++++ reading in tem field' 629 jfld = jfld + 1 630 blf_i(jfld) = bn_tem 631 ibdy(jfld) = ib_bdy 632 igrid(jfld) = 1 633 ilen1(jfld) = nblen(igrid(jfld)) 634 ilen3(jfld) = jpk 635 ENDIF 636 637 IF( dta%ll_sal ) THEN 638 if(lwp) write(numout,*) '++++++ reading in sal field' 639 jfld = jfld + 1 640 blf_i(jfld) = bn_sal 641 ibdy(jfld) = ib_bdy 642 igrid(jfld) = 1 643 ilen1(jfld) = nblen(igrid(jfld)) 644 ilen3(jfld) = jpk 645 ENDIF 529 646 530 647 ENDIF … … 532 649 #if defined key_lim2 533 650 ! sea ice 534 IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1 ) THEN 535 536 jfld = jfld + 1 537 blf_i(jfld) = bn_frld 538 ibdy(jfld) = ib_bdy 539 igrid(jfld) = 1 540 ilen1(jfld) = nblen(igrid(jfld)) 541 ilen3(jfld) = 1 542 543 jfld = jfld + 1 544 blf_i(jfld) = bn_hicif 545 ibdy(jfld) = ib_bdy 546 igrid(jfld) = 1 547 ilen1(jfld) = nblen(igrid(jfld)) 548 ilen3(jfld) = 1 549 550 jfld = jfld + 1 551 blf_i(jfld) = bn_hsnif 552 ibdy(jfld) = ib_bdy 553 igrid(jfld) = 1 554 ilen1(jfld) = nblen(igrid(jfld)) 555 ilen3(jfld) = 1 556 557 ENDIF 651 IF( nn_ice_lim2_dta(ib_bdy) .eq. 1 ) THEN 652 653 IF( dta%ll_frld ) THEN 654 jfld = jfld + 1 655 blf_i(jfld) = bn_frld 656 ibdy(jfld) = ib_bdy 657 igrid(jfld) = 1 658 ilen1(jfld) = nblen(igrid(jfld)) 659 ilen3(jfld) = 1 660 ENDIF 661 662 IF( dta%ll_hicif ) THEN 663 jfld = jfld + 1 664 blf_i(jfld) = bn_hicif 665 ibdy(jfld) = ib_bdy 666 igrid(jfld) = 1 667 ilen1(jfld) = nblen(igrid(jfld)) 668 ilen3(jfld) = 1 669 ENDIF 670 671 IF( dta%ll_hsnif ) THEN 672 jfld = jfld + 1 673 blf_i(jfld) = bn_hsnif 674 ibdy(jfld) = ib_bdy 675 igrid(jfld) = 1 676 ilen1(jfld) = nblen(igrid(jfld)) 677 ilen3(jfld) = 1 678 ENDIF 679 680 ENDIF 681 #elif defined key_lim3 682 ! sea ice 683 IF( nn_ice_lim_dta(ib_bdy) .eq. 1 ) THEN 684 685 ! Test for types of ice input (lim2 or lim3) 686 CALL iom_open ( bn_a_i%clname, inum ) 687 id1 = iom_varid ( inum, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. ) 688 CALL iom_close ( inum ) 689 !CALL fld_clopn ( bn_a_i, nyear, nmonth, nday, ldstop=.TRUE. ) 690 !CALL iom_open ( bn_a_i %clname, inum ) 691 !id1 = iom_varid ( bn_a_i%num, bn_a_i%clvar, kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE. ) 692 IF ( zndims == 4 ) THEN 693 ll_bdylim3 = .TRUE. ! lim3 input 694 ELSE 695 ll_bdylim3 = .FALSE. ! lim2 input 696 ENDIF 697 ! End test 698 699 IF( dta%ll_a_i ) THEN 700 jfld = jfld + 1 701 blf_i(jfld) = bn_a_i 702 ibdy(jfld) = ib_bdy 703 igrid(jfld) = 1 704 ilen1(jfld) = nblen(igrid(jfld)) 705 IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF 706 ENDIF 707 708 IF( dta%ll_ht_i ) THEN 709 jfld = jfld + 1 710 blf_i(jfld) = bn_ht_i 711 ibdy(jfld) = ib_bdy 712 igrid(jfld) = 1 713 ilen1(jfld) = nblen(igrid(jfld)) 714 IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF 715 ENDIF 716 717 IF( dta%ll_ht_s ) THEN 718 jfld = jfld + 1 719 blf_i(jfld) = bn_ht_s 720 ibdy(jfld) = ib_bdy 721 igrid(jfld) = 1 722 ilen1(jfld) = nblen(igrid(jfld)) 723 IF ( ll_bdylim3 ) THEN ; ilen3(jfld)=jpl ; ELSE ; ilen3(jfld)=1 ; ENDIF 724 ENDIF 725 558 726 #endif 559 727 ! Recalculate field counts … … 568 736 ENDIF 569 737 738 dta%nread(1) = nb_bdy_fld(ib_bdy) 739 570 740 ENDIF ! nn_dta .eq. 1 571 741 ENDDO ! ib_bdy … … 596 766 597 767 nblen => idx_bdy(ib_bdy)%nblen 598 nblenrim => idx_bdy(ib_bdy)%nblenrim 599 600 IF (nn_dyn2d(ib_bdy) .gt. 0) THEN 601 IF( nn_dyn2d_dta(ib_bdy) .eq. 0 .or. nn_dyn2d_dta(ib_bdy) .eq. 2 .or. ln_full_vel_array(ib_bdy) ) THEN 602 ilen0(1:3) = nblen(1:3) 603 ALLOCATE( dta_bdy(ib_bdy)%u2d(ilen0(2)) ) 604 ALLOCATE( dta_bdy(ib_bdy)%v2d(ilen0(3)) ) 605 IF ( nn_dyn2d(ib_bdy) .ne. jp_frs .and. (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3) ) THEN 606 jfld = jfld + 1 607 dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1) 768 dta => dta_bdy(ib_bdy) 769 770 if(lwp) then 771 write(numout,*) '++++++ dta%ll_ssh = ',dta%ll_ssh 772 write(numout,*) '++++++ dta%ll_u2d = ',dta%ll_u2d 773 write(numout,*) '++++++ dta%ll_v2d = ',dta%ll_v2d 774 write(numout,*) '++++++ dta%ll_u3d = ',dta%ll_u3d 775 write(numout,*) '++++++ dta%ll_v3d = ',dta%ll_v3d 776 write(numout,*) '++++++ dta%ll_tem = ',dta%ll_tem 777 write(numout,*) '++++++ dta%ll_sal = ',dta%ll_sal 778 endif 779 780 IF ( nn_dyn2d_dta(ib_bdy) .eq. 0 .or. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN 781 if(lwp) write(numout,*) '++++++ dta%ssh/u2d/u3d allocated space' 782 IF( dta%ll_ssh ) ALLOCATE( dta%ssh(nblen(1)) ) 783 IF( dta%ll_u2d ) ALLOCATE( dta%u2d(nblen(2)) ) 784 IF( dta%ll_v2d ) ALLOCATE( dta%v2d(nblen(3)) ) 785 ENDIF 786 IF ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN 787 IF( dta%ll_ssh ) THEN 788 if(lwp) write(numout,*) '++++++ dta%ssh pointing to fnow' 789 jfld = jfld + 1 790 dta%ssh => bf(jfld)%fnow(:,1,1) 791 ENDIF 792 IF ( dta%ll_u2d ) THEN 793 IF ( ln_full_vel_array(ib_bdy) ) THEN 794 if(lwp) write(numout,*) '++++++ dta%u2d allocated space' 795 ALLOCATE( dta%u2d(nblen(2)) ) 608 796 ELSE 609 ALLOCATE( dta_bdy(ib_bdy)%ssh(nblen(1)) ) 610 ENDIF 611 ELSE 612 IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN 613 jfld = jfld + 1 614 dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1) 615 ENDIF 797 if(lwp) write(numout,*) '++++++ dta%u2d pointing to fnow' 798 jfld = jfld + 1 799 dta%u2d => bf(jfld)%fnow(:,1,1) 800 ENDIF 801 ENDIF 802 IF ( dta%ll_v2d ) THEN 803 IF ( ln_full_vel_array(ib_bdy) ) THEN 804 if(lwp) write(numout,*) '++++++ dta%v2d allocated space' 805 ALLOCATE( dta%v2d(nblen(3)) ) 806 ELSE 807 if(lwp) write(numout,*) '++++++ dta%v2d pointing to fnow' 808 jfld = jfld + 1 809 dta%v2d => bf(jfld)%fnow(:,1,1) 810 ENDIF 811 ENDIF 812 ENDIF 813 814 IF ( nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 815 if(lwp) write(numout,*) '++++++ dta%u3d/v3d allocated space' 816 IF( dta%ll_u3d ) ALLOCATE( dta_bdy(ib_bdy)%u3d(nblen(2),jpk) ) 817 IF( dta%ll_v3d ) ALLOCATE( dta_bdy(ib_bdy)%v3d(nblen(3),jpk) ) 818 ENDIF 819 IF ( nn_dyn3d_dta(ib_bdy) .eq. 1 .or. & 820 & ( ln_full_vel_array(ib_bdy) .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN 821 IF ( dta%ll_u3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_u2d ) ) THEN 822 if(lwp) write(numout,*) '++++++ dta%u3d pointing to fnow' 616 823 jfld = jfld + 1 617 dta_bdy(ib_bdy)%u2d => bf(jfld)%fnow(:,1,1) 824 dta_bdy(ib_bdy)%u3d => bf(jfld)%fnow(:,1,:) 825 ENDIF 826 IF ( dta%ll_v3d .or. ( ln_full_vel_array(ib_bdy) .and. dta%ll_v2d ) ) THEN 827 if(lwp) write(numout,*) '++++++ dta%v3d pointing to fnow' 618 828 jfld = jfld + 1 619 dta_bdy(ib_bdy)%v2d => bf(jfld)%fnow(:,1,1) 620 ENDIF 621 ENDIF 622 623 IF ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 624 ilen0(1:3) = nblen(1:3) 625 ALLOCATE( dta_bdy(ib_bdy)%u3d(ilen0(2),jpk) ) 626 ALLOCATE( dta_bdy(ib_bdy)%v3d(ilen0(3),jpk) ) 627 ENDIF 628 IF ( ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ).or. & 629 & ( ln_full_vel_array(ib_bdy) .and. nn_dyn2d(ib_bdy) .gt. 0 .and. & 630 & ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN 631 jfld = jfld + 1 632 dta_bdy(ib_bdy)%u3d => bf(jfld)%fnow(:,1,:) 633 jfld = jfld + 1 634 dta_bdy(ib_bdy)%v3d => bf(jfld)%fnow(:,1,:) 635 ENDIF 636 637 IF (nn_tra(ib_bdy) .gt. 0) THEN 638 IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN 639 ilen0(1:3) = nblen(1:3) 640 ALLOCATE( dta_bdy(ib_bdy)%tem(ilen0(1),jpk) ) 641 ALLOCATE( dta_bdy(ib_bdy)%sal(ilen0(1),jpk) ) 642 ELSE 829 dta_bdy(ib_bdy)%v3d => bf(jfld)%fnow(:,1,:) 830 ENDIF 831 ENDIF 832 833 IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN 834 if(lwp) write(numout,*) '++++++ dta%tem/sal allocated space' 835 IF( dta%ll_tem ) ALLOCATE( dta_bdy(ib_bdy)%tem(nblen(1),jpk) ) 836 IF( dta%ll_sal ) ALLOCATE( dta_bdy(ib_bdy)%sal(nblen(1),jpk) ) 837 ELSE 838 IF( dta%ll_tem ) THEN 839 if(lwp) write(numout,*) '++++++ dta%tem pointing to fnow' 643 840 jfld = jfld + 1 644 841 dta_bdy(ib_bdy)%tem => bf(jfld)%fnow(:,1,:) 842 ENDIF 843 IF( dta%ll_sal ) THEN 844 if(lwp) write(numout,*) '++++++ dta%sal pointing to fnow' 645 845 jfld = jfld + 1 646 846 dta_bdy(ib_bdy)%sal => bf(jfld)%fnow(:,1,:) … … 649 849 650 850 #if defined key_lim2 651 IF (nn_ice_lim 2(ib_bdy) .gt. 0) THEN851 IF (nn_ice_lim(ib_bdy) .gt. 0) THEN 652 852 IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 653 ilen0(1:3) = nblen(1:3) 654 ALLOCATE( dta_bdy(ib_bdy)%frld(ilen0(1)) ) 655 ALLOCATE( dta_bdy(ib_bdy)%hicif(ilen0(1)) ) 656 ALLOCATE( dta_bdy(ib_bdy)%hsnif(ilen0(1)) ) 853 ALLOCATE( dta_bdy(ib_bdy)%frld(nblen(1)) ) 854 ALLOCATE( dta_bdy(ib_bdy)%hicif(nblen(1)) ) 855 ALLOCATE( dta_bdy(ib_bdy)%hsnif(nblen(1)) ) 657 856 ELSE 658 857 jfld = jfld + 1 … … 662 861 jfld = jfld + 1 663 862 dta_bdy(ib_bdy)%hsnif => bf(jfld)%fnow(:,1,1) 863 ENDIF 864 ENDIF 865 #elif defined key_lim3 866 IF (nn_ice_lim(ib_bdy) .gt. 0) THEN 867 IF( nn_ice_lim_dta(ib_bdy) .eq. 0 ) THEN 868 ALLOCATE( dta_bdy(ib_bdy)%a_i (nblen(1),jpl) ) 869 ALLOCATE( dta_bdy(ib_bdy)%ht_i(nblen(1),jpl) ) 870 ALLOCATE( dta_bdy(ib_bdy)%ht_s(nblen(1),jpl) ) 871 ELSE 872 IF ( ll_bdylim3 ) THEN ! case input is lim3 type 873 jfld = jfld + 1 874 dta_bdy(ib_bdy)%a_i => bf(jfld)%fnow(:,1,:) 875 jfld = jfld + 1 876 dta_bdy(ib_bdy)%ht_i => bf(jfld)%fnow(:,1,:) 877 jfld = jfld + 1 878 dta_bdy(ib_bdy)%ht_s => bf(jfld)%fnow(:,1,:) 879 ELSE ! case input is lim2 type 880 jfld_ai = jfld + 1 881 jfld_hti = jfld + 2 882 jfld_hts = jfld + 3 883 jfld = jfld + 3 884 ALLOCATE( dta_bdy(ib_bdy)%a_i (nblen(1),jpl) ) 885 ALLOCATE( dta_bdy(ib_bdy)%ht_i(nblen(1),jpl) ) 886 ALLOCATE( dta_bdy(ib_bdy)%ht_s(nblen(1),jpl) ) 887 dta_bdy(ib_bdy)%a_i (:,:) = 0.0 888 dta_bdy(ib_bdy)%ht_i(:,:) = 0.0 889 dta_bdy(ib_bdy)%ht_s(:,:) = 0.0 890 ENDIF 891 664 892 ENDIF 665 893 ENDIF -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90
r4153 r4292 30 30 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 31 31 USE in_out_manager ! 32 USE domvvl ! variable volume32 USE domvvl 33 33 34 34 IMPLICIT NONE … … 57 57 LOGICAL, INTENT( in ), OPTIONAL :: dyn3d_only ! T => only update baroclinic velocities 58 58 !! 59 INTEGER :: jk,ii,ij,ib,igrd ! Loop counter 60 LOGICAL :: ll_dyn2d, ll_dyn3d 61 !! 59 INTEGER :: jk,ii,ij,ib_bdy,ib,igrd ! Loop counter 60 LOGICAL :: ll_dyn2d, ll_dyn3d, ll_orlanski 61 !! 62 REAL(wp), POINTER, DIMENSION(:,:) :: phur1, phvr1 ! inverse depth at u and v points 62 63 63 64 IF( nn_timing == 1 ) CALL timing_start('bdy_dyn') … … 70 71 ENDIF 71 72 73 ll_orlanski = .false. 74 DO ib_bdy = 1, nb_bdy 75 IF ( cn_dyn2d(ib_bdy) == 'orlanski' .or. cn_dyn2d(ib_bdy) == 'orlanski_npo' & 76 & .or. cn_dyn3d(ib_bdy) == 'orlanski' .or. cn_dyn3d(ib_bdy) == 'orlanski_npo') ll_orlanski = .true. 77 ENDDO 78 72 79 !------------------------------------------------------- 73 80 ! Set pointers … … 77 84 phur => hur 78 85 phvr => hvr 79 CALL wrk_alloc(jpi,jpj,pu2d,pv2d) 86 CALL wrk_alloc(jpi,jpj,pua2d,pva2d) 87 IF ( ll_orlanski ) CALL wrk_alloc(jpi,jpj,pub2d,pvb2d,phur1,phvr1) 80 88 81 89 !------------------------------------------------------- … … 83 91 !------------------------------------------------------- 84 92 85 pu2d(:,:) = 0.e0 86 pv2d(:,:) = 0.e0 93 ! "After" velocities: 94 95 pua2d(:,:) = 0.e0 96 pva2d(:,:) = 0.e0 97 87 98 IF (lk_vvl) THEN 88 DO jk = 1, jpkm1 !! Vertically integrated momentum trends 89 pu2d(:,:) = pu2d(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 90 pv2d(:,:) = pv2d(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 91 END DO 92 pu2d(:,:) = pu2d(:,:) / ( hu_0(:,:) + sshu_a(:,:) + 1._wp - umask(:,:,1) ) 93 pv2d(:,:) = pv2d(:,:) / ( hv_0(:,:) + sshv_a(:,:) + 1._wp - vmask(:,:,1) ) 99 phur1(:,:) = 0. 100 phvr1(:,:) = 0. 101 DO jk = 1, jpkm1 102 phur1(:,:) = phur1(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) 103 phvr1(:,:) = phvr1(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) 104 pua2d(:,:) = pua2d(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 105 pva2d(:,:) = pva2d(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 106 END DO 107 phur1(:,:) = umask(:,:,1) / ( phur1(:,:) + 1. - umask(:,:,1) ) 108 phvr1(:,:) = vmask(:,:,1) / ( phvr1(:,:) + 1. - vmask(:,:,1) ) 109 pua2d(:,:) = pua2d(:,:) * phur1(:,:) 110 pva2d(:,:) = pva2d(:,:) * phvr1(:,:) 94 111 ELSE 95 DO jk = 1, jpkm1 !! Vertically integrated momentum trends96 pu 2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)97 pv 2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)98 END DO 99 pu 2d(:,:) = pu2d(:,:) * phur(:,:)100 pv 2d(:,:) = pv2d(:,:) * phvr(:,:)112 DO jk = 1, jpkm1 113 pua2d(:,:) = pua2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 114 pva2d(:,:) = pva2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 115 END DO 116 pua2d(:,:) = pua2d(:,:) * phur(:,:) 117 pva2d(:,:) = pva2d(:,:) * phvr(:,:) 101 118 ENDIF 119 102 120 DO jk = 1 , jpkm1 103 ua(:,:,jk) = ua(:,:,jk) - pu 2d(:,:) * umask(:,:,jk)104 va(:,:,jk) = va(:,:,jk) - pv 2d(:,:) * vmask(:,:,jk)121 ua(:,:,jk) = ua(:,:,jk) - pua2d(:,:) 122 va(:,:,jk) = va(:,:,jk) - pva2d(:,:) 105 123 END DO 124 125 ! "Before" velocities (required for Orlanski condition): 126 127 IF ( ll_orlanski ) THEN 128 pub2d(:,:) = 0.e0 129 pvb2d(:,:) = 0.e0 130 131 IF (lk_vvl) THEN 132 phur1(:,:) = 0. 133 phvr1(:,:) = 0. 134 DO jk = 1, jpkm1 !! Vertically integrated momentum trends 135 phur1(:,:) = phur1(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 136 phvr1(:,:) = phvr1(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 137 pub2d(:,:) = pub2d(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) * ub(:,:,jk) 138 pvb2d(:,:) = pvb2d(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) * vb(:,:,jk) 139 END DO 140 phur1(:,:) = umask(:,:,1) / ( phur1(:,:) + 1. - umask(:,:,1) ) 141 phvr1(:,:) = vmask(:,:,1) / ( phvr1(:,:) + 1. - vmask(:,:,1) ) 142 pub2d(:,:) = pub2d(:,:) * phur1(:,:) 143 pvb2d(:,:) = pvb2d(:,:) * phvr1(:,:) 144 ELSE 145 DO jk = 1, jpkm1 !! Vertically integrated momentum trends 146 pub2d(:,:) = pub2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ub(:,:,jk) 147 pvb2d(:,:) = pvb2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * vb(:,:,jk) 148 END DO 149 pub2d(:,:) = pub2d(:,:) * phur(:,:) 150 pvb2d(:,:) = pvb2d(:,:) * phvr(:,:) 151 ENDIF 152 153 DO jk = 1 , jpkm1 154 ub(:,:,jk) = ub(:,:,jk) - pub2d(:,:) 155 vb(:,:,jk) = vb(:,:,jk) - pvb2d(:,:) 156 END DO 157 END IF 106 158 107 159 !------------------------------------------------------- … … 119 171 120 172 DO jk = 1 , jpkm1 121 ua(:,:,jk) = ( ua(:,:,jk) + pu 2d(:,:) ) * umask(:,:,jk)122 va(:,:,jk) = ( va(:,:,jk) + pv 2d(:,:) ) * vmask(:,:,jk)173 ua(:,:,jk) = ( ua(:,:,jk) + pua2d(:,:) ) * umask(:,:,jk) 174 va(:,:,jk) = ( va(:,:,jk) + pva2d(:,:) ) * vmask(:,:,jk) 123 175 END DO 124 176 125 CALL wrk_dealloc(jpi,jpj,pu2d,pv2d) 177 IF ( ll_orlanski ) THEN 178 DO jk = 1 , jpkm1 179 ub(:,:,jk) = ( ub(:,:,jk) + pub2d(:,:) ) * umask(:,:,jk) 180 vb(:,:,jk) = ( vb(:,:,jk) + pvb2d(:,:) ) * vmask(:,:,jk) 181 END DO 182 END IF 183 184 CALL wrk_dealloc(jpi,jpj,pua2d,pva2d) 185 IF ( ll_orlanski ) CALL wrk_dealloc(jpi,jpj,pub2d,pvb2d) 126 186 127 187 IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn') -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90
r3680 r4292 6 6 !! History : 3.4 ! 2011 (D. Storkey) new module as part of BDY rewrite 7 7 !! 3.5 ! 2012 (S. Mocavero, I. Epicoco) Optimization of BDY communications 8 !! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes 8 9 !!---------------------------------------------------------------------- 9 10 #if defined key_bdy … … 11 12 !! 'key_bdy' : Unstructured Open Boundary Condition 12 13 !!---------------------------------------------------------------------- 13 !! bdy_dyn2d : Apply open boundary conditions to barotropic variables. 14 !! bdy_dyn2d_fla : Apply Flather condition 14 !! bdy_dyn2d : Apply open boundary conditions to barotropic variables. 15 !! bdy_dyn2d_frs : Apply Flow Relaxation Scheme 16 !! bdy_dyn2d_fla : Apply Flather condition 17 !! bdy_dyn2d_orlanski : Orlanski Radiation 18 !! bdy_ssh : Duplicate sea level across open boundaries 15 19 !!---------------------------------------------------------------------- 16 20 USE timing ! Timing … … 18 22 USE dom_oce ! ocean space and time domain 19 23 USE bdy_oce ! ocean open boundary conditions 24 USE bdylib ! BDY library routines 20 25 USE dynspg_oce ! for barotropic variables 21 26 USE phycst ! physical constants … … 26 31 PRIVATE 27 32 28 PUBLIC bdy_dyn2d ! routine called in dynspg_ts and bdy_dyn 33 PUBLIC bdy_dyn2d ! routine called in dynspg_ts and bdy_dyn 34 PUBLIC bdy_ssh ! routine called in dynspg_ts or sshwzv 29 35 30 36 !!---------------------------------------------------------------------- … … 48 54 DO ib_bdy=1, nb_bdy 49 55 50 SELECT CASE( nn_dyn2d(ib_bdy) )51 CASE( jp_none)56 SELECT CASE( cn_dyn2d(ib_bdy) ) 57 CASE('none') 52 58 CYCLE 53 CASE( jp_frs)59 CASE('frs') 54 60 CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy ) 55 CASE( jp_flather)61 CASE('flather') 56 62 CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy ) 63 CASE('orlanski') 64 CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. ) 65 CASE('orlanski_npo') 66 CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. ) 57 67 CASE DEFAULT 58 68 CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' ) … … 89 99 ij = idx%nbj(jb,igrd) 90 100 zwgt = idx%nbw(jb,igrd) 91 pu 2d(ii,ij) = ( pu2d(ii,ij) + zwgt * ( dta%u2d(jb) - pu2d(ii,ij) ) ) * umask(ii,ij,1)101 pua2d(ii,ij) = ( pua2d(ii,ij) + zwgt * ( dta%u2d(jb) - pua2d(ii,ij) ) ) * umask(ii,ij,1) 92 102 END DO 93 103 ! … … 97 107 ij = idx%nbj(jb,igrd) 98 108 zwgt = idx%nbw(jb,igrd) 99 pv 2d(ii,ij) = ( pv2d(ii,ij) + zwgt * ( dta%v2d(jb) - pv2d(ii,ij) ) ) * vmask(ii,ij,1)109 pva2d(ii,ij) = ( pva2d(ii,ij) + zwgt * ( dta%v2d(jb) - pva2d(ii,ij) ) ) * vmask(ii,ij,1) 100 110 END DO 101 CALL lbc_bdy_lnk( pu 2d, 'U', -1., ib_bdy )102 CALL lbc_bdy_lnk( pv 2d, 'V', -1., ib_bdy) ! Boundary points should be updated111 CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy ) 112 CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy) ! Boundary points should be updated 103 113 ! 104 114 IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_frs') … … 133 143 INTEGER :: jb, igrd ! dummy loop indices 134 144 INTEGER :: ii, ij, iim1, iip1, ijm1, ijp1 ! 2D addresses 145 REAL(wp), POINTER :: flagu, flagv ! short cuts 135 146 REAL(wp) :: zcorr ! Flather correction 136 147 REAL(wp) :: zforc ! temporary scalar 148 REAL(wp) :: zflag, z1_2 ! " " 137 149 !!---------------------------------------------------------------------- 138 150 139 151 IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_fla') 152 153 z1_2 = 0.5_wp 140 154 141 155 ! ---------------------------------! … … 160 174 ii = idx%nbi(jb,igrd) 161 175 ij = idx%nbj(jb,igrd) 162 iim1 = ii + MAX( 0, INT( idx%flagu(jb) ) ) ! T pts i-indice inside the boundary 163 iip1 = ii - MIN( 0, INT( idx%flagu(jb) ) ) ! T pts i-indice outside the boundary 176 flagu => idx%flagu(jb,igrd) 177 iim1 = ii + MAX( 0, INT( flagu ) ) ! T pts i-indice inside the boundary 178 iip1 = ii - MIN( 0, INT( flagu ) ) ! T pts i-indice outside the boundary 164 179 ! 165 zcorr = - idx%flagu(jb) * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 166 zforc = dta%u2d(jb) 167 pu2d(ii,ij) = zforc + zcorr * umask(ii,ij,1) 180 zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 181 182 ! jchanut tschanges: Set zflag to 0 below to revert to Flather scheme 183 ! Use characteristics method instead 184 zflag = ABS(flagu) 185 zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij) 186 pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 168 187 END DO 169 188 ! … … 173 192 ii = idx%nbi(jb,igrd) 174 193 ij = idx%nbj(jb,igrd) 175 ijm1 = ij + MAX( 0, INT( idx%flagv(jb) ) ) ! T pts j-indice inside the boundary 176 ijp1 = ij - MIN( 0, INT( idx%flagv(jb) ) ) ! T pts j-indice outside the boundary 194 flagv => idx%flagv(jb,igrd) 195 ijm1 = ij + MAX( 0, INT( flagv ) ) ! T pts j-indice inside the boundary 196 ijp1 = ij - MIN( 0, INT( flagv ) ) ! T pts j-indice outside the boundary 177 197 ! 178 zcorr = - idx%flagv(jb) * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 179 zforc = dta%v2d(jb) 180 pv2d(ii,ij) = zforc + zcorr * vmask(ii,ij,1) 181 END DO 182 CALL lbc_bdy_lnk( pu2d, 'U', -1., ib_bdy ) ! Boundary points should be updated 183 CALL lbc_bdy_lnk( pv2d, 'V', -1., ib_bdy ) ! 198 zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 199 200 ! jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme 201 ! Use characteristics method instead 202 zflag = ABS(flagv) 203 zforc = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1) 204 pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1) 205 END DO 206 CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy ) ! Boundary points should be updated 207 CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy ) ! 184 208 ! 185 209 IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_fla') 186 210 ! 187 211 END SUBROUTINE bdy_dyn2d_fla 212 213 214 SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, ll_npo ) 215 !!---------------------------------------------------------------------- 216 !! *** SUBROUTINE bdy_dyn2d_orlanski *** 217 !! 218 !! - Apply Orlanski radiation condition adaptively: 219 !! - radiation plus weak nudging at outflow points 220 !! - no radiation and strong nudging at inflow points 221 !! 222 !! 223 !! References: Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001) 224 !!---------------------------------------------------------------------- 225 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 226 TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data 227 INTEGER, INTENT(in) :: ib_bdy ! number of current open boundary set 228 LOGICAL, INTENT(in) :: ll_npo ! flag for NPO version 229 230 INTEGER :: ib, igrd ! dummy loop indices 231 INTEGER :: ii, ij, iibm1, ijbm1 ! indices 232 !!---------------------------------------------------------------------- 233 234 IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_orlanski') 235 ! 236 igrd = 2 ! Orlanski bc on u-velocity; 237 ! 238 CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, ll_npo ) 239 240 igrd = 3 ! Orlanski bc on v-velocity 241 ! 242 CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, ll_npo ) 243 ! 244 IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski') 245 ! 246 CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy ) ! Boundary points should be updated 247 CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy ) ! 248 ! 249 IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski') 250 ! 251 END SUBROUTINE bdy_dyn2d_orlanski 252 253 SUBROUTINE bdy_ssh( zssh ) 254 !!---------------------------------------------------------------------- 255 !! *** SUBROUTINE bdy_ssh *** 256 !! 257 !! ** Purpose : Duplicate sea level across open boundaries 258 !! 259 !!---------------------------------------------------------------------- 260 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zssh ! Sea level 261 !! 262 INTEGER :: ib_bdy, ib, igrd ! local integers 263 INTEGER :: ii, ij, zcoef, zcoef1, zcoef2, ip, jp ! " " 264 265 igrd = 1 ! Everything is at T-points here 266 267 DO ib_bdy = 1, nb_bdy 268 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 269 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 270 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 271 ! Set gradient direction: 272 zcoef1 = bdytmask(ii-1,ij ) + bdytmask(ii+1,ij ) 273 zcoef2 = bdytmask(ii ,ij-1) + bdytmask(ii ,ij+1) 274 IF ( zcoef1+zcoef2 == 0 ) THEN 275 ! corner 276 ! zcoef = tmask(ii-1,ij,1) + tmask(ii+1,ij,1) + tmask(ii,ij-1,1) + tmask(ii,ij+1,1) 277 ! zssh(ii,ij) = zssh(ii-1,ij ) * tmask(ii-1,ij ,1) + & 278 ! & zssh(ii+1,ij ) * tmask(ii+1,ij ,1) + & 279 ! & zssh(ii ,ij-1) * tmask(ii ,ij-1,1) + & 280 ! & zssh(ii ,ij+1) * tmask(ii ,ij+1,1) 281 zcoef = bdytmask(ii-1,ij) + bdytmask(ii+1,ij) + bdytmask(ii,ij-1) + bdytmask(ii,ij+1) 282 zssh(ii,ij) = zssh(ii-1,ij ) * bdytmask(ii-1,ij ) + & 283 & zssh(ii+1,ij ) * bdytmask(ii+1,ij ) + & 284 & zssh(ii ,ij-1) * bdytmask(ii ,ij-1) + & 285 & zssh(ii ,ij+1) * bdytmask(ii ,ij+1) 286 zssh(ii,ij) = ( zssh(ii,ij) / MAX( 1, zcoef) ) * tmask(ii,ij,1) 287 ELSE 288 ip = bdytmask(ii+1,ij ) - bdytmask(ii-1,ij ) 289 jp = bdytmask(ii ,ij+1) - bdytmask(ii ,ij-1) 290 zssh(ii,ij) = zssh(ii+ip,ij+jp) * tmask(ii+ip,ij+jp,1) 291 ENDIF 292 END DO 293 294 ! Boundary points should be updated 295 CALL lbc_bdy_lnk( zssh(:,:), 'T', 1., ib_bdy ) 296 END DO 297 298 END SUBROUTINE bdy_ssh 299 188 300 #else 189 301 !!---------------------------------------------------------------------- … … 192 304 CONTAINS 193 305 SUBROUTINE bdy_dyn2d( kt ) ! Empty routine 194 WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt 306 INTEGER, intent(in) :: kt 307 WRITE(*,*) 'bdy_dyn2d: You should not have seen this print! error?', kt 195 308 END SUBROUTINE bdy_dyn2d 309 196 310 #endif 197 311 198 312 !!====================================================================== 199 313 END MODULE bdydyn2d 314 -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90
r3703 r4292 19 19 USE dom_oce ! ocean space and time domain 20 20 USE bdy_oce ! ocean open boundary conditions 21 USE bdylib ! for orlanski library routines 21 22 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 22 23 USE in_out_manager ! … … 52 53 DO ib_bdy=1, nb_bdy 53 54 54 !!$ IF ( using Orlanski radiation conditions ) THEN 55 !!$ CALL bdy_rad( kt, bdyidx(ib_bdy) ) 56 !!$ ENDIF 57 58 SELECT CASE( nn_dyn3d(ib_bdy) ) 59 CASE(jp_none) 55 SELECT CASE( cn_dyn3d(ib_bdy) ) 56 CASE('none') 60 57 CYCLE 61 CASE( jp_frs)58 CASE('frs') 62 59 CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 63 CASE( 2)60 CASE('specified') 64 61 CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 65 CASE( 3)62 CASE('zero') 66 63 CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 64 CASE('orlanski') 65 CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. ) 66 CASE('orlanski_npo') 67 CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. ) 67 68 CASE DEFAULT 68 69 CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' ) … … 109 110 END DO 110 111 END DO 111 CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy ) ; CALL lbc_bdy_lnk( va, 'V', -1.,ib_bdy ) ! Boundary points should be updated 112 CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy ) ! Boundary points should be updated 113 CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy ) 112 114 ! 113 115 IF( kt .eq. nit000 ) CLOSE( unit = 102 ) … … 204 206 END DO 205 207 END DO 206 CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy ) ; CALL lbc_bdy_lnk( va, 'V', -1.,ib_bdy ) ! Boundary points should be updated 208 CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy ) ! Boundary points should be updated 209 CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy ) 207 210 ! 208 211 IF( kt .eq. nit000 ) CLOSE( unit = 102 ) … … 211 214 212 215 END SUBROUTINE bdy_dyn3d_frs 216 217 SUBROUTINE bdy_dyn3d_orlanski( idx, dta, ib_bdy, ll_npo ) 218 !!---------------------------------------------------------------------- 219 !! *** SUBROUTINE bdy_dyn3d_orlanski *** 220 !! 221 !! - Apply Orlanski radiation to baroclinic velocities. 222 !! - Wrapper routine for bdy_orlanski_3d 223 !! 224 !! 225 !! References: Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001) 226 !!---------------------------------------------------------------------- 227 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 228 TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data 229 INTEGER, INTENT(in) :: ib_bdy ! BDY set index 230 LOGICAL, INTENT(in) :: ll_npo ! switch for NPO version 231 232 INTEGER :: jb, igrd ! dummy loop indices 233 !!---------------------------------------------------------------------- 234 235 IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_orlanski') 236 ! 237 !! Note that at this stage the ub and ua arrays contain the baroclinic velocities. 238 ! 239 igrd = 2 ! Orlanski bc on u-velocity; 240 ! 241 CALL bdy_orlanski_3d( idx, igrd, ub, ua, dta%u3d, ll_npo ) 242 243 igrd = 3 ! Orlanski bc on v-velocity 244 ! 245 CALL bdy_orlanski_3d( idx, igrd, vb, va, dta%v3d, ll_npo ) 246 ! 247 CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy ) ! Boundary points should be updated 248 CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy ) 249 ! 250 IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_orlanski') 251 ! 252 END SUBROUTINE bdy_dyn3d_orlanski 253 213 254 214 255 SUBROUTINE bdy_dyn3d_dmp( kt ) … … 225 266 REAL(wp) :: zwgt ! boundary weight 226 267 INTEGER :: ib_bdy ! loop index 268 REAL(wp), POINTER, DIMENSION(:,:) :: phur1, phvr1 ! inverse depth at u and v points 227 269 !!---------------------------------------------------------------------- 228 270 ! … … 232 274 ! Remove barotropic part from before velocity 233 275 !------------------------------------------------------- 234 CALL wrk_alloc(jpi,jpj,pu2d,pv2d) 235 236 pu2d(:,:) = 0.e0 237 pv2d(:,:) = 0.e0 238 276 CALL wrk_alloc(jpi,jpj,pub2d,pvb2d,phur1,phvr1) 277 278 pub2d(:,:) = 0.e0 279 pvb2d(:,:) = 0.e0 280 281 phur1(:,:) = 0. 282 phvr1(:,:) = 0. 239 283 DO jk = 1, jpkm1 240 284 #if defined key_vvl 241 pu2d(:,:) = pu2d(:,:) + fse3u_b(:,:,jk)* ub(:,:,jk) *umask(:,:,jk) 242 pv2d(:,:) = pv2d(:,:) + fse3v_b(:,:,jk)* vb(:,:,jk) *vmask(:,:,jk) 285 phur1(:,:) = phur1(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) 286 phvr1(:,:) = phvr1(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) 287 pub2d(:,:) = pub2d(:,:) + fse3u_b(:,:,jk)* ub(:,:,jk) *umask(:,:,jk) 288 pvb2d(:,:) = pvb2d(:,:) + fse3v_b(:,:,jk)* vb(:,:,jk) *vmask(:,:,jk) 243 289 #else 244 pu 2d(:,:) = pu2d(:,:) + fse3u_0(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)245 pv 2d(:,:) = pv2d(:,:) + fse3v_0(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)290 pub2d(:,:) = pub2d(:,:) + fse3u_0(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 291 pvb2d(:,:) = pvb2d(:,:) + fse3v_0(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 246 292 #endif 247 293 END DO 248 294 249 295 IF( lk_vvl ) THEN 250 pu2d(:,:) = pu2d(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 251 pv2d(:,:) = pv2d(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 296 phur1(:,:) = umask(:,:,1) / ( phur1(:,:) + 1. - umask(:,:,1) ) 297 phvr1(:,:) = vmask(:,:,1) / ( phvr1(:,:) + 1. - vmask(:,:,1) ) 298 pub2d(:,:) = pub2d(:,:) * umask(:,:,1) * phur1(:,:) 299 pvb2d(:,:) = pvb2d(:,:) * vmask(:,:,1) * phvr1(:,:) 252 300 ELSE 253 pu 2d(:,:) = pv2d(:,:) * hur(:,:)254 pv 2d(:,:) = pu2d(:,:) * hvr(:,:)301 pub2d(:,:) = pvb2d(:,:) * hur(:,:) 302 pvb2d(:,:) = pub2d(:,:) * hvr(:,:) 255 303 ENDIF 256 304 257 305 DO ib_bdy=1, nb_bdy 258 IF ( ln_dyn3d_dmp(ib_bdy) .and.nn_dyn3d(ib_bdy).gt.0) THEN306 IF ( ln_dyn3d_dmp(ib_bdy) .and. cn_dyn3d(ib_bdy) /= 'none' ) THEN 259 307 igrd = 2 ! Relaxation of zonal velocity 260 308 DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd) … … 264 312 DO jk = 1, jpkm1 265 313 ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - & 266 ub(ii,ij,jk) + pu 2d(ii,ij)) ) * umask(ii,ij,jk)314 ub(ii,ij,jk) + pub2d(ii,ij)) ) * umask(ii,ij,jk) 267 315 END DO 268 316 END DO … … 275 323 DO jk = 1, jpkm1 276 324 va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) - & 277 vb(ii,ij,jk) + pv 2d(ii,ij)) ) * vmask(ii,ij,jk)325 vb(ii,ij,jk) + pvb2d(ii,ij)) ) * vmask(ii,ij,jk) 278 326 END DO 279 327 END DO … … 281 329 ENDDO 282 330 ! 283 CALL wrk_dealloc(jpi,jpj,pu 2d,pv2d)331 CALL wrk_dealloc(jpi,jpj,pub2d,pvb2d) 284 332 ! 285 333 CALL lbc_lnk( ua, 'U', -1. ) ; CALL lbc_lnk( va, 'V', -1. ) ! Boundary points should be updated -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
r4148 r4292 21 21 !! bdy_init : Initialization of unstructured open boundaries 22 22 !!---------------------------------------------------------------------- 23 USE wrk_nemo ! Memory Allocation 23 24 USE timing ! Timing 24 25 USE oce ! ocean dynamics and tracers variables … … 79 80 INTEGER :: jpbdtau, jpbdtas ! - - 80 81 INTEGER :: ib_bdy1, ib_bdy2, ib1, ib2 ! - - 82 INTEGER :: i_offset, j_offset ! - - 81 83 INTEGER, POINTER :: nbi, nbj, nbr ! short cuts 82 REAL , POINTER :: flagu, flagv ! - - 84 REAL(wp), POINTER :: flagu, flagv ! - - 85 REAL(wp), POINTER, DIMENSION(:,:) :: pmask ! pointer to 2D mask fields 83 86 REAL(wp) :: zefl, zwfl, znfl, zsfl ! local scalars 84 87 INTEGER, DIMENSION (2) :: kdimsz … … 90 93 INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b ! Flags for boundaries receiving 91 94 INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4) ! Arrays for neighbours coordinates 95 REAL(wp), POINTER, DIMENSION(:,:) :: zfmask ! temporary fmask array excluding coastal boundary condition (shlat) 92 96 93 97 !! 94 NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file, &95 & ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta,&96 & nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta,&97 & ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, 98 #if defined key_lim299 & nn_ice_lim2, nn_ice_lim2_dta,&98 NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file, & 99 & ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta, & 100 & cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta, & 101 & ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, & 102 #if ( defined key_lim2 || defined key_lim3 ) 103 & cn_ice_lim, nn_ice_lim_dta, & 100 104 #endif 101 105 & ln_vol, nn_volctl, nn_rimwidth … … 156 160 157 161 IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution: ' 158 SELECT CASE( nn_dyn2d(ib_bdy) ) 159 CASE(jp_none) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 160 CASE(jp_frs) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 161 CASE(jp_flather) ; IF(lwp) WRITE(numout,*) ' Flather radiation condition' 162 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_dyn2d' ) 162 SELECT CASE( cn_dyn2d(ib_bdy) ) 163 CASE('none') 164 IF(lwp) WRITE(numout,*) ' no open boundary condition' 165 dta_bdy(ib_bdy)%ll_ssh = .false. 166 dta_bdy(ib_bdy)%ll_u2d = .false. 167 dta_bdy(ib_bdy)%ll_v2d = .false. 168 CASE('frs') 169 IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 170 dta_bdy(ib_bdy)%ll_ssh = .false. 171 dta_bdy(ib_bdy)%ll_u2d = .true. 172 dta_bdy(ib_bdy)%ll_v2d = .true. 173 CASE('flather') 174 IF(lwp) WRITE(numout,*) ' Flather radiation condition' 175 dta_bdy(ib_bdy)%ll_ssh = .true. 176 dta_bdy(ib_bdy)%ll_u2d = .true. 177 dta_bdy(ib_bdy)%ll_v2d = .true. 178 CASE('orlanski') 179 IF(lwp) WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging' 180 dta_bdy(ib_bdy)%ll_ssh = .false. 181 dta_bdy(ib_bdy)%ll_u2d = .true. 182 dta_bdy(ib_bdy)%ll_v2d = .true. 183 CASE('orlanski_npo') 184 IF(lwp) WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging' 185 dta_bdy(ib_bdy)%ll_ssh = .false. 186 dta_bdy(ib_bdy)%ll_u2d = .true. 187 dta_bdy(ib_bdy)%ll_v2d = .true. 188 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_dyn2d' ) 163 189 END SELECT 164 IF( nn_dyn2d(ib_bdy) .gt. 0) THEN190 IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN 165 191 SELECT CASE( nn_dyn2d_dta(ib_bdy) ) ! 166 192 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data' … … 177 203 178 204 IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities: ' 179 SELECT CASE( nn_dyn3d(ib_bdy) ) 180 CASE(jp_none) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 181 CASE(jp_frs) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 182 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Specified value' 183 CASE( 3 ) ; IF(lwp) WRITE(numout,*) ' Zero baroclinic velocities (runoff case)' 184 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_dyn3d' ) 205 SELECT CASE( cn_dyn3d(ib_bdy) ) 206 CASE('none') 207 IF(lwp) WRITE(numout,*) ' no open boundary condition' 208 dta_bdy(ib_bdy)%ll_u3d = .false. 209 dta_bdy(ib_bdy)%ll_v3d = .false. 210 CASE('frs') 211 IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 212 dta_bdy(ib_bdy)%ll_u3d = .true. 213 dta_bdy(ib_bdy)%ll_v3d = .true. 214 CASE('specified') 215 IF(lwp) WRITE(numout,*) ' Specified value' 216 dta_bdy(ib_bdy)%ll_u3d = .true. 217 dta_bdy(ib_bdy)%ll_v3d = .true. 218 CASE('zero') 219 IF(lwp) WRITE(numout,*) ' Zero baroclinic velocities (runoff case)' 220 dta_bdy(ib_bdy)%ll_u3d = .false. 221 dta_bdy(ib_bdy)%ll_v3d = .false. 222 CASE('orlanski') 223 IF(lwp) WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging' 224 dta_bdy(ib_bdy)%ll_u3d = .true. 225 dta_bdy(ib_bdy)%ll_v3d = .true. 226 CASE('orlanski_npo') 227 IF(lwp) WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging' 228 dta_bdy(ib_bdy)%ll_u3d = .true. 229 dta_bdy(ib_bdy)%ll_v3d = .true. 230 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_dyn3d' ) 185 231 END SELECT 186 IF( nn_dyn3d(ib_bdy) .gt. 0) THEN232 IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 187 233 SELECT CASE( nn_dyn3d_dta(ib_bdy) ) ! 188 234 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data' … … 193 239 194 240 IF ( ln_dyn3d_dmp(ib_bdy) ) THEN 195 IF ( nn_dyn3d(ib_bdy).EQ.0) THEN241 IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN 196 242 IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.' 197 243 ln_dyn3d_dmp(ib_bdy)=.false. 198 ELSEIF ( nn_dyn3d(ib_bdy).EQ.1) THEN244 ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN 199 245 CALL ctl_stop( 'Use FRS OR relaxation' ) 200 246 ELSE … … 202 248 IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days' 203 249 IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 250 dta_bdy(ib_bdy)%ll_u3d = .true. 251 dta_bdy(ib_bdy)%ll_v3d = .true. 204 252 ENDIF 205 253 ELSE … … 209 257 210 258 IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity: ' 211 SELECT CASE( nn_tra(ib_bdy) ) 212 CASE(jp_none) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 213 CASE(jp_frs) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 214 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Specified value' 215 CASE( 3 ) ; IF(lwp) WRITE(numout,*) ' Neumann conditions' 216 CASE( 4 ) ; IF(lwp) WRITE(numout,*) ' Runoff conditions : Neumann for T and specified to 0.1 for salinity' 217 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_tra' ) 259 SELECT CASE( cn_tra(ib_bdy) ) 260 CASE('none') 261 IF(lwp) WRITE(numout,*) ' no open boundary condition' 262 dta_bdy(ib_bdy)%ll_tem = .false. 263 dta_bdy(ib_bdy)%ll_sal = .false. 264 CASE('frs') 265 IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 266 dta_bdy(ib_bdy)%ll_tem = .true. 267 dta_bdy(ib_bdy)%ll_sal = .true. 268 CASE('specified') 269 IF(lwp) WRITE(numout,*) ' Specified value' 270 dta_bdy(ib_bdy)%ll_tem = .true. 271 dta_bdy(ib_bdy)%ll_sal = .true. 272 CASE('neumann') 273 IF(lwp) WRITE(numout,*) ' Neumann conditions' 274 dta_bdy(ib_bdy)%ll_tem = .false. 275 dta_bdy(ib_bdy)%ll_sal = .false. 276 CASE('runoff') 277 IF(lwp) WRITE(numout,*) ' Runoff conditions : Neumann for T and specified to 0.1 for salinity' 278 dta_bdy(ib_bdy)%ll_tem = .false. 279 dta_bdy(ib_bdy)%ll_sal = .false. 280 CASE('orlanski') 281 IF(lwp) WRITE(numout,*) ' Orlanski (fully oblique) radiation condition with adaptive nudging' 282 dta_bdy(ib_bdy)%ll_tem = .true. 283 dta_bdy(ib_bdy)%ll_sal = .true. 284 CASE('orlanski_npo') 285 IF(lwp) WRITE(numout,*) ' Orlanski (NPO) radiation condition with adaptive nudging' 286 dta_bdy(ib_bdy)%ll_tem = .true. 287 dta_bdy(ib_bdy)%ll_sal = .true. 288 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_tra' ) 218 289 END SELECT 219 IF( nn_tra(ib_bdy) .gt. 0) THEN290 IF( cn_tra(ib_bdy) /= 'none' ) THEN 220 291 SELECT CASE( nn_tra_dta(ib_bdy) ) ! 221 292 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data' … … 226 297 227 298 IF ( ln_tra_dmp(ib_bdy) ) THEN 228 IF ( nn_tra(ib_bdy).EQ.0) THEN299 IF ( cn_tra(ib_bdy) == 'none' ) THEN 229 300 IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.' 230 301 ln_tra_dmp(ib_bdy)=.false. 231 ELSEIF ( nn_tra(ib_bdy).EQ.1) THEN302 ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN 232 303 CALL ctl_stop( 'Use FRS OR relaxation' ) 233 304 ELSE 234 305 IF(lwp) WRITE(numout,*) ' + T/S relaxation zone' 235 306 IF(lwp) WRITE(numout,*) ' Damping time scale: ',rn_time_dmp(ib_bdy),' days' 307 IF(lwp) WRITE(numout,*) ' Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days' 236 308 IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 309 dta_bdy(ib_bdy)%ll_tem = .true. 310 dta_bdy(ib_bdy)%ll_sal = .true. 237 311 ENDIF 238 312 ELSE … … 243 317 #if defined key_lim2 244 318 IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice: ' 245 SELECT CASE( nn_ice_lim2(ib_bdy) ) 246 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 247 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 248 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_tra' ) 319 SELECT CASE( cn_ice_lim(ib_bdy) ) 320 CASE('none') 321 IF(lwp) WRITE(numout,*) ' no open boundary condition' 322 dta_bdy(ib_bdy)%ll_frld = .false. 323 dta_bdy(ib_bdy)%ll_hicif = .false. 324 dta_bdy(ib_bdy)%ll_hsnif = .false. 325 CASE('frs') 326 IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 327 dta_bdy(ib_bdy)%ll_frld = .true. 328 dta_bdy(ib_bdy)%ll_hicif = .true. 329 dta_bdy(ib_bdy)%ll_hsnif = .true. 330 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_ice_lim' ) 249 331 END SELECT 250 IF( nn_ice_lim2(ib_bdy) .gt. 0) THEN251 SELECT CASE( nn_ice_lim 2_dta(ib_bdy) ) !332 IF( cn_ice_lim(ib_bdy) /= 'none' ) THEN 333 SELECT CASE( nn_ice_lim_dta(ib_bdy) ) ! 252 334 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data' 253 335 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file' 254 CASE DEFAULT ; CALL ctl_stop( 'nn_ice_lim2_dta must be 0 or 1' ) 336 CASE DEFAULT ; CALL ctl_stop( 'nn_ice_lim_dta must be 0 or 1' ) 337 END SELECT 338 ENDIF 339 IF(lwp) WRITE(numout,*) 340 #elif defined key_lim3 341 IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice: ' 342 SELECT CASE( cn_ice_lim(ib_bdy) ) 343 CASE('none') 344 IF(lwp) WRITE(numout,*) ' no open boundary condition' 345 dta_bdy(ib_bdy)%ll_a_i = .false. 346 dta_bdy(ib_bdy)%ll_ht_i = .false. 347 dta_bdy(ib_bdy)%ll_ht_s = .false. 348 CASE('frs') 349 IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 350 dta_bdy(ib_bdy)%ll_a_i = .true. 351 dta_bdy(ib_bdy)%ll_ht_i = .true. 352 dta_bdy(ib_bdy)%ll_ht_s = .true. 353 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for cn_ice_lim' ) 354 END SELECT 355 IF( cn_ice_lim(ib_bdy) /= 'none' ) THEN 356 SELECT CASE( nn_ice_lim_dta(ib_bdy) ) ! 357 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' initial state used for bdy data' 358 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' boundary data taken from file' 359 CASE DEFAULT ; CALL ctl_stop( 'nn_ice_lim_dta must be 0 or 1' ) 255 360 END SELECT 256 361 ENDIF … … 740 845 IF(lwp) THEN ! Since all procs read global data only need to do this check on one proc... 741 846 IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN 742 CALL ctl_stop('bdy_init : ERROR : boundary data in file & 743 must be defined in order of distance from edge nbr.', & 744 'A utility for re-ordering boundary coordinates and data & 745 files exists in the TOOLS/OBC directory') 847 CALL ctl_stop('bdy_init : ERROR : boundary data in file must be defined in order of distance from edge nbr.', & 848 'A utility for re-ordering boundary coordinates and data files exists in the TOOLS/OBC directory') 746 849 ENDIF 747 850 ENDIF … … 766 869 ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) ) 767 870 ALLOCATE( idx_bdy(ib_bdy)%nbd(ilen1,jpbgrd) ) 871 ALLOCATE( idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) ) 768 872 ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) ) 769 873 ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) ) 770 ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1 ) )771 ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1 ) )874 ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1,jpbgrd) ) 875 ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1,jpbgrd) ) 772 876 773 877 ! Dispatch mapping indices and discrete distances on each processor … … 937 1041 ENDDO 938 1042 ENDDO 1043 939 1044 ! definition of the i- and j- direction local boundaries arrays 940 1045 ! used for sending the boudaries … … 990 1095 nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 991 1096 idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) & 1097 & *(FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2. ! quadratic 1098 idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) & 992 1099 & *(FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2. ! quadratic 993 1100 END DO … … 1092 1199 ENDDO 1093 1200 1201 ! For the flagu/flagv calculation below we require a version of fmask without 1202 ! the land boundary condition (shlat) included: 1203 CALL wrk_alloc(jpi,jpj,zfmask) 1204 DO ij = 2, jpjm1 1205 DO ii = 2, jpim1 1206 zfmask(ii,ij) = tmask(ii,ij ,1) * tmask(ii+1,ij ,1) & 1207 & * tmask(ii,ij+1,1) * tmask(ii+1,ij+1,1) 1208 END DO 1209 END DO 1210 1094 1211 ! Lateral boundary conditions 1212 CALL lbc_lnk( zfmask , 'F', 1. ) 1095 1213 CALL lbc_lnk( fmask , 'F', 1. ) ; CALL lbc_lnk( bdytmask(:,:), 'T', 1. ) 1096 1214 CALL lbc_lnk( bdyumask(:,:), 'U', 1. ) ; CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) … … 1098 1216 DO ib_bdy = 1, nb_bdy ! Indices and directions of rim velocity components 1099 1217 1100 idx_bdy(ib_bdy)%flagu(: ) = 0.e01101 idx_bdy(ib_bdy)%flagv(: ) = 0.e01218 idx_bdy(ib_bdy)%flagu(:,:) = 0.e0 1219 idx_bdy(ib_bdy)%flagv(:,:) = 0.e0 1102 1220 icount = 0 1103 1221 1104 !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward 1105 !flagu = 0 : u is tangential 1106 !flagu = 1 : u is normal to the boundary and is direction is inward 1222 ! Calculate relationship of U direction to the local orientation of the boundary 1223 ! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward 1224 ! flagu = 0 : u is tangential 1225 ! flagu = 1 : u is normal to the boundary and is direction is inward 1107 1226 1108 igrd = 2 ! u-component 1109 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 1110 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 1111 nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 1112 zefl = bdytmask(nbi ,nbj) 1113 zwfl = bdytmask(nbi+1,nbj) 1114 IF( zefl + zwfl == 2 ) THEN 1115 icount = icount + 1 1116 ELSE 1117 idx_bdy(ib_bdy)%flagu(ib)=-zefl+zwfl 1118 ENDIF 1227 DO igrd = 1,jpbgrd 1228 SELECT CASE( igrd ) 1229 CASE( 1 ) 1230 pmask => umask(:,:,1) 1231 i_offset = 0 1232 CASE( 2 ) 1233 pmask => bdytmask 1234 i_offset = 1 1235 CASE( 3 ) 1236 pmask => zfmask(:,:) 1237 i_offset = 0 1238 END SELECT 1239 icount = 0 1240 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 1241 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 1242 nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 1243 zefl = pmask(nbi+i_offset-1,nbj) 1244 zwfl = pmask(nbi+i_offset,nbj) 1245 ! This error check only works if you are using the bdyXmask arrays 1246 IF( i_offset == 1 .and. zefl + zwfl == 2 ) THEN 1247 icount = icount + 1 1248 IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj) 1249 ELSE 1250 idx_bdy(ib_bdy)%flagu(ib,igrd) = -zefl + zwfl 1251 ENDIF 1252 END DO 1253 IF( icount /= 0 ) THEN 1254 IF(lwp) WRITE(numout,*) 1255 IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid(igrd),' grid points,', & 1256 ' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy 1257 IF(lwp) WRITE(numout,*) ' ========== ' 1258 IF(lwp) WRITE(numout,*) 1259 nstop = nstop + 1 1260 ENDIF 1119 1261 END DO 1120 1262 1121 !flagv = -1 : u component is normal to the dynamical boundary but its direction is outward 1122 !flagv = 0 : u is tangential 1123 !flagv = 1 : u is normal to the boundary and is direction is inward 1124 1125 igrd = 3 ! v-component 1126 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 1127 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 1128 nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 1129 znfl = bdytmask(nbi,nbj ) 1130 zsfl = bdytmask(nbi,nbj+1) 1131 IF( znfl + zsfl == 2 ) THEN 1132 icount = icount + 1 1133 ELSE 1134 idx_bdy(ib_bdy)%flagv(ib) = -znfl + zsfl 1135 END IF 1263 ! Calculate relationship of V direction to the local orientation of the boundary 1264 ! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward 1265 ! flagv = 0 : v is tangential 1266 ! flagv = 1 : v is normal to the boundary and is direction is inward 1267 1268 DO igrd = 1,jpbgrd 1269 SELECT CASE( igrd ) 1270 CASE( 1 ) 1271 pmask => vmask(:,:,1) 1272 j_offset = 0 1273 CASE( 2 ) 1274 pmask => zfmask(:,:) 1275 j_offset = 0 1276 CASE( 3 ) 1277 pmask => bdytmask 1278 j_offset = 1 1279 END SELECT 1280 icount = 0 1281 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 1282 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 1283 nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 1284 znfl = pmask(nbi,nbj+j_offset-1 ) 1285 zsfl = pmask(nbi,nbj+j_offset) 1286 ! This error check only works if you are using the bdyXmask arrays 1287 IF( j_offset == 1 .and. znfl + zsfl == 2 ) THEN 1288 IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(nbi),mjg(nbj) 1289 icount = icount + 1 1290 ELSE 1291 idx_bdy(ib_bdy)%flagv(ib,igrd) = -znfl + zsfl 1292 END IF 1293 END DO 1294 IF( icount /= 0 ) THEN 1295 IF(lwp) WRITE(numout,*) 1296 IF(lwp) WRITE(numout,*) ' E R R O R : Some ',cgrid(igrd),' grid points,', & 1297 ' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy 1298 IF(lwp) WRITE(numout,*) ' ========== ' 1299 IF(lwp) WRITE(numout,*) 1300 nstop = nstop + 1 1301 ENDIF 1136 1302 END DO 1137 1303 1138 IF( icount /= 0 ) THEN 1139 IF(lwp) WRITE(numout,*) 1140 IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,', & 1141 ' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_bdy 1142 IF(lwp) WRITE(numout,*) ' ========== ' 1143 IF(lwp) WRITE(numout,*) 1144 nstop = nstop + 1 1145 ENDIF 1146 1147 ENDDO 1304 END DO 1148 1305 1149 1306 ! Compute total lateral surface for volume correction: … … 1157 1314 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 1158 1315 nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 1159 flagu => idx_bdy(ib_bdy)%flagu(ib )1316 flagu => idx_bdy(ib_bdy)%flagu(ib,igrd) 1160 1317 bdysurftot = bdysurftot + hu (nbi , nbj) & 1161 1318 & * e2u (nbi , nbj) * ABS( flagu ) & … … 1170 1327 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 1171 1328 nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 1172 flagv => idx_bdy(ib_bdy)%flagv(ib )1329 flagv => idx_bdy(ib_bdy)%flagv(ib,igrd) 1173 1330 bdysurftot = bdysurftot + hv (nbi, nbj ) & 1174 1331 & * e1v (nbi, nbj ) * ABS( flagv ) & … … 1186 1343 DEALLOCATE(nbidta, nbjdta, nbrdta) 1187 1344 ENDIF 1345 1346 CALL wrk_dealloc(jpi,jpj,zfmask) 1188 1347 1189 1348 IF( nn_timing == 1 ) CALL timing_stop('bdy_init') … … 1580 1739 itest = 0 1581 1740 1582 IF ( nn_dyn2d(ib1)/=nn_dyn2d(ib2)) itest = itest + 11583 IF ( nn_dyn3d(ib1)/=nn_dyn3d(ib2)) itest = itest + 11584 IF ( nn_tra(ib1)/=nn_tra(ib2)) itest = itest + 11741 IF (cn_dyn2d(ib1)/=cn_dyn2d(ib2)) itest = itest + 1 1742 IF (cn_dyn3d(ib1)/=cn_dyn3d(ib2)) itest = itest + 1 1743 IF (cn_tra(ib1)/=cn_tra(ib2)) itest = itest + 1 1585 1744 ! 1586 1745 IF (nn_dyn2d_dta(ib1)/=nn_dyn2d_dta(ib2)) itest = itest + 1 -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
r4147 r4292 9 9 !! 3.3 ! 2010-09 (D.Storkey and E.O'Dea) bug fixes 10 10 !! 3.4 ! 2012-09 (G. Reffray and J. Chanut) New inputs + mods 11 !! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes 11 12 !!---------------------------------------------------------------------- 12 13 #if defined key_bdy … … 32 33 ! USE tide_mod ! Useless ?? 33 34 USE fldread, ONLY: fld_map 35 USE dynspg_oce, ONLY: lk_dynspg_ts 34 36 35 37 IMPLICIT NONE … … 38 40 PUBLIC bdytide_init ! routine called in bdy_init 39 41 PUBLIC bdytide_update ! routine called in bdy_dta 42 PUBLIC bdy_dta_tides ! routine called in dyn_spg_ts 40 43 41 44 TYPE, PUBLIC :: TIDES_DATA !: Storage for external tidal harmonics data … … 49 52 50 53 TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides !: External tidal harmonics data 54 TYPE(OBC_DATA) , PRIVATE, DIMENSION(jp_bdy) :: dta_bdy_s !: bdy external data (slow component) 51 55 52 56 !!---------------------------------------------------------------------- … … 131 135 ! JC: If FRS scheme is used, we assume that tidal is needed over the whole 132 136 ! relaxation area 133 IF( nn_dyn2d(ib_bdy) .eq. jp_frs) THEN137 IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN 134 138 ilen0(:)=nblen(:) 135 139 ELSE … … 146 150 ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) ) 147 151 148 td%ssh0(:,:,:) = 0. e0149 td%ssh (:,:,:) = 0.e0150 td%u0 (:,:,:) = 0.e0151 td%u (:,:,:) = 0.e0152 td%v0 (:,:,:) = 0.e0153 td%v (:,:,:) = 0.e0152 td%ssh0(:,:,:) = 0._wp 153 td%ssh (:,:,:) = 0._wp 154 td%u0 (:,:,:) = 0._wp 155 td%u (:,:,:) = 0._wp 156 td%v0 (:,:,:) = 0._wp 157 td%v (:,:,:) = 0._wp 154 158 155 159 IF (ln_bdytide_2ddta) THEN … … 255 259 ENDIF 256 260 ! 261 IF ( lk_dynspg_ts ) THEN ! Allocate arrays to save slowly varying boundary data during 262 ! time splitting integration 263 ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) ) 264 ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) ) 265 ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) ) 266 dta_bdy_s(ib_bdy)%ssh(:) = 0.e0 267 dta_bdy_s(ib_bdy)%u2d(:) = 0.e0 268 dta_bdy_s(ib_bdy)%v2d(:) = 0.e0 269 ENDIF 270 ! 257 271 ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2 258 272 ! … … 300 314 ENDIF 301 315 302 IF ( nsec_day == NINT(0.5 * rdttra(1)) .AND. zflag==1 ) THEN316 IF ( nsec_day == NINT(0.5_wp * rdttra(1)) .AND. zflag==1 ) THEN 303 317 ! 304 318 kt_tide = kt … … 321 335 322 336 IF( PRESENT(jit) ) THEN 323 z_arg = ( ((kt-kt_tide)-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) )337 z_arg = ((kt-kt_tide) * rdt + (jit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) ) 324 338 ELSE 325 339 z_arg = ((kt-kt_tide)+time_add) * rdt … … 327 341 328 342 ! Linear ramp on tidal component at open boundaries 329 zramp = 1. 330 IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0. ),1.)343 zramp = 1._wp 344 IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0._wp),1._wp) 331 345 332 346 DO itide = 1, nb_harmo … … 354 368 ! 355 369 END SUBROUTINE bdytide_update 370 371 SUBROUTINE bdy_dta_tides( kt, kit, time_offset ) 372 !!---------------------------------------------------------------------- 373 !! *** SUBROUTINE bdy_dta_tides *** 374 !! 375 !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays. 376 !! 377 !!---------------------------------------------------------------------- 378 INTEGER, INTENT( in ) :: kt ! Main timestep counter 379 INTEGER, INTENT( in ),OPTIONAL :: kit ! Barotropic timestep counter (for timesplitting option) 380 INTEGER, INTENT( in ),OPTIONAL :: time_offset ! time offset in units of timesteps. NB. if kit 381 ! is present then units = subcycle timesteps. 382 ! time_offset = 0 => get data at "now" time level 383 ! time_offset = -1 => get data at "before" time level 384 ! time_offset = +1 => get data at "after" time level 385 ! etc. 386 !! 387 LOGICAL :: lk_first_btstp ! =.TRUE. if time splitting and first barotropic step 388 INTEGER, DIMENSION(jpbgrd) :: ilen0 389 INTEGER, DIMENSION(1:jpbgrd) :: nblen, nblenrim ! short cuts 390 INTEGER :: itide, ib_bdy, ib, igrd ! loop indices 391 INTEGER :: time_add ! time offset in units of timesteps 392 REAL(wp) :: z_arg, z_sarg, zramp, zoff, z_cost, z_sist 393 !!---------------------------------------------------------------------- 394 395 IF( nn_timing == 1 ) CALL timing_start('bdy_dta_tides') 396 397 lk_first_btstp=.TRUE. 398 IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF 399 400 time_add = 0 401 IF( PRESENT(time_offset) ) THEN 402 time_add = time_offset 403 ENDIF 404 405 ! Absolute time from model initialization: 406 IF( PRESENT(kit) ) THEN 407 z_arg = ( kt + (kit+0.5_wp*(time_add-1)) / REAL(nn_baro,wp) ) * rdt 408 ELSE 409 z_arg = ( kt + time_add ) * rdt 410 ENDIF 411 412 ! Linear ramp on tidal component at open boundaries 413 zramp = 1. 414 IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.) 415 416 DO ib_bdy = 1,nb_bdy 417 418 ! line below should be simplified (runoff case) 419 !! CHANUT: TO BE SORTED OUT 420 !! IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(nn_tra(ib_bdy).NE.4)) THEN 421 IF ( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN 422 423 nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd) 424 nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd) 425 426 IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN 427 ilen0(:)=nblen(:) 428 ELSE 429 ilen0(:)=nblenrim(:) 430 ENDIF 431 432 ! We refresh nodal factors every day below 433 ! This should be done somewhere else 434 IF ( nsec_day == NINT(0.5_wp * rdttra(1)) .AND. lk_first_btstp ) THEN 435 ! 436 kt_tide = kt 437 ! 438 IF(lwp) THEN 439 WRITE(numout,*) 440 WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt 441 WRITE(numout,*) '~~~~~~~~~~~~~~ ' 442 ENDIF 443 ! 444 CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) ) 445 CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) ) 446 ! 447 ENDIF 448 zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time 449 ! 450 ! If time splitting, save data at first barotropic iteration 451 IF ( PRESENT(kit) ) THEN 452 IF ( lk_first_btstp ) THEN ! Save slow varying open boundary data: 453 dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy(ib_bdy)%ssh(1:ilen0(1)) 454 dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy(ib_bdy)%u2d(1:ilen0(2)) 455 dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy(ib_bdy)%v2d(1:ilen0(3)) 456 457 ELSE ! Initialize arrays from slow varying open boundary data: 458 dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) 459 dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) 460 dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) 461 ENDIF 462 ENDIF 463 ! 464 ! Update open boundary data arrays: 465 DO itide = 1, nb_harmo 466 ! 467 z_sarg = (z_arg + zoff) * omega_tide(itide) 468 z_cost = zramp * COS( z_sarg ) 469 z_sist = zramp * SIN( z_sarg ) 470 ! 471 igrd=1 ! SSH on tracer grid 472 DO ib = 1, ilen0(igrd) 473 dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + & 474 & ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + & 475 & tides(ib_bdy)%ssh(ib,itide,2)*z_sist ) 476 END DO 477 ! 478 igrd=2 ! U grid 479 DO ib = 1, ilen0(igrd) 480 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + & 481 & ( tides(ib_bdy)%u(ib,itide,1)*z_cost + & 482 & tides(ib_bdy)%u(ib,itide,2)*z_sist ) 483 END DO 484 ! 485 igrd=3 ! V grid 486 DO ib = 1, ilen0(igrd) 487 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + & 488 & ( tides(ib_bdy)%v(ib,itide,1)*z_cost + & 489 & tides(ib_bdy)%v(ib,itide,2)*z_sist ) 490 END DO 491 END DO 492 END IF 493 END DO 494 ! 495 IF( nn_timing == 1 ) CALL timing_stop('bdy_dta_tides') 496 ! 497 END SUBROUTINE bdy_dta_tides 356 498 357 499 SUBROUTINE tide_init_elevation( idx, td ) … … 460 602 WRITE(*,*) 'bdytide_update: You should not have seen this print! error?', kt, jit 461 603 END SUBROUTINE bdytide_update 604 SUBROUTINE bdy_dta_tides( kt, kit, time_offset ) ! Empty routine 605 INTEGER, INTENT( in ) :: kt ! Dummy argument empty routine 606 INTEGER, INTENT( in ),OPTIONAL :: kit ! Dummy argument empty routine 607 INTEGER, INTENT( in ),OPTIONAL :: time_offset ! Dummy argument empty routine 608 WRITE(*,*) 'bdy_dta_tides: You should not have seen this print! error?', kt, jit 609 END SUBROUTINE bdy_dta_tides 462 610 #endif 463 611 464 612 !!====================================================================== 465 613 END MODULE bdytides 614 -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90
r3777 r4292 20 20 USE dom_oce ! ocean space and time domain variables 21 21 USE bdy_oce ! ocean open boundary conditions 22 USE bdylib ! for orlanski library routines 22 23 USE bdydta, ONLY: bf 23 24 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 51 52 DO ib_bdy=1, nb_bdy 52 53 53 SELECT CASE( nn_tra(ib_bdy) )54 CASE( jp_none)54 SELECT CASE( cn_tra(ib_bdy) ) 55 CASE('none') 55 56 CYCLE 56 CASE( jp_frs)57 CASE('frs') 57 58 CALL bdy_tra_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 58 CASE( 2)59 CASE('specified') 59 60 CALL bdy_tra_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 60 CASE( 3)61 CASE('neumann') 61 62 CALL bdy_tra_nmn( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 62 CASE(4) 63 CASE('orlanski') 64 CALL bdy_tra_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ll_npo=.false. ) 65 CASE('orlanski_npo') 66 CALL bdy_tra_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ll_npo=.true. ) 67 CASE('runoff') 63 68 CALL bdy_tra_rnf( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 64 69 CASE DEFAULT … … 196 201 ! 197 202 END SUBROUTINE bdy_tra_nmn 203 204 205 SUBROUTINE bdy_tra_orlanski( idx, dta, ll_npo ) 206 !!---------------------------------------------------------------------- 207 !! *** SUBROUTINE bdy_tra_orlanski *** 208 !! 209 !! - Apply Orlanski radiation to temperature and salinity. 210 !! - Wrapper routine for bdy_orlanski_3d 211 !! 212 !! 213 !! References: Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001) 214 !!---------------------------------------------------------------------- 215 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 216 TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data 217 LOGICAL, INTENT(in) :: ll_npo ! switch for NPO version 218 219 INTEGER :: igrd ! grid index 220 !!---------------------------------------------------------------------- 221 222 IF( nn_timing == 1 ) CALL timing_start('bdy_tra_orlanski') 223 ! 224 igrd = 1 ! Orlanski bc on temperature; 225 ! 226 CALL bdy_orlanski_3d( idx, igrd, tsb(:,:,:,jp_tem), tsa(:,:,:,jp_tem), dta%tem, ll_npo ) 227 228 igrd = 1 ! Orlanski bc on salinity; 229 ! 230 CALL bdy_orlanski_3d( idx, igrd, tsb(:,:,:,jp_sal), tsa(:,:,:,jp_sal), dta%sal, ll_npo ) 231 ! 232 IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_orlanski') 233 ! 234 235 END SUBROUTINE bdy_tra_orlanski 236 198 237 199 238 SUBROUTINE bdy_tra_rnf( idx, dta, kt ) -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90
r3294 r4292 104 104 ii = idx%nbi(jb,jgrd) 105 105 ij = idx%nbj(jb,jgrd) 106 zubtpecor = zubtpecor + idx%flagu(jb ) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk)106 zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 107 107 END DO 108 108 END DO … … 112 112 ii = idx%nbi(jb,jgrd) 113 113 ij = idx%nbj(jb,jgrd) 114 zubtpecor = zubtpecor + idx%flagv(jb ) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk)114 zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 115 115 END DO 116 116 END DO … … 136 136 ii = idx%nbi(jb,jgrd) 137 137 ij = idx%nbj(jb,jgrd) 138 ua(ii,ij,jk) = ua(ii,ij,jk) - idx%flagu(jb ) * zubtpecor * umask(ii,ij,jk)139 ztranst = ztranst + idx%flagu(jb ) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk)138 ua(ii,ij,jk) = ua(ii,ij,jk) - idx%flagu(jb,jgrd) * zubtpecor * umask(ii,ij,jk) 139 ztranst = ztranst + idx%flagu(jb,jgrd) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk) 140 140 END DO 141 141 END DO … … 145 145 ii = idx%nbi(jb,jgrd) 146 146 ij = idx%nbj(jb,jgrd) 147 va(ii,ij,jk) = va(ii,ij,jk) -idx%flagv(jb ) * zubtpecor * vmask(ii,ij,jk)148 ztranst = ztranst + idx%flagv(jb ) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk)147 va(ii,ij,jk) = va(ii,ij,jk) -idx%flagv(jb,jgrd) * zubtpecor * vmask(ii,ij,jk) 148 ztranst = ztranst + idx%flagv(jb,jgrd) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 149 149 END DO 150 150 END DO -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
r3294 r4292 196 196 thick0(:,:) = 0._wp 197 197 DO jk = 1, jpkm1 198 vol0 = vol0 + SUM( area (:,:) * tmask(:,:,jk) * fse3t_0(:,:,jk) )199 thick0(:,:) = thick0(:,:) + tmask_i(:,:) * tmask(:,:,jk) * fse3t_0(:,:,jk)198 vol0 = vol0 + SUM( area (:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) ) 199 thick0(:,:) = thick0(:,:) + tmask_i(:,:) * tmask(:,:,jk) * e3t_0(:,:,jk) 200 200 END DO 201 201 IF( lk_mpp ) CALL mpp_sum( vol0 ) … … 212 212 ik = mbkt(ji,jj) 213 213 IF( ik > 1 ) THEN 214 zztmp = ( gdept_ 0(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_0(ik) - gdept_0(ik-1) )214 zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) ) 215 215 sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1) 216 216 ENDIF -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diadimg.F90
r3294 r4292 112 112 113 113 CASE ( 'T') 114 z4dep(:)=gdept_ 0(:)114 z4dep(:)=gdept_1d(:) 115 115 116 116 CASE ( 'W' ) 117 z4dep(:)=gdepw_ 0(:)117 z4dep(:)=gdepw_1d(:) 118 118 119 119 CASE ( '2' ) -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90
r4147 r4292 1 1 MODULE diaharm 2 3 #if defined key_diaharm && defined key_tide 4 !!================================================================================= 2 !!====================================================================== 5 3 !! *** MODULE diaharm *** 6 4 !! Harmonic analysis of tidal constituents 7 !!================================================================================= 8 !! * Modules used 5 !!====================================================================== 6 !! History : 3.1 ! 2007 (O. Le Galloudec, J. Chanut) Original code 7 !!---------------------------------------------------------------------- 8 #if defined key_diaharm && defined key_tide 9 !!---------------------------------------------------------------------- 10 !! 'key_diaharm' 11 !! 'key_tide' 12 !!---------------------------------------------------------------------- 9 13 USE oce ! ocean dynamics and tracers variables 10 14 USE dom_oce ! ocean space and time domain 11 USE in_out_manager ! I/O units12 USE lbclnk ! ocean lateral boundary conditions (or mpp link)13 USE ioipsl ! NetCDF IPSL library14 USE diadimg ! To write dimg15 15 USE phycst 16 16 USE dynspg_oce … … 18 18 USE daymod 19 19 USE tide_mod 20 USE iom 20 USE in_out_manager ! I/O units 21 USE iom ! I/0 library 22 USE ioipsl ! NetCDF IPSL library 23 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 24 USE diadimg ! To write dimg 21 25 USE timing ! preformance summary 22 26 USE wrk_nemo ! working arrays … … 30 34 INTEGER, PARAMETER :: jpdimsparse = jpincomax*300*24 31 35 32 INTEGER :: & !! namelist variables 33 nit000_han , & ! First time step used for harmonic analysis 34 nitend_han , & ! Last time step used for harmonic analysis 35 nstep_han , & ! Time step frequency for harmonic analysis 36 nb_ana ! Number of harmonics to analyse 37 38 INTEGER , ALLOCATABLE, DIMENSION(:) :: name 39 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp 40 REAL(wp), ALLOCATABLE, DIMENSION(:) :: ana_freq, vt, ut, ft 41 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: out_eta, & 42 out_u , & 43 out_v 44 45 INTEGER :: ninco, nsparse 46 INTEGER , DIMENSION(jpdimsparse) :: njsparse, nisparse 47 INTEGER , SAVE, DIMENSION(jpincomax) :: ipos1 48 REAL(wp), DIMENSION(jpdimsparse) :: valuesparse 49 REAL(wp), DIMENSION(jpincomax) :: ztmp4 , ztmp7 50 REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) :: ztmp3 , zpilier 51 REAL(wp), SAVE, DIMENSION(jpincomax) :: zpivot 52 53 CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: & 54 tname ! Names of tidal constituents ('M2', 'K1',...) 55 56 57 !! * Routine accessibility 58 PUBLIC dia_harm ! routine called by step.F90 59 60 !!--------------------------------------------------------------------------------- 61 !! 62 !!--------------------------------------------------------------------------------- 63 36 ! !!!namelist variables 37 INTEGER :: nit000_han ! First time step used for harmonic analysis 38 INTEGER :: nitend_han ! Last time step used for harmonic analysis 39 INTEGER :: nstep_han ! Time step frequency for harmonic analysis 40 INTEGER :: nb_ana ! Number of harmonics to analyse 41 42 INTEGER , ALLOCATABLE, DIMENSION(:) :: name 43 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp 44 REAL(wp), ALLOCATABLE, DIMENSION(:) :: ana_freq, ut , vt , ft 45 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: out_eta , out_u, out_v 46 47 INTEGER :: ninco, nsparse 48 INTEGER , DIMENSION(jpdimsparse) :: njsparse, nisparse 49 INTEGER , SAVE, DIMENSION(jpincomax) :: ipos1 50 REAL(wp), DIMENSION(jpdimsparse) :: valuesparse 51 REAL(wp), DIMENSION(jpincomax) :: ztmp4 , ztmp7 52 REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) :: ztmp3 , zpilier 53 REAL(wp), SAVE, DIMENSION(jpincomax) :: zpivot 54 55 CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: tname ! Names of tidal constituents ('M2', 'K1',...) 56 57 PUBLIC dia_harm ! routine called by step.F90 58 59 !!---------------------------------------------------------------------- 60 !! NEMO/OPA 3.5 , NEMO Consortium (2013) 61 !! $Id:$ 62 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 63 !!---------------------------------------------------------------------- 64 64 CONTAINS 65 65 … … 67 67 !!---------------------------------------------------------------------- 68 68 !! *** ROUTINE dia_harm_init *** 69 !!----------------------------------------------------------------------70 69 !! 71 70 !! ** Purpose : Initialization of tidal harmonic analysis … … 73 72 !! ** Method : Initialize frequency array and nodal factor for nit000_han 74 73 !! 75 !! History : 76 !! 9.0 O. Le Galloudec and J. Chanut (Original) 77 !!-------------------------------------------------------------------- 78 !! * Local declarations 74 !!-------------------------------------------------------------------- 79 75 INTEGER :: jh, nhan, jk, ji 80 76 INTEGER :: ios ! Local integer output status for namelist read … … 108 104 ! Basic checks on harmonic analysis time window: 109 105 ! ---------------------------------------------- 110 IF (nit000 > nit000_han) THEN 111 IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : nit000_han must be greater than nit000, stop' 112 IF(lwp) WRITE(numout,*) ' restart capability not implemented' 113 nstop = nstop + 1 114 ENDIF 115 IF (nitend < nitend_han) THEN 116 IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : nitend_han must be lower than nitend, stop' 117 IF(lwp) WRITE(numout,*) ' restart capability not implemented' 118 nstop = nstop + 1 119 ENDIF 120 121 IF (MOD(nitend_han-nit000_han+1,nstep_han).NE.0) THEN 122 IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : analysis time span must be a multiple of nstep_han, stop' 123 nstop = nstop + 1 124 END IF 125 126 nb_ana=0 106 IF( nit000 > nit000_han ) CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000', & 107 & ' restart capability not implemented' ) 108 IF( nitend < nitend_han ) CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend', & 109 & 'restart capability not implemented' ) 110 111 IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 ) & 112 & CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' ) 113 114 nb_ana = 0 127 115 DO jk=1,jpmax_harmo 128 116 DO ji=1,jpmax_harmo … … 157 145 ! Initialize frequency array: 158 146 ! --------------------------- 159 ALLOCATE(ana_freq(nb_ana)) 160 ALLOCATE(vt (nb_ana)) 161 ALLOCATE(ut (nb_ana)) 162 ALLOCATE(ft (nb_ana)) 163 164 CALL tide_harmo(ana_freq, vt, ut , ft, name ,nb_ana) 147 ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) ) 148 149 CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana ) 165 150 166 151 IF(lwp) WRITE(numout,*) 'Analysed frequency : ',nb_ana ,'Frequency ' … … 172 157 ! Initialize temporary arrays: 173 158 ! ---------------------------- 174 ALLOCATE( ana_temp(jpi,jpj, nb_ana*2,3))159 ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) ) 175 160 ana_temp(:,:,:,:) = 0.e0 176 161 177 162 END SUBROUTINE dia_harm_init 178 163 164 179 165 SUBROUTINE dia_harm ( kt ) 180 166 !!---------------------------------------------------------------------- 181 167 !! *** ROUTINE dia_harm *** 182 !!----------------------------------------------------------------------183 168 !! 184 169 !! ** Purpose : Tidal harmonic analysis main routine … … 186 171 !! ** Action : Sums ssh/u/v over time analysis [nit000_han,nitend_han] 187 172 !! 188 !! History : 189 !! 9.0 O. Le Galloudec and J. Chanut (Original) 190 !!-------------------------------------------------------------------- 191 !! * Argument: 173 !!-------------------------------------------------------------------- 192 174 INTEGER, INTENT( IN ) :: kt 193 194 !! * Local declarations 175 ! 195 176 INTEGER :: ji, jj, jh, jc, nhc 196 177 REAL(wp) :: ztime, ztemp … … 198 179 IF( nn_timing == 1 ) CALL timing_start('dia_harm') 199 180 200 IF ( kt .EQ.nit000 ) CALL dia_harm_init181 IF ( kt == nit000 ) CALL dia_harm_init 201 182 202 183 IF ( ((kt.GE.nit000_han).AND.(kt.LE.nitend_han)).AND. & … … 215 196 DO ji = 1,jpi 216 197 ! Elevation 217 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) & 218 + ztemp*sshn(ji,jj)*tmask(ji,jj,1) 198 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj) *tmask(ji,jj,1) 219 199 #if defined key_dynspg_ts 220 ! ubar 221 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) & 222 + ztemp*un_b(ji,jj)*hur(ji,jj)*umask(ji,jj,1) 223 ! vbar 224 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) & 225 + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask(ji,jj,1) 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) 226 202 #endif 227 203 END DO … … 233 209 END IF 234 210 235 IF ( kt .EQ. nitend_han )CALL dia_harm_end211 IF ( kt == nitend_han ) CALL dia_harm_end 236 212 237 213 IF( nn_timing == 1 ) CALL timing_stop('dia_harm') … … 239 215 END SUBROUTINE dia_harm 240 216 217 241 218 SUBROUTINE dia_harm_end 242 219 !!---------------------------------------------------------------------- 243 220 !! *** ROUTINE diaharm_end *** 244 !!----------------------------------------------------------------------245 221 !! 246 222 !! ** Purpose : Compute the Real and Imaginary part of tidal constituents … … 248 224 !! ** Action : Decompose the signal on the harmonic constituents 249 225 !! 250 !! History : 251 !! 9.0 O. Le Galloudec and J. Chanut (Original) 252 !!-------------------------------------------------------------------- 253 254 !! * Local declarations 226 !!-------------------------------------------------------------------- 255 227 INTEGER :: ji, jj, jh, jc, jn, nhan, jl 256 228 INTEGER :: ksp, kun, keq … … 283 255 nisparse(ksp) = keq 284 256 njsparse(ksp) = kun 285 valuesparse(ksp)= & 286 +( MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) & 287 +(1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh))) 257 valuesparse(ksp) = ( MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) & 258 & + (1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)) ) 288 259 END DO 289 260 END DO 290 261 END DO 291 262 292 nsparse =ksp263 nsparse = ksp 293 264 294 265 ! Elevation: … … 296 267 DO ji = 1, jpi 297 268 ! Fill input array 298 kun =0299 DO jh = 1, nb_ana300 DO jc = 1, 2269 kun = 0 270 DO jh = 1, nb_ana 271 DO jc = 1, 2 301 272 kun = kun + 1 302 273 ztmp4(kun)=ana_temp(ji,jj,kun,1) 303 END DO304 END DO274 END DO 275 END DO 305 276 306 277 CALL SUR_DETERMINE(jj) … … 314 285 END DO 315 286 316 ALLOCATE( out_eta(jpi,jpj,2*nb_ana))317 ALLOCATE(out_u (jpi,jpj,2*nb_ana))318 ALLOCATE(out_v (jpi,jpj,2*nb_ana))287 ALLOCATE( out_eta(jpi,jpj,2*nb_ana), & 288 & out_u (jpi,jpj,2*nb_ana), & 289 & out_v (jpi,jpj,2*nb_ana) ) 319 290 320 291 DO jj = 1, jpj 321 292 DO ji = 1, jpi 322 293 DO jh = 1, nb_ana 323 X1 =ana_amp(ji,jj,jh,1)324 X2 =-ana_amp(ji,jj,jh,2)325 out_eta(ji,jj,jh )=X1 * tmask(ji,jj,1)326 out_eta(ji,jj, nb_ana+jh)=X2 * tmask(ji,jj,1)294 X1 = ana_amp(ji,jj,jh,1) 295 X2 =-ana_amp(ji,jj,jh,2) 296 out_eta(ji,jj,jh ) = X1 * tmask(ji,jj,1) 297 out_eta(ji,jj,jh+nb_ana) = X2 * tmask(ji,jj,1) 327 298 ENDDO 328 299 ENDDO … … 402 373 END SUBROUTINE dia_harm_end 403 374 375 404 376 SUBROUTINE dia_wri_harm 405 377 !!-------------------------------------------------------------------- 406 378 !! *** ROUTINE dia_wri_harm *** 407 !!--------------------------------------------------------------------408 379 !! 409 380 !! ** Purpose : Write tidal harmonic analysis results in a netcdf file 410 !! 411 !! 412 !! History : 413 !! 9.0 O. Le Galloudec and J. Chanut (Original) 414 !!-------------------------------------------------------------------- 415 416 !! * Local declarations 381 !!-------------------------------------------------------------------- 417 382 CHARACTER(LEN=lc) :: cltext 418 383 CHARACTER(LEN=lc) :: & … … 472 437 #else 473 438 DO jh = 1, nb_ana 474 CALL iom_put( TRIM(tname(jh))//'x_v', out_u(:,:,jh) )475 CALL iom_put( TRIM(tname(jh))//'y_v', out_u(:,:,nb_ana+jh) )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) ) 476 441 END DO 477 442 #endif 478 443 479 444 END SUBROUTINE dia_wri_harm 445 480 446 481 447 SUBROUTINE SUR_DETERMINE(init) … … 486 452 !! 487 453 !!--------------------------------------------------------------------------------- 488 INTEGER, INTENT(in) :: init489 454 INTEGER, INTENT(in) :: init 455 ! 490 456 INTEGER :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd 491 457 REAL(wp) :: zval1, zval2, zx1 … … 496 462 CALL wrk_alloc( jpincomax , ipos2 , ipivot ) 497 463 498 IF( init==1 )THEN 499 500 IF( nsparse .GT. jpdimsparse ) & 501 CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse') 502 503 IF( ninco .GT. jpincomax ) & 504 CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax') 505 506 ztmp3(:,:)=0.e0 507 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 ! 508 470 DO jk1_sd = 1, nsparse 509 471 DO jk2_sd = 1, nsparse 510 511 nisparse(jk2_sd)=nisparse(jk2_sd) 512 njsparse(jk2_sd)=njsparse(jk2_sd) 513 472 nisparse(jk2_sd) = nisparse(jk2_sd) 473 njsparse(jk2_sd) = njsparse(jk2_sd) 514 474 IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN 515 475 ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) & 516 476 + valuesparse(jk1_sd)*valuesparse(jk2_sd) 517 477 ENDIF 518 519 ENDDO 520 ENDDO 478 END DO 479 END DO 521 480 522 481 DO jj_sd = 1 ,ninco … … 588 547 ENDDO 589 548 590 591 549 CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 ) 592 550 CALL wrk_dealloc( jpincomax , ipos2 , ipivot ) … … 594 552 END SUBROUTINE SUR_DETERMINE 595 553 596 597 554 #else 598 555 !!---------------------------------------------------------------------- … … 601 558 LOGICAL, PUBLIC, PARAMETER :: lk_diaharm = .FALSE. 602 559 CONTAINS 603 604 560 SUBROUTINE dia_harm ( kt ) ! Empty routine 605 561 INTEGER, INTENT( IN ) :: kt 606 562 WRITE(*,*) 'dia_harm: you should not have seen this print' 607 563 END SUBROUTINE dia_harm 608 609 610 #endif 564 #endif 565 611 566 !!====================================================================== 612 567 END MODULE diaharm -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90
r3764 r4292 304 304 ! ----------------------------- ! 305 305 306 ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_ 0to do this search...)306 ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_1d to do this search...) 307 307 ilevel = 0 308 308 zthick_0 = 0._wp 309 309 DO jk = 1, jpkm1 310 zthick_0 = zthick_0 + e3t_ 0(jk)310 zthick_0 = zthick_0 + e3t_1d(jk) 311 311 IF( zthick_0 < 300. ) ilevel = jk 312 312 END DO … … 326 326 htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) * MIN( fse3t(ji,jj,ilevel+1), zthick(ji,jj) ) & 327 327 * tmask(ji,jj,ilevel+1) 328 htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) * MIN( fse3t(ji,jj,ilevel+1), zthick(ji,jj) ) &329 & * tmask(ji,jj,ilevel+1)330 328 END DO 331 329 END DO -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90
r4148 r4292 259 259 ! 260 260 #if defined key_mpp_mpi 261 ijpjjpk = jpj*jpk 261 262 ish(1) = ijpjjpk ; ish2(1) = jpj ; ish2(2) = jpk 262 263 zwork(1:ijpjjpk) = RESHAPE( p_fval, ish ) … … 314 315 END DO 315 316 #if defined key_mpp_mpi 317 ijpjjpk = jpj*jpk 316 318 ish(1) = jpj*jpk ; ish2(1) = jpj ; ish2(2) = jpk 317 319 zwork(1:ijpjjpk)= RESHAPE( p_fval, ish ) … … 670 672 CALL histbeg(clhstnam, 1, zfoo, jpj, zphi, & 671 673 1, 1, 1, jpj, niter, zjulian, zdt*nn_fptr, nhoridz, numptr, domain_id=nidom_ptr) 672 ! Vertical grids : gdept_ 0, gdepw_0674 ! Vertical grids : gdept_1d, gdepw_1d 673 675 CALL histvert( numptr, "deptht", "Vertical T levels", & 674 & "m", jpk, gdept_ 0, ndepidzt, "down" )676 & "m", jpk, gdept_1d, ndepidzt, "down" ) 675 677 CALL histvert( numptr, "depthw", "Vertical W levels", & 676 & "m", jpk, gdepw_ 0, ndepidzw, "down" )678 & "m", jpk, gdepw_1d, ndepidzw, "down" ) 677 679 ! 678 680 CALL wheneq ( jpj*jpk, MIN(sjk(:,:,1), 1._wp), 1, 1., ndex , ndim ) ! Lat-Depth -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
r4161 r4292 25 25 USE oce ! ocean dynamics and tracers 26 26 USE dom_oce ! ocean space and time domain 27 USE dynadv, ONLY: ln_dynadv_vec 27 28 USE zdf_oce ! ocean vertical physics 28 29 USE ldftra_oce ! ocean active tracers: lateral physics … … 44 45 USE diadimg ! dimg direct access file format output 45 46 USE diaar5, ONLY : lk_diaar5 47 USE dynadv, ONLY : ln_dynadv_vec 46 48 USE iom 47 49 USE ioipsl … … 144 146 ENDIF 145 147 146 CALL iom_put( "toce" , tsn(:,:,:,jp_tem) ) ! temperature 147 CALL iom_put( "soce" , tsn(:,:,:,jp_sal) ) ! salinity 148 CALL iom_put( "sst" , tsn(:,:,1,jp_tem) ) ! sea surface temperature 149 CALL iom_put( "sst2" , tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) ) ! square of sea surface temperature 150 CALL iom_put( "sss" , tsn(:,:,1,jp_sal) ) ! sea surface salinity 151 CALL iom_put( "sss2" , tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) ) ! square of sea surface salinity 152 CALL iom_put( "uoce" , un ) ! i-current 153 CALL iom_put( "suoce" , un(:,:,1) ) ! surface i-current 154 CALL iom_put( "voce" , vn ) ! j-current 155 CALL iom_put( "svoce" , vn(:,:,1) ) ! surface j-current 156 157 CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef. 158 CALL iom_put( "avm" , avmu ) ! T vert. eddy visc. coef. 148 IF( lk_vvl ) THEN 149 z3d(:,:,:) = tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) 150 CALL iom_put( "toce" , z3d ) ! heat content 151 CALL iom_put( "sst" , z3d(:,:,1) ) ! sea surface heat content 152 z3d(:,:,1) = tsn(:,:,1,jp_tem) * z3d(:,:,1) 153 CALL iom_put( "sst2" , z3d(:,:,1) ) ! sea surface content of squared temperature 154 z3d(:,:,:) = tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) 155 CALL iom_put( "soce" , z3d ) ! salinity content 156 CALL iom_put( "sss" , z3d(:,:,1) ) ! sea surface salinity content 157 z3d(:,:,1) = tsn(:,:,1,jp_sal) * z3d(:,:,1) 158 CALL iom_put( "sss2" , z3d(:,:,1) ) ! sea surface content of squared salinity 159 ELSE 160 CALL iom_put( "toce" , tsn(:,:,:,jp_tem) ) ! temperature 161 CALL iom_put( "sst" , tsn(:,:,1,jp_tem) ) ! sea surface temperature 162 CALL iom_put( "sst2" , tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) ) ! square of sea surface temperature 163 CALL iom_put( "soce" , tsn(:,:,:,jp_sal) ) ! salinity 164 CALL iom_put( "sss" , tsn(:,:,1,jp_sal) ) ! sea surface salinity 165 CALL iom_put( "sss2" , tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) ) ! square of sea surface salinity 166 END IF 167 IF( lk_vvl .AND. (.NOT. ln_dynadv_vec) ) THEN 168 CALL iom_put( "uoce" , un(:,:,:) * fse3u_n(:,:,:) ) ! i-transport 169 CALL iom_put( "voce" , vn(:,:,:) * fse3v_n(:,:,:) ) ! j-transport 170 ELSE 171 CALL iom_put( "uoce" , un ) ! i-current 172 CALL iom_put( "voce" , vn ) ! j-current 173 END IF 174 CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef. 175 CALL iom_put( "avm" , avmu ) ! T vert. eddy visc. coef. 159 176 IF( lk_zdfddm ) THEN 160 177 CALL iom_put( "avs" , fsavs(:,:,:) ) ! S vert. eddy diff. coef. … … 252 269 ! 253 270 CALL wrk_alloc( jpi , jpj , zw2d ) 254 IF ( ln_traldf_gdia ) call wrk_alloc( jpi , jpj , jpk , zw3d )271 IF ( ln_traldf_gdia .OR. lk_vvl ) call wrk_alloc( jpi , jpj , jpk , zw3d ) 255 272 ! 256 273 ! Output the initial state and forcings … … 325 342 & nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set ) 326 343 CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept 327 & "m", ipk, gdept_ 0, nz_T, "down" )344 & "m", ipk, gdept_1d, nz_T, "down" ) 328 345 ! ! Index of ocean points 329 346 CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume … … 361 378 & nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set ) 362 379 CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept 363 & "m", ipk, gdept_ 0, nz_U, "down" )380 & "m", ipk, gdept_1d, nz_U, "down" ) 364 381 ! ! Index of ocean points 365 382 CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U ) ! volume … … 374 391 & nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set ) 375 392 CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept 376 & "m", ipk, gdept_ 0, nz_V, "down" )393 & "m", ipk, gdept_1d, nz_V, "down" ) 377 394 ! ! Index of ocean points 378 395 CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume … … 387 404 & nit000-1, zjulian, zdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set ) 388 405 CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw 389 & "m", ipk, gdepw_ 0, nz_W, "down" )406 & "m", ipk, gdepw_1d, nz_W, "down" ) 390 407 391 408 … … 397 414 CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn 398 415 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 416 IF( lk_vvl ) THEN 417 CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t_n 418 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 419 CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t_n 420 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 421 CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t_n 422 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 423 ENDIF 399 424 ! !!! nid_T : 2D 400 425 CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst … … 408 433 CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx 409 434 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 410 #if ! defined key_vvl 411 CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"& ! emp * tsn(:,:,1,jp_tem)435 IF( .NOT. lk_vvl ) THEN 436 CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * tsn(:,:,1,jp_tem) 412 437 & , "KgC/m2/s", & ! sosst_cd 413 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )414 CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"& ! emp * tsn(:,:,1,jp_sal)438 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 439 CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * tsn(:,:,1,jp_sal) 415 440 & , "KgPSU/m2/s",& ! sosss_cd 416 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )417 #endif 441 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 442 ENDIF 418 443 CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr 419 444 & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) … … 587 612 ! --------------------- 588 613 589 ! ndex(1) est utilise ssi l'avant dernier argument est diff ferent de614 ! ndex(1) est utilise ssi l'avant dernier argument est different de 590 615 ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument 591 616 ! donne le nombre d'elements, et ndex la liste des indices a sortir … … 597 622 598 623 ! Write fields on T grid 599 CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem), ndim_T , ndex_T ) ! temperature 600 CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal), ndim_T , ndex_T ) ! salinity 601 CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem), ndim_hT, ndex_hT ) ! sea surface temperature 602 CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal), ndim_hT, ndex_hT ) ! sea surface salinity 624 IF( lk_vvl ) THEN 625 CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content 626 CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content 627 CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content 628 CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content 629 ELSE 630 CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T ) ! temperature 631 CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T ) ! salinity 632 CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT ) ! sea surface temperature 633 CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT ) ! sea surface salinity 634 635 ENDIF 636 IF( lk_vvl ) THEN 637 zw3d(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 638 CALL histwrite( nid_T, "vovvle3t", it, fse3t_n (:,:,:) , ndim_T , ndex_T ) ! level thickness 639 CALL histwrite( nid_T, "vovvldep", it, fsdept_n(:,:,:) , ndim_T , ndex_T ) ! t-point depth 640 CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation 641 ENDIF 603 642 CALL histwrite( nid_T, "sossheig", it, sshn , ndim_hT, ndex_hT ) ! sea surface height 604 643 CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux … … 606 645 ! (includes virtual salt flux beneath ice 607 646 ! in linear free surface case) 608 #if ! defined key_vvl 609 zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)610 CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )! c/d term on sst611 zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)612 CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )! c/d term on sss613 #endif 647 IF( .NOT. lk_vvl ) THEN 648 zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem) 649 CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst 650 zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal) 651 CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss 652 ENDIF 614 653 CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux 615 654 CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux … … 752 791 ! 753 792 CALL wrk_dealloc( jpi , jpj , zw2d ) 754 IF ( ln_traldf_gdia ) call wrk_dealloc( jpi , jpj , jpk , zw3d )793 IF ( ln_traldf_gdia .OR. lk_vvl ) call wrk_dealloc( jpi , jpj , jpk , zw3d ) 755 794 ! 756 795 IF( nn_timing == 1 ) CALL timing_stop('dia_wri') … … 813 852 1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit 814 853 CALL histvert( id_i, "deptht", "Vertical T levels", & ! Vertical grid : gdept 815 "m", jpk, gdept_ 0, nz_i, "down")854 "m", jpk, gdept_1d, nz_i, "down") 816 855 817 856 ! Declare all the output fields as NetCDF variables … … 841 880 CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2" , & ! j-wind stress 842 881 & jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout ) 882 IF( lk_vvl ) THEN 883 CALL histdef( id_i, "vovvldep", "T point depth" , "m" , & ! t-point depth 884 & jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 885 END IF 843 886 844 887 #if defined key_lim2 -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r4247 r4292 152 152 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gphit, gphiu !: latitude of t-, u-, v- and f-points (degre) 153 153 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gphiv, gphif !: 154 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1t, e2t !: horizontal scale factors at t-point (m)155 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1u, e2u !: horizontal scale factors at u-point (m)156 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1v, e2v !: horizontal scale factors at v-point (m)157 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1f, e2f !: horizontal scale factors at f-point (m)154 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1t, e2t !: horizontal scale factors at t-point (m) 155 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1u, e2u !: horizontal scale factors at u-point (m) 156 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1v, e2v !: horizontal scale factors at v-point (m) 157 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) :: e1f, e2f !: horizontal scale factors at f-point (m) 158 158 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e1e2t !: surface at t-point (m2) 159 159 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ff !: coriolis factor (2.*omega*sin(yphi) ) (s-1) … … 169 169 !! All coordinates 170 170 !! --------------- 171 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdep3w !: depth of T-points (sum of e3w) (m)172 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept , gdepw !: analytical depth at T-Wpoints (m)173 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3v , e3f !: analytical vertical scale factors at V--F174 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t , e3u !: T--Upoints (m)175 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3vw !: analytical vertical scale factors at VW--176 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w , e3uw !: W--UWpoints (m)171 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdep3w_0 !: depth of t-points (sum of e3w) (m) 172 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_0, gdepw_0 !: analytical (time invariant) depth at t-w points (m) 173 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3v_0 , e3f_0 !: analytical (time invariant) vertical scale factors at v-f 174 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_0 , e3u_0 !: t-u points (m) 175 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3vw_0 !: analytical (time invariant) vertical scale factors at vw 176 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w_0 , e3uw_0 !: w-uw points (m) 177 177 #if defined key_vvl 178 178 LOGICAL, PUBLIC, PARAMETER :: lk_vvl = .TRUE. !: variable grid flag … … 180 180 !! All coordinates 181 181 !! --------------- 182 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdep3w_1 !: depth of T-points (sum of e3w) (m) 183 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_1, gdepw_1 !: analytical depth at T-W points (m) 184 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3v_1 , e3f_1 !: analytical vertical scale factors at V--F 185 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_1 , e3u_1 !: T--U points (m) 186 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3vw_1 !: analytical vertical scale factors at VW-- 187 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w_1 , e3uw_1 !: W--UW points (m) 188 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_b !: before - - - - T points (m) 189 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_b , e3v_b !: - - - - - U--V points (m) 182 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdep3w_n !: now depth of T-points (sum of e3w) (m) 183 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_n, gdepw_n !: now depth at T-W points (m) 184 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_n !: now vertical scale factors at t point (m) 185 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_n , e3v_n !: - - - - u --v points (m) 186 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w_n , e3f_n !: - - - - w --f points (m) 187 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3uw_n , e3vw_n !: - - - - uw--vw points (m) 188 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_b !: before - - - - t points (m) 189 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_b , e3v_b !: - - - - - u --v points (m) 190 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3uw_b , e3vw_b !: - - - - - uw--vw points (m) 191 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_a !: after - - - - t point (m) 192 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3u_a , e3v_a !: - - - - - u --v points (m) 190 193 #else 191 194 LOGICAL, PUBLIC, PARAMETER :: lk_vvl = .FALSE. !: fixed grid flag 192 195 #endif 193 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET :: hur , hvr !: inverse of u and v-points ocean depth (1/m) 194 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu , hv !: depth at u- and v-points (meters) 195 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0 , hv_0 !: refernce depth at u- and v-points (meters) 196 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET :: hur , hvr !: inverse of u and v-points ocean depth (1/m) 197 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu , hv !: depth at u- and v-points (meters) 198 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0 !: reference depth at t- points (meters) 199 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0 , hv_0 !: reference depth at u- and v-points (meters) 200 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: re2u_e1u !: scale factor coeffs at u points (e2u/e1u) 201 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: re1v_e2v !: scale factor coeffs at v points (e1v/e2v) 202 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e12t , r1_e12t !: horizontal cell surface and inverse at t points 203 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e12u , r1_e12u !: horizontal cell surface and inverse at u points 204 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e12v , r1_e12v !: horizontal cell surface and inverse at v points 205 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e12f , r1_e12f !: horizontal cell surface and inverse at f points 196 206 197 207 INTEGER, PUBLIC :: nla10 !: deepest W level Above ~10m (nlb10 - 1) … … 200 210 !! z-coordinate with full steps (also used in the other cases as reference z-coordinate) 201 211 !! =-----------------====------ 202 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gdept_0, gdepw_0!: reference depth of t- and w-points (m)203 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: e3t_0 , e3w_0!: reference vertical scale factors at T- and W-pts (m)204 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3tp , e3wp!: ocean bottom level thickness at T and W points212 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gdept_1d, gdepw_1d !: reference depth of t- and w-points (m) 213 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: e3t_1d , e3w_1d !: reference vertical scale factors at T- and W-pts (m) 214 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3tp , e3wp !: ocean bottom level thickness at T and W points 205 215 206 216 !! s-coordinate and hybrid z-s-coordinate 207 217 !! =----------------======--------------- 208 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gsigt, gsigw!: model level depth coefficient at t-, w-levels (analytic)209 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gsi3w!: model level depth coefficient at w-level (sum of gsigw)210 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: esigt, esigw!: vertical scale factor coef. at t-, w-levels211 212 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbatv , hbatf !: ocean depth at the vertical of V--F213 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbatt , hbatu !: T--Upoints (m)214 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: scosrf, scobot !: ocean surface and bottom topographies215 ! ! (if deviating from coordinate surfaces in HYBRID)216 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hifv , hiff !: interface depth between stretching at V--F217 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hift , hifu !: and quasi-uniform spacing T--Upoints (m)218 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rx1 !: Maximum grid stiffness ratio218 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gsigt, gsigw !: model level depth coefficient at t-, w-levels (analytic) 219 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: gsi3w !: model level depth coefficient at w-level (sum of gsigw) 220 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: esigt, esigw !: vertical scale factor coef. at t-, w-levels 221 222 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbatv , hbatf !: ocean depth at the vertical of v--f 223 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbatt , hbatu !: t--u points (m) 224 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: scosrf, scobot !: ocean surface and bottom topographies 225 ! ! (if deviating from coordinate surfaces in HYBRID) 226 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hifv , hiff !: interface depth between stretching at v--f 227 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hift , hifu !: and quasi-uniform spacing t--u points (m) 228 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rx1 !: Maximum grid stiffness ratio 219 229 220 230 !!---------------------------------------------------------------------- 221 231 !! masks, bathymetry 222 232 !! --------------------------------------------------------------------- 223 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbathy !: number of ocean level (=0, 1, ... , jpk-1)224 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt !: vertical index of the bottom last T- ocean level225 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbku, mbkv !: vertical index of the bottom last U- and W- ocean level226 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters)227 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i !: interior domain T-point mask228 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bmask !: land/ocean mask of barotropic stream function229 230 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tmask, umask, vmask, fmask !: land/ocean mask at T-, U-, V- and F-pts233 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbathy !: number of ocean level (=0, 1, ... , jpk-1) 234 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt !: vertical index of the bottom last T- ocean level 235 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbku, mbkv !: vertical index of the bottom last U- and W- ocean level 236 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bathy !: ocean depth (meters) 237 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tmask_i !: interior domain T-point mask 238 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bmask !: land/ocean mask of barotropic stream function 239 240 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask !: land/ocean mask at T-, U-, V- and F-pts 231 241 232 242 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: tpol, fpol !: north fold mask (jperio= 3 or 4) 233 243 234 244 #if defined key_noslip_accurate 235 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,: ) :: npcoa!: ???236 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: nicoa, njcoa!: ???245 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,: ) :: npcoa !: ??? 246 INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: nicoa, njcoa !: ??? 237 247 #endif 238 248 … … 316 326 & glamf(jpi,jpj) , gphif(jpi,jpj) , e1f(jpi,jpj) , e2f(jpi,jpj) , ff (jpi,jpj) , STAT=ierr(3) ) 317 327 ! 318 ALLOCATE( gdep3w (jpi,jpj,jpk) , e3v(jpi,jpj,jpk) , e3f(jpi,jpj,jpk) , &319 & gdept (jpi,jpj,jpk) , e3t(jpi,jpj,jpk) , e3u(jpi,jpj,jpk) , &320 & gdepw (jpi,jpj,jpk) , e3w(jpi,jpj,jpk) , e3vw(jpi,jpj,jpk) , e3uw(jpi,jpj,jpk) , STAT=ierr(4) )328 ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) , & 329 & gdept_0 (jpi,jpj,jpk) , e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) , & 330 & gdepw_0 (jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , STAT=ierr(4) ) 321 331 ! 322 332 #if defined key_vvl 323 ALLOCATE( gdep3w_1(jpi,jpj,jpk) , e3v_1(jpi,jpj,jpk) , e3f_1 (jpi,jpj,jpk) , & 324 & gdept_1 (jpi,jpj,jpk) , e3t_1(jpi,jpj,jpk) , e3u_1 (jpi,jpj,jpk) , & 325 & gdepw_1 (jpi,jpj,jpk) , e3w_1(jpi,jpj,jpk) , e3vw_1(jpi,jpj,jpk) , e3uw_1(jpi,jpj,jpk) , & 326 & e3t_b (jpi,jpj,jpk) , e3u_b(jpi,jpj,jpk) , e3v_b (jpi,jpj,jpk) , STAT=ierr(5) ) 327 #endif 328 ! 329 ALLOCATE( hu(jpi,jpj) , hur(jpi,jpj) , hu_0(jpi,jpj) , & 330 & hv(jpi,jpj) , hvr(jpi,jpj) , hv_0(jpi,jpj) , STAT=ierr(6) ) 331 ! 332 ALLOCATE( gdept_0(jpk) , gdepw_0(jpk) , & 333 & e3t_0 (jpk) , e3w_0 (jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) , & 334 & gsigt (jpk) , gsigw (jpk) , gsi3w(jpk) , & 335 & esigt (jpk) , esigw (jpk) , STAT=ierr(7) ) 333 ALLOCATE( gdep3w_n(jpi,jpj,jpk) , e3t_n (jpi,jpj,jpk) , e3u_n (jpi,jpj,jpk) , & 334 & gdept_n (jpi,jpj,jpk) , e3v_n (jpi,jpj,jpk) , e3w_n (jpi,jpj,jpk) , & 335 & gdepw_n (jpi,jpj,jpk) , e3f_n (jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) , e3uw_n(jpi,jpj,jpk) , & 336 & e3t_b (jpi,jpj,jpk) , e3u_b (jpi,jpj,jpk) , e3v_b (jpi,jpj,jpk) , & 337 & e3uw_b (jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) , & 338 & e3t_a (jpi,jpj,jpk) , e3u_a (jpi,jpj,jpk) , e3v_a (jpi,jpj,jpk) , STAT=ierr(5) ) 339 #endif 340 ! 341 ALLOCATE( hu (jpi,jpj) , hur (jpi,jpj) , hu_0(jpi,jpj) , ht_0 (jpi,jpj) , & 342 & hv (jpi,jpj) , hvr (jpi,jpj) , hv_0(jpi,jpj) , & 343 & re2u_e1u(jpi,jpj) , re1v_e2v(jpi,jpj) , & 344 & e12t (jpi,jpj) , r1_e12t (jpi,jpj) , & 345 & e12u (jpi,jpj) , r1_e12u (jpi,jpj) , & 346 & e12v (jpi,jpj) , r1_e12v (jpi,jpj) , & 347 & e12f (jpi,jpj) , r1_e12f (jpi,jpj) , STAT=ierr(6) ) 348 ! 349 ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , & 350 & e3t_1d (jpk) , e3w_1d (jpk) , e3tp (jpi,jpj), e3wp(jpi,jpj) , & 351 & gsigt (jpk) , gsigw (jpk) , gsi3w(jpk) , & 352 & esigt (jpk) , esigw (jpk) , STAT=ierr(7) ) 336 353 ! 337 354 ALLOCATE( hbatv (jpi,jpj) , hbatf (jpi,jpj) , & -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r4245 r4292 87 87 CALL dom_msk ! Masks 88 88 IF( ln_sco ) CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency 89 IF( lk_vvl ) CALL dom_vvl! Vertical variable mesh89 IF( lk_vvl ) CALL dom_vvl_init ! Vertical variable mesh 90 90 ! 91 91 IF( lk_c1d ) CALL cor_c1d ! 1D configuration: Coriolis set at T-point 92 ! 93 ! - ML - Used in dom_vvl_sf_nxt and lateral diffusion routines 94 ! but could be usefull in many other routines 95 e12t (:,:) = e1t(:,:) * e2t(:,:) 96 e12u (:,:) = e1u(:,:) * e2u(:,:) 97 e12v (:,:) = e1v(:,:) * e2v(:,:) 98 e12f (:,:) = e1f(:,:) * e2f(:,:) 99 r1_e12t (:,:) = 1._wp / e12t(:,:) 100 r1_e12u (:,:) = 1._wp / e12u(:,:) 101 r1_e12v (:,:) = 1._wp / e12v(:,:) 102 r1_e12f (:,:) = 1._wp / e12f(:,:) 103 re2u_e1u(:,:) = e2u(:,:) / e1u(:,:) 104 re1v_e2v(:,:) = e1v(:,:) / e2v(:,:) 92 105 ! 93 106 hu(:,:) = 0._wp ! Ocean depth at U- and V-points 94 107 hv(:,:) = 0._wp 95 108 DO jk = 1, jpk 96 hu(:,:) = hu(:,:) + fse3u (:,:,jk) * umask(:,:,jk)97 hv(:,:) = hv(:,:) + fse3v (:,:,jk) * vmask(:,:,jk)109 hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 110 hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 98 111 END DO 99 112 ! ! Inverse of the local depth … … 407 420 DO jj = 2, jpjm1 408 421 DO jk = 1, jpkm1 409 zr1(1) = umask(ji-1,jj ,jk) *abs( (gdepw (ji ,jj ,jk )-gdepw(ji-1,jj ,jk ) &410 & +gdepw (ji ,jj ,jk+1)-gdepw(ji-1,jj ,jk+1)) &411 & /(gdepw (ji ,jj ,jk )+gdepw(ji-1,jj ,jk ) &412 & -gdepw (ji ,jj ,jk+1)-gdepw(ji-1,jj ,jk+1) + rsmall) )413 zr1(2) = umask(ji ,jj ,jk) *abs( (gdepw (ji+1,jj ,jk )-gdepw(ji ,jj ,jk ) &414 & +gdepw (ji+1,jj ,jk+1)-gdepw(ji ,jj ,jk+1)) &415 & /(gdepw (ji+1,jj ,jk )+gdepw(ji ,jj ,jk ) &416 & -gdepw (ji+1,jj ,jk+1)-gdepw(ji ,jj ,jk+1) + rsmall) )417 zr1(3) = vmask(ji ,jj ,jk) *abs( (gdepw (ji ,jj+1,jk )-gdepw(ji ,jj ,jk ) &418 & +gdepw (ji ,jj+1,jk+1)-gdepw(ji ,jj ,jk+1)) &419 & /(gdepw (ji ,jj+1,jk )+gdepw(ji ,jj ,jk ) &420 & -gdepw (ji ,jj+1,jk+1)-gdepw(ji ,jj ,jk+1) + rsmall) )421 zr1(4) = vmask(ji ,jj-1,jk) *abs( (gdepw (ji ,jj ,jk )-gdepw(ji ,jj-1,jk ) &422 & +gdepw (ji ,jj ,jk+1)-gdepw(ji ,jj-1,jk+1)) &423 & /(gdepw (ji ,jj ,jk )+gdepw(ji ,jj-1,jk ) &424 & -gdepw (ji, jj ,jk+1)-gdepw(ji ,jj-1,jk+1) + rsmall) )422 zr1(1) = umask(ji-1,jj ,jk) *abs( (gdepw_0(ji ,jj ,jk )-gdepw_0(ji-1,jj ,jk ) & 423 & +gdepw_0(ji ,jj ,jk+1)-gdepw_0(ji-1,jj ,jk+1)) & 424 & /(gdepw_0(ji ,jj ,jk )+gdepw_0(ji-1,jj ,jk ) & 425 & -gdepw_0(ji ,jj ,jk+1)-gdepw_0(ji-1,jj ,jk+1) + rsmall) ) 426 zr1(2) = umask(ji ,jj ,jk) *abs( (gdepw_0(ji+1,jj ,jk )-gdepw_0(ji ,jj ,jk ) & 427 & +gdepw_0(ji+1,jj ,jk+1)-gdepw_0(ji ,jj ,jk+1)) & 428 & /(gdepw_0(ji+1,jj ,jk )+gdepw_0(ji ,jj ,jk ) & 429 & -gdepw_0(ji+1,jj ,jk+1)-gdepw_0(ji ,jj ,jk+1) + rsmall) ) 430 zr1(3) = vmask(ji ,jj ,jk) *abs( (gdepw_0(ji ,jj+1,jk )-gdepw_0(ji ,jj ,jk ) & 431 & +gdepw_0(ji ,jj+1,jk+1)-gdepw_0(ji ,jj ,jk+1)) & 432 & /(gdepw_0(ji ,jj+1,jk )+gdepw_0(ji ,jj ,jk ) & 433 & -gdepw_0(ji ,jj+1,jk+1)-gdepw_0(ji ,jj ,jk+1) + rsmall) ) 434 zr1(4) = vmask(ji ,jj-1,jk) *abs( (gdepw_0(ji ,jj ,jk )-gdepw_0(ji ,jj-1,jk ) & 435 & +gdepw_0(ji ,jj ,jk+1)-gdepw_0(ji ,jj-1,jk+1)) & 436 & /(gdepw_0(ji ,jj ,jk )+gdepw_0(ji ,jj-1,jk ) & 437 & -gdepw_0(ji, jj ,jk+1)-gdepw_0(ji ,jj-1,jk+1) + rsmall) ) 425 438 zrxmax = MAXVAL(zr1(1:4)) 426 439 rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax) -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domstp.F90
r2715 r4292 90 90 91 91 DO jk = 1, jpk 92 IF( gdept_ 0(jk) <= rdth ) rdttra(jk) = rdtmin93 IF( gdept_ 0(jk) > rdth ) THEN92 IF( gdept_1d(jk) <= rdth ) rdttra(jk) = rdtmin 93 IF( gdept_1d(jk) > rdth ) THEN 94 94 rdttra(jk) = rdtmin + ( rdtmax - rdtmin ) & 95 * ( EXP( ( gdept_ 0(jk ) - rdth ) / rdth ) - 1. ) &96 / ( EXP( ( gdept_ 0(jpk) - rdth ) / rdth ) - 1. )95 * ( EXP( ( gdept_1d(jk ) - rdth ) / rdth ) - 1. ) & 96 / ( EXP( ( gdept_1d(jpk) - rdth ) / rdth ) - 1. ) 97 97 ENDIF 98 98 IF(lwp) WRITE(numout,"(36x,f5.2,5x,i3)") rdttra(jk)/3600., jk -
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r4153 r4292 6 6 !! History : 2.0 ! 2006-06 (B. Levier, L. Marie) original code 7 7 !! 3.1 ! 2009-02 (G. Madec, M. Leclair, R. Benshila) pure z* coordinate 8 !! ----------------------------------------------------------------------9 #if defined key_vvl 8 !! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: 9 !! vvl option includes z_star and z_tilde coordinates 10 10 !!---------------------------------------------------------------------- 11 11 !! 'key_vvl' variable volume 12 12 !!---------------------------------------------------------------------- 13 !! dom_vvl : defined coefficients to distribute ssh on each layers14 13 !!---------------------------------------------------------------------- 14 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness 15 !! dom_vvl_sf_nxt : Compute next vertical scale factors 16 !! dom_vvl_sf_swp : Swap vertical scale factors and update the vertical grid 17 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 18 !! dom_vvl_rst : read/write restart file 19 !! dom_vvl_ctl : Check the vvl options 20 !! dom_vvl_orca_fix : Recompute some area-weighted interpolations of vertical scale factors 21 !! : to account for manual changes to e[1,2][u,v] in some Straits 22 !!---------------------------------------------------------------------- 23 !! * Modules used 15 24 USE oce ! ocean dynamics and tracers 16 25 USE dom_oce ! ocean space and time domain 17 USE sbc_oce ! surface boundary condition: ocean 18 USE phycst ! physical constants 26 USE sbc_oce ! ocean surface boundary condition 19 27 USE in_out_manager ! I/O manager 28 USE iom ! I/O manager library 29 USE restart ! ocean restart 20 30 USE lib_mpp ! distributed memory computing library 21 31 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 26 36 PRIVATE 27 37 28 PUBLIC dom_vvl ! called by domain.F90 29 PUBLIC dom_vvl_2 ! called by domain.F90 30 PUBLIC dom_vvl_alloc ! called by nemogcm.F90 31 32 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: mut , muu , muv , muf !: 1/H_0 at t-,u-,v-,f-points 33 34 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: r2dt ! vertical profile time-step, = 2 rdttra 35 ! ! except at nit000 (=rdttra) if neuler=0 38 !! * Routine accessibility 39 PUBLIC dom_vvl_init ! called by domain.F90 40 PUBLIC dom_vvl_sf_nxt ! called by step.F90 41 PUBLIC dom_vvl_sf_swp ! called by step.F90 42 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 43 PRIVATE dom_vvl_orca_fix ! called by dom_vvl_interpol 44 45 !!* Namelist nam_vvl 46 LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate 47 LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate 48 LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate 49 LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate 50 LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate 51 LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer 52 ! ! conservation: not used yet 53 REAL(wp) :: rn_ahe3 = 0.0_wp ! thickness diffusion coefficient 54 REAL(wp) :: rn_rst_e3t = 30._wp ! ztilde to zstar restoration timescale [days] 55 REAL(wp) :: rn_lf_cutoff = 5.0_wp ! cutoff frequency for low-pass filter [days] 56 REAL(wp) :: rn_zdef_max = 0.9_wp ! maximum fractional e3t deformation 57 LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints 58 59 !! * Module variables 60 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport 61 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence 62 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors 63 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a ! baroclinic scale factors 64 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors 65 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 36 66 37 67 !! * Substitutions … … 39 69 # include "vectopt_loop_substitute.h90" 40 70 !!---------------------------------------------------------------------- 41 !! NEMO/OPA 4.0 , NEMO Consortium (2011)71 !! NEMO/OPA 3.3 , NEMO-Consortium (2010) 42 72 !! $Id$ 43 73 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 44 74 !!---------------------------------------------------------------------- 45 CONTAINS 75 76 CONTAINS 46 77 47 78 INTEGER FUNCTION dom_vvl_alloc() 48 79 !!---------------------------------------------------------------------- 49 !! *** ROUTINE dom_vvl_alloc *** 50 !!---------------------------------------------------------------------- 80 !! *** FUNCTION dom_vvl_alloc *** 81 !!---------------------------------------------------------------------- 82 IF( ln_vvl_zstar ) dom_vvl_alloc = 0 83 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 84 ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , & 85 & un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , STAT = dom_vvl_alloc ) 86 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 87 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 88 un_td = 0.0_wp 89 vn_td = 0.0_wp 90 ENDIF 91 IF( ln_vvl_ztilde ) THEN 92 ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc ) 93 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 94 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 95 ENDIF 96 97 END FUNCTION dom_vvl_alloc 98 99 100 SUBROUTINE dom_vvl_init 101 !!---------------------------------------------------------------------- 102 !! *** ROUTINE dom_vvl_init *** 103 !! 104 !! ** Purpose : Initialization of all scale factors, depths 105 !! and water column heights 106 !! 107 !! ** Method : - use restart file and/or initialize 108 !! - interpolate scale factors 109 !! 110 !! ** Action : - fse3t_(n/b) and tilde_e3t_(n/b) 111 !! - Regrid: fse3(u/v)_n 112 !! fse3(u/v)_b 113 !! fse3w_n 114 !! fse3(u/v)w_b 115 !! fse3(u/v)w_n 116 !! fsdept_n, fsdepw_n and fsde3w_n 117 !! - h(t/u/v)_0 118 !! - frq_rst_e3t and frq_rst_hdv 119 !! 120 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 121 !!---------------------------------------------------------------------- 122 USE phycst, ONLY : rpi, rsmall, rad 123 !! * Local declarations 124 INTEGER :: ji,jj,jk 125 INTEGER :: ii0, ii1, ij0, ij1 126 !!---------------------------------------------------------------------- 127 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 128 129 IF(lwp) WRITE(numout,*) 130 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 131 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 132 133 ! choose vertical coordinate (z_star, z_tilde or layer) 134 ! ========================== 135 CALL dom_vvl_ctl 136 137 ! Allocate module arrays 138 ! ====================== 139 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 140 141 ! Read or initialize fse3t_(b/n), tilde_e3t_(b/n) and hdiv_lf (and e3t_a(jpk)) 142 ! ============================================================================ 143 CALL dom_vvl_rst( nit000, 'READ' ) 144 fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 145 146 ! Reconstruction of all vertical scale factors at now and before time steps 147 ! ============================================================================= 148 ! Horizontal scale factor interpolations 149 ! -------------------------------------- 150 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 151 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 152 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,