Changeset 4292


Ignore:
Timestamp:
2013-11-20T17:28:04+01:00 (7 years ago)
Author:
cetlod
Message:

dev_MERGE_2013 : 1st step of the merge, see ticket #1185

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  
    507507!----------------------------------------------------------------------- 
    508508   ln_tide_pot   = .true.   !  use tidal potential forcing 
     509   nb_harmo      =    11    !  number of constituents used 
    509510   clname(1)     =   'M2'   !  name of constituent 
    510511   clname(2)     =   'S2' 
     
    700701   ln_dynadv_cen2= .false. !  flux form - 2nd order centered scheme 
    701702   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) 
    702716/ 
    703717!----------------------------------------------------------------------- 
  • branches/2013/dev_MERGE_2013/NEMOGCM/CONFIG/SHARED/namelist_ref

    r4245 r4292  
    513513!----------------------------------------------------------------------- 
    514514   ln_tide_pot   = .true.   !  use tidal potential forcing 
     515   nb_harmo      =    11    !  number of constituents used 
    515516   clname(1)     =   'M2'   !  name of constituent 
    516517   clname(2)     =   'S2' 
     
    712713   ln_dynadv_cen2= .false. !  flux form - 2nd order centered scheme 
    713714   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) 
    714728/ 
    715729!----------------------------------------------------------------------- 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90

    r4148 r4292  
    2424   USE phycst           ! physical constants 
    2525   USE dom_oce          ! ocean domain 
     26   USE domvvl           ! ocean vertical scale factors 
    2627   USE dom_ice_2        ! LIM-2: ice domain 
    2728   USE ice_2            ! LIM-2: ice variables 
     
    5960   !! * Substitutions 
    6061#  include "vectopt_loop_substitute.h90" 
     62#  include "domzgr_substitute.h90" 
    6163   !!---------------------------------------------------------------------- 
    6264   !! NEMO/LIM2 4.0 , UCL - NEMO Consortium (2011) 
     
    446448      !!------------------------------------------------------------------- 
    447449      ! 
     450      INTEGER :: jk           ! local integer 
     451      ! 
    448452      IF(lwp) WRITE(numout,*) 
    449453      IF(lwp) WRITE(numout,*) 'lim_sbc_init_2 : LIM-2 sea-ice - surface boundary condition' 
     
    472476         snwice_mass  (:,:) = 0.e0           ! no mass exchanges 
    473477         snwice_mass_b(:,:) = 0.e0           ! no mass exchanges 
     478         snwice_fmass (:,:) = 0.e0           ! no mass exchanges 
    474479      ENDIF 
    475480      IF( nn_ice_embd == 2 .AND.          &  ! full embedment (case 2) & no restart :  
     
    477482         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 
    478483         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 
    479515      ENDIF 
    480516      ! 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90

    r4147 r4292  
    237237             
    238238            !  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 ) 
    240240             
    241241            !  calculate oceanic heat flux. 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90

    r4147 r4292  
    4040   INTEGER :: e1u_id, e2v_id, sshn_id, gcb_id 
    4141   INTEGER :: trn_id, trb_id, tra_id 
     42   INTEGER :: unb_id, vnb_id 
    4243 
    4344   !!---------------------------------------------------------------------- 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90

    r3294 r4292  
    2727   USE agrif_opa_sponge 
    2828   USE lib_mpp 
    29    USE wrk_nemo   
     29   USE wrk_nemo 
     30   USE dynspg_oce   
    3031 
    3132   IMPLICIT NONE 
    3233   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 
    3341     
    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 
    3544 
    3645#  include "domzgr_substitute.h90"   
     
    169178      REAL(wp) :: timeref 
    170179      REAL(wp) :: z2dt, znugdt 
    171       REAL(wp) :: zrhox, rhoy 
     180      REAL(wp) :: zrhox, zrhoy 
    172181      REAL(wp), POINTER, DIMENSION(:,:,:) :: zua, zva 
    173182      REAL(wp), POINTER, DIMENSION(:,:) :: spgv1, spgu1, zua2d, zva2d 
     
    180189 
    181190      zrhox = Agrif_Rhox() 
    182       rhoy = Agrif_Rhoy() 
     191      zrhoy = Agrif_Rhoy() 
    183192 
    184193      timeref = 1. 
     
    201210      zva2d = 0. 
    202211 
     212#if defined key_dynspg_flt 
    203213      Agrif_SpecialValue=0. 
    204214      Agrif_UseSpecialValue = ln_spc_dyn 
    205215      CALL Agrif_Bc_variable(zua2d,e1u_id,calledweight=1.,procname=interpu2d) 
    206216      CALL Agrif_Bc_variable(zva2d,e2v_id,calledweight=1.,procname=interpv2d) 
     217#endif 
    207218      Agrif_UseSpecialValue = .FALSE. 
    208219 
     
    210221      IF((nbondi == -1).OR.(nbondi == 2)) THEN 
    211222 
     223#if defined key_dynspg_flt 
    212224         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))) 
    219232               ua(1:2,jj,jk) = ua(1:2,jj,jk) / fse3u(1:2,jj,jk) 
    220233            END DO 
    221234         END DO 
    222235 
     236#if defined key_dynspg_flt 
    223237         DO jk=1,jpkm1 
    224238            DO jj=1,jpj 
     
    240254            ENDIF 
    241255         END DO 
     256#else 
     257         spgu(2,:) = ua_b(2,:) 
     258#endif 
    242259 
    243260         DO jk=1,jpkm1 
     
    278295 
    279296      IF((nbondi == 1).OR.(nbondi == 2)) THEN 
    280  
     297#if defined key_dynspg_flt 
    281298         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))) 
    288306 
    289307               ua(nlci-2:nlci-1,jj,jk) = ua(nlci-2:nlci-1,jj,jk) / fse3u(nlci-2:nlci-1,jj,jk) 
     
    292310         END DO 
    293311 
     312#if defined key_dynspg_flt 
    294313         DO jk=1,jpkm1 
    295314            DO jj=1,jpj 
     
    312331            ENDIF 
    313332         END DO 
     333#else 
     334         spgu(nlci-2,:) = ua_b(nlci-2,:) 
     335#endif 
    314336 
    315337         DO jk=1,jpkm1 
     
    353375      IF((nbondj == -1).OR.(nbondj == 2)) THEN 
    354376 
     377#if defined key_dynspg_flt 
    355378         DO ji=1,jpi 
    356379            laplacv(ji,2) = timeref * (zva2d(ji,2)/(zrhox*e1v(ji,2))) 
    357380         END DO 
     381#endif 
    358382 
    359383         DO jk=1,jpkm1 
     
    364388         END DO 
    365389 
     390#if defined key_dynspg_flt 
    366391         DO jk=1,jpkm1 
    367392            DO ji=1,jpi 
     
    383408            ENDIF 
    384409         END DO 
     410#else 
     411         spgv(:,2)=va_b(:,2) 
     412#endif 
    385413 
    386414         DO jk=1,jpkm1 
     
    413441         DO jk=1,jpkm1 
    414442            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)  
    416444               ua(ji,2,jk) = ua(ji,2,jk) / fse3u(ji,2,jk) 
    417445            END DO 
     
    422450      IF((nbondj == 1).OR.(nbondj == 2)) THEN 
    423451 
     452#if defined key_dynspg_flt 
    424453         DO ji=1,jpi 
    425454            laplacv(ji,nlcj-2) = timeref * (zva2d(ji,nlcj-2)/(zrhox*e1v(ji,nlcj-2))) 
    426455         END DO 
     456#endif 
    427457 
    428458         DO jk=1,jpkm1 
     
    433463         END DO 
    434464 
     465#if defined key_dynspg_flt 
    435466         DO jk=1,jpkm1 
    436467            DO ji=1,jpi 
     
    438469            END DO 
    439470         END DO 
    440  
    441471 
    442472         spgv(:,nlcj-2)=0. 
     
    453483            ENDIF 
    454484         END DO 
     485#else 
     486         spgv(:,nlcj-2)=va_b(:,nlcj-2) 
     487#endif 
    455488 
    456489         DO jk=1,jpkm1 
     
    483516         DO jk=1,jpkm1 
    484517            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) 
    486519               ua(ji,nlcj-1,jk) = ua(ji,nlcj-1,jk) / fse3u(ji,nlcj-1,jk) 
    487520            END DO 
     
    495528   END SUBROUTINE Agrif_dyn 
    496529 
     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 
    497654 
    498655   SUBROUTINE Agrif_ssh( kt ) 
     
    518675 
    519676      IF((nbondj == -1).OR.(nbondj == 2)) THEN 
    520          ssha(:,2)=sshn(:,3) 
    521          sshn(:,2)=sshb(:,3) 
     677         ssha(:,2)=ssha(:,3) 
     678         sshn(:,2)=sshn(:,3) 
    522679      ENDIF 
    523680 
    524681      IF((nbondj == 1).OR.(nbondj == 2)) THEN 
    525682         ssha(:,nlcj-1)=ssha(:,nlcj-2) 
    526          ssha(:,nlcj-1)=sshn(:,nlcj-2)                 
     683         sshn(:,nlcj-1)=sshn(:,nlcj-2)                 
    527684      ENDIF 
    528685 
    529686   END SUBROUTINE Agrif_ssh 
    530687 
     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 
    531727 
    532728   SUBROUTINE interpu(tabres,i1,i2,j1,j2,k1,k2) 
     
    611807 
    612808   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 
    613853 
    614854#else 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OFF_SRC/domain.F90

    r4248 r4292  
    295295      !!      vertical scale factors. 
    296296      !! 
    297       !! ** Method  : - reference 1D vertical coordinate (gdep._0, e3._0) 
     297      !! ** Method  : - reference 1D vertical coordinate (gdep._1d, e3._1d) 
    298298      !!              - read/set ocean depth and ocean levels (bathy, mbathy) 
    299299      !!              - vertical coordinate (gdep., e3.) depending on the  
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OFF_SRC/domrea.F90

    r3680 r4292  
    2525 
    2626   PUBLIC   dom_rea    ! routine called by inidom.F90 
     27  !! * Substitutions 
     28#  include "domzgr_substitute.h90" 
    2729   !!---------------------------------------------------------------------- 
    2830   !! NEMO/OFF 3.3 , NEMO Consortium (2010) 
     
    173175            CALL iom_get( inum4, jpdom_unknown, 'esigw', esigw ) 
    174176 
    175             CALL iom_get( inum4, jpdom_data, 'e3t', e3t ) ! scale factors 
    176             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 ) ! depth 
    181             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 ) 
    182184         ENDIF 
    183185 
    184186  
    185187         IF( ln_zps ) THEN                                           ! z-coordinate - partial steps 
    186             CALL iom_get( inum4, jpdom_unknown, 'gdept_0', gdept_0 )    ! reference depth 
    187             CALL iom_get( inum4, jpdom_unknown, 'gdepw_0', gdepw_0 ) 
    188             CALL iom_get( inum4, jpdom_unknown, 'e3t_0'  , e3t_0   )    ! reference scale factors 
    189             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   ) 
    190192            ! 
    191193            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(:,:,:) ) 
    196198            ELSE                                                        ! 2D bottom scale factors 
    197199               CALL iom_get( inum4, jpdom_data, 'e3t_ps', e3tp ) 
     
    199201               !                                                        ! deduces the 3D scale factors 
    200202               DO jk = 1, jpk 
    201                   e3t(:,:,jk) = e3t_0(jk)                                     ! set to the ref. factors 
    202                   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) 
    205207               END DO 
    206208               DO jj = 1,jpj                                                  ! adjust the deepest values 
    207209                  DO ji = 1,jpi 
    208210                     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) ) 
    211213                  END DO 
    212214               END DO 
     
    214216                  DO jj = 1, jpjm1 
    215217                     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) ) 
    218220                     END DO 
    219221                  END DO 
    220222               END DO 
    221                CALL lbc_lnk( e3u , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw, 'U', 1._wp )   ! lateral boundary conditions 
    222                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 ) 
    223225               ! 
    224226               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) 
    227229               END DO 
    228230            END IF 
    229231 
    230232            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(:,:,:) ) 
    233235            ELSE                                                           ! 2D bottom depth 
    234236               CALL iom_get( inum4, jpdom_data, 'hdept', zprt ) 
     
    236238               ! 
    237239               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) 
    240242               END DO 
    241243               DO jj = 1, jpj 
     
    243245                     ik = mbkt(ji,jj) 
    244246                     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) 
    248250                     ENDIF 
    249251                  END DO 
     
    254256 
    255257         IF( ln_zco ) THEN           ! Vertical coordinates and scales factors 
    256             CALL iom_get( inum4, jpdom_unknown, 'gdept_0', gdept_0 ) ! depth 
    257             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   ) 
    260262            DO jk = 1, jpk 
    261                e3t  (:,:,jk) = e3t_0(jk)                                     ! set to the ref. factors 
    262                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) 
    267269            END DO 
    268270         ENDIF 
     
    270272!!gm BUG in s-coordinate this does not work! 
    271273      ! 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 ~10m 
     274      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 
    274276      nla10 = nlb10 - 1                                              ! deepest    W level Above ~10m 
    275277!!gm end bug 
     
    312314         WRITE(numout,*) '              Reference z-coordinate depth and scale factors:' 
    313315         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 ) 
    315317      ENDIF 
    316318 
    317319      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 ' ) 
    320322      END DO 
    321323      !                                     ! ============================ 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r3651 r4292  
    88   !!            3.3  !  2010-09  (D. Storkey) add ice boundary conditions 
    99   !!            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 
    1011   !!---------------------------------------------------------------------- 
    1112#if defined key_bdy  
     
    2728      INTEGER, POINTER, DIMENSION(:,:)   ::  nbr 
    2829      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 
    3335   END TYPE OBC_INDEX 
    3436 
     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 
    3542   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 
    4358#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 
    4772#endif 
    4873   END TYPE OBC_DATA 
     
    6388   INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P  
    6489   !                                                        !  = 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 ;  
    6792                                                            !: = 1 read it in a NetCDF file 
    6893                                                            !: = 2 read tidal harmonic forcing from a NetCDF file 
    6994                                                            !: = 3 read external data AND tidal harmonic forcing from NetCDF files 
    70    INTEGER, DIMENSION(jp_bdy) ::   nn_dyn3d                 ! Choice of boundary condition for baroclinic velocities  
    71    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 ;  
    7297                                                            !: = 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 ;  
    75100                                                            !: = 1 read it in a NetCDF file 
    76101   LOGICAL, DIMENSION(jp_bdy) ::   ln_tra_dmp               !: =T Tracer damping 
    77102   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 
    79105 
    80106#if defined key_lim2 
    81    INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim2              ! Choice of boundary condition for sea ice variables  
    82    INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim2_dta          !: = 0 use the initial state as bdy dta ;  
    83                                                             !: = 1 read it in a NetCDF file 
     107   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 
    84110#endif 
    85111   ! 
     
    88114   !! Global variables 
    89115   !!---------------------------------------------------------------------- 
    90    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdytmask   !: Mask defining computational domain at T-points 
    91    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdyumask   !: Mask defining computational domain at U-points 
    92    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdyvmask   !: Mask defining computational domain at V-points 
     116   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 
    93119 
    94120   REAL(wp)                                    ::   bdysurftot !: Lateral surface of unstructured open boundary 
    95121 
    96    REAL(wp), POINTER, DIMENSION(:,:)           ::   pssh       !:  
    97    REAL(wp), POINTER, DIMENSION(:,:)           ::   phur       !:  
    98    REAL(wp), POINTER, DIMENSION(:,:)           ::   phvr       !: Pointers for barotropic fields  
    99    REAL(wp), POINTER, DIMENSION(:,:)           ::   pu2d       !:  
    100    REAL(wp), POINTER, DIMENSION(:,:)           ::   pv2d       !:  
     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   !:  
    101127 
    102128   !!---------------------------------------------------------------------- 
     
    109135   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2       !: workspace for reading in global data arrays (struct. bdy) 
    110136   TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process) 
    111    TYPE(OBC_DATA) , DIMENSION(jp_bdy)              ::   dta_bdy           !: bdy external data (local process) 
     137   TYPE(OBC_DATA) , DIMENSION(jp_bdy), TARGET      ::   dta_bdy           !: bdy external data (local process) 
    112138 
    113139   !!---------------------------------------------------------------------- 
     
    125151      !!---------------------------------------------------------------------- 
    126152      ! 
    127       ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj),                     
     153      ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj),      
    128154         &      STAT=bdy_oce_alloc ) 
    129155         ! 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_par.F90

    r3294 r4292  
    2323# endif 
    2424   INTEGER, PUBLIC, PARAMETER ::   jp_bdy  = 10       !: Maximum number of bdy sets 
    25    INTEGER, PUBLIC, PARAMETER ::   jpbtime = 1000     !: Max number of time dumps per file 
    2625   INTEGER, PUBLIC, PARAMETER ::   jpbgrd  = 3        !: Number of horizontal grid types used  (T, U, V) 
    2726 
    28    !! Flags for choice of schemes 
    29    INTEGER, PUBLIC, PARAMETER ::   jp_none         = 0        !: Flag for no open boundary condition 
    30    INTEGER, PUBLIC, PARAMETER ::   jp_frs          = 1        !: Flag for Flow Relaxation Scheme 
    31    INTEGER, PUBLIC, PARAMETER ::   jp_flather      = 2        !: Flag for Flather 
    3227#else 
    3328   !!---------------------------------------------------------------------- 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r4230 r4292  
    1111   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions 
    1212   !!            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 
    1314   !!---------------------------------------------------------------------- 
    1415#if defined key_bdy 
     
    2930   USE iom             ! IOM library 
    3031   USE in_out_manager  ! I/O logical units 
     32   USE dynspg_oce, ONLY: lk_dynspg_ts ! Split-explicit free surface flag 
    3133#if defined key_lim2 
    3234   USE ice_2 
     35#elif defined key_lim3 
     36   USE par_ice 
     37   USE ice 
     38   USE limcat_1D          ! redistribute ice input into categories 
    3339#endif 
    3440   USE sbcapr 
     
    4955 
    5056   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 
    5162 
    5263#  include "domzgr_substitute.h90" 
     
    7788                                                        ! etc. 
    7889      !! 
    79       INTEGER     ::  ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd  ! local indices 
     90      INTEGER     ::  ib_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd, jl  ! local indices 
    8091      INTEGER,          DIMENSION(jpbgrd) ::   ilen1  
    8192      INTEGER, POINTER, DIMENSION(:)      ::   nblen, nblenrim  ! short cuts 
     93      TYPE(OBC_DATA), POINTER             ::   dta              ! short cut 
    8294      !! 
    8395      !!--------------------------------------------------------------------------- 
     
    91103         ! Calculate depth-mean currents 
    92104         !----------------------------- 
    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 
    98109         DO ik = 1, jpkm1   !! Vertically integrated momentum trends 
    99              pu2d(:,:) = pu2d(:,:) + fse3u(:,:,ik) * umask(:,:,ik) * un(:,:,ik) 
    100              pv2d(:,:) = 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) 
    101112         END DO 
    102          pu2d(:,:) = pu2d(:,:) * hur(:,:) 
    103          pv2d(:,:) = pv2d(:,:) * hvr(:,:) 
     113         pun2d(:,:) = pun2d(:,:) * hur(:,:) 
     114         pvn2d(:,:) = pvn2d(:,:) * hvr(:,:) 
    104115          
    105116         DO ib_bdy = 1, nb_bdy 
     
    107118            nblen => idx_bdy(ib_bdy)%nblen 
    108119            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  
    111123               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) 
    137127                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    138128                     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 DO  
    142                igrd = 3  
    143                DO ib = 1, ilen1(igrd) 
    144                   DO ik = 1, jpkm1 
     129                     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) 
    145135                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    146136                     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) 
    157143                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    158144                     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 
    161237                  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 
    176259            ENDIF 
    177260#endif 
     
    179262         ENDDO ! ib_bdy 
    180263 
    181          CALL wrk_dealloc(jpi,jpj,pu2d,pv2d)  
     264         CALL wrk_dealloc(jpi,jpj,pun2d,pvn2d)  
    182265 
    183266      ENDIF ! kt .eq. nit000 
     
    188271      jstart = 1 
    189272      DO ib_bdy = 1, nb_bdy    
     273         dta => dta_bdy(ib_bdy) 
    190274         IF( nn_dta(ib_bdy) .eq. 1 ) THEN ! skip this bit if no external data required 
    191275       
     
    193277               ! Update barotropic boundary conditions only 
    194278               ! jit is optional argument for fld_read and bdytide_update 
    195                IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN 
     279               IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN 
    196280                  IF( nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
    197                      dta_bdy(ib_bdy)%ssh(:) = 0.0 
    198                      dta_bdy(ib_bdy)%u2d(:) = 0.0 
    199                      dta_bdy(ib_bdy)%v2d(:) = 0.0 
     281                     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 
    200284                  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 
    208289                        CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  & 
    209290                                     & 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 
    213293                        IF( ln_full_vel_array(ib_bdy) .AND.                                             & 
    214294                          &    ( nn_dyn2d_dta(ib_bdy) .EQ. 1 .OR. nn_dyn2d_dta(ib_bdy) .EQ. 3 .OR.  & 
     
    216296 
    217297                           igrd = 2                      ! zonal velocity 
    218                            dta_bdy(ib_bdy)%u2d(:) = 0.0 
     298                           dta%u2d(:) = 0.0 
    219299                           DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    220300                              ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    221301                              ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    222302                              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) 
    225305                              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) 
    230307                           END DO 
    231308                           igrd = 3                      ! meridional velocity 
    232                            dta_bdy(ib_bdy)%v2d(:) = 0.0 
     309                           dta%v2d(:) = 0.0 
    233310                           DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    234311                              ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    235312                              ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    236313                              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) 
    239316                              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) 
    244318                           END DO 
    245319                        ENDIF                     
    246320                     ENDIF 
    247321                     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),   &  
    249323                          &                 jit=jit, time_offset=time_offset ) 
    250324                     ENDIF 
     
    252326               ENDIF 
    253327            ELSE 
    254                IF (nn_tra(ib_bdy).eq.4) then      ! runoff condition 
     328               IF (cn_tra(ib_bdy) == 'runoff') then      ! runoff condition 
    255329                  jend = nb_bdy_fld(ib_bdy) 
    256330                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend),  & 
     
    261335                     ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    262336                     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) ) 
    264338                  END DO 
    265339                  ! 
     
    268342                     ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    269343                     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) ) 
    271345                  END DO 
    272346               ELSE 
    273                   IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
    274                      dta_bdy(ib_bdy)%ssh(:) = 0.0 
    275                      dta_bdy(ib_bdy)%u2d(:) = 0.0 
    276                      dta_bdy(ib_bdy)%v2d(:) = 0.0 
     347                  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 
    277351                  ENDIF 
    278                   IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data 
    279                      jend = nb_bdy_fld(ib_bdy) 
     352                  IF( dta%nread(1) .gt. 0 ) THEN ! update external data 
     353                     jend = jstart + dta%nread(1) - 1 
    280354                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), & 
    281355                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset ) 
     
    286360                    &   nn_dyn3d_dta(ib_bdy) .EQ. 1 ) ) THEN 
    287361                     igrd = 2                      ! zonal velocity 
    288                      dta_bdy(ib_bdy)%u2d(:) = 0.0 
     362                     dta%u2d(:) = 0.0 
    289363                     DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    290364                        ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    291365                        ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    292366                        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) 
    295369                        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) 
    297371                        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) 
    299373                        END DO 
    300374                     END DO 
    301375                     igrd = 3                      ! meridional velocity 
    302                      dta_bdy(ib_bdy)%v2d(:) = 0.0 
     376                     dta%v2d(:) = 0.0 
    303377                     DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    304378                        ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    305379                        ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    306380                        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) 
    309383                        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) 
    311385                        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) 
    313387                        END DO 
    314388                     END DO 
    315389                  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) 
    323400         END IF ! nn_dta(ib_bdy) = 1 
    324401      END DO  ! ib_bdy 
    325402 
     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 
    326410      IF ( ln_apr_obc ) THEN 
    327411         DO ib_bdy = 1, nb_bdy 
    328             IF (nn_tra(ib_bdy).NE.4)THEN 
     412            IF (cn_tra(ib_bdy) /= 'runoff')THEN 
    329413               igrd = 1                      ! meridional velocity 
    330414               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     
    349433      !!                for open boundary conditions 
    350434      !! 
    351       !! ** Method  :   Use fldread.F90 
     435      !! ** Method  :    
    352436      !!                 
    353437      !!---------------------------------------------------------------------- 
     
    362446                                                                ! =F => baroclinic velocities in 3D boundary data 
    363447      INTEGER                                ::   ilen_global   ! Max length required for global bdy dta arrays 
    364       INTEGER,              DIMENSION(jpbgrd) ::  ilen0         ! size of local arrays 
    365448      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ilen1, ilen3  ! size of 1st and 3rd dimensions of local arrays 
    366449      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ibdy           ! bdy set for a particular jfld 
    367450      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   igrid         ! index for grid type (1,2,3 = T,U,V) 
    368451      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 
    369458      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   blf_i         !  array of namelist information structures 
    370459      TYPE(FLD_N) ::   bn_tem, bn_sal, bn_u3d, bn_v3d   !  
     
    372461#if defined key_lim2 
    373462      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       
    374465#endif 
    375466      NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d  
    376467#if defined key_lim2 
    377468      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 
    378471#endif 
    379472      NAMELIST/nambdy_dta/ ln_full_vel 
     
    392485                               ,nn_dyn3d_dta(ib_bdy)       & 
    393486                               ,nn_tra_dta(ib_bdy)         & 
    394 #if defined key_lim2 
    395                                ,nn_ice_lim2_dta(ib_bdy)    & 
     487#if ( defined key_lim2 || defined key_lim3 ) 
     488                              ,nn_ice_lim_dta(ib_bdy)    & 
    396489#endif 
    397490                              ) 
     
    404497      nb_bdy_fld(:) = 0 
    405498      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 ) ) THEN 
     499         IF( cn_dyn2d(ib_bdy) /= 'none' .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN 
    407500            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3 
    408501         ENDIF 
    409          IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) THEN 
     502         IF( cn_dyn3d(ib_bdy) /= 'none' .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) THEN 
    410503            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2 
    411504         ENDIF 
    412          IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 1  ) THEN 
     505         IF( cn_tra(ib_bdy) /= 'none' .and. nn_tra_dta(ib_bdy) .eq. 1  ) THEN 
    413506            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2 
    414507         ENDIF 
    415 #if defined key_lim2 
    416          IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1  ) THEN 
     508#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 
    417510            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3 
    418511         ENDIF 
     
    458551            nblen => idx_bdy(ib_bdy)%nblen 
    459552            nblenrim => idx_bdy(ib_bdy)%nblenrim 
     553            dta => dta_bdy(ib_bdy) 
     554            dta%nread(2) = 0 
    460555 
    461556            ! Only read in necessary fields for this set. 
    462557            ! 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' 
    466562                  jfld = jfld + 1 
    467563                  blf_i(jfld) = bn_ssh 
     
    470566                  ilen1(jfld) = nblen(igrid(jfld)) 
    471567                  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' 
    475573                  jfld = jfld + 1 
    476574                  blf_i(jfld) = bn_u2d 
     
    479577                  ilen1(jfld) = nblen(igrid(jfld)) 
    480578                  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' 
    482584                  jfld = jfld + 1 
    483585                  blf_i(jfld) = bn_v2d 
     
    486588                  ilen1(jfld) = nblen(igrid(jfld)) 
    487589                  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 
    510621 
    511622            ENDIF 
    512623 
    513624            ! 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 
    529646 
    530647            ENDIF 
     
    532649#if defined key_lim2 
    533650            ! 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 
    558726#endif 
    559727            ! Recalculate field counts 
     
    568736            ENDIF 
    569737 
     738            dta%nread(1) = nb_bdy_fld(ib_bdy) 
     739 
    570740         ENDIF ! nn_dta .eq. 1 
    571741      ENDDO ! ib_bdy 
     
    596766 
    597767         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)) ) 
    608796               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' 
    616823               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' 
    618828               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' 
    643840               jfld = jfld + 1 
    644841               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' 
    645845               jfld = jfld + 1 
    646846               dta_bdy(ib_bdy)%sal => bf(jfld)%fnow(:,1,:) 
     
    649849 
    650850#if defined key_lim2 
    651          IF (nn_ice_lim2(ib_bdy) .gt. 0) THEN 
     851         IF (nn_ice_lim(ib_bdy) .gt. 0) THEN 
    652852            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)) ) 
    657856            ELSE 
    658857               jfld = jfld + 1 
     
    662861               jfld = jfld + 1 
    663862               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 
    664892            ENDIF 
    665893         ENDIF 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90

    r4153 r4292  
    3030   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    3131   USE in_out_manager  ! 
    32    USE domvvl          ! variable volume 
     32   USE domvvl 
    3333 
    3434   IMPLICIT NONE 
     
    5757      LOGICAL, INTENT( in ), OPTIONAL :: dyn3d_only       ! T => only update baroclinic velocities 
    5858      !! 
    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 
    6263 
    6364      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn') 
     
    7071      ENDIF 
    7172 
     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 
    7279      !------------------------------------------------------- 
    7380      ! Set pointers 
     
    7784      phur => hur 
    7885      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)  
    8088 
    8189      !------------------------------------------------------- 
     
    8391      !------------------------------------------------------- 
    8492 
    85       pu2d(:,:) = 0.e0 
    86       pv2d(:,:) = 0.e0 
     93      ! "After" velocities:  
     94 
     95      pua2d(:,:) = 0.e0 
     96      pva2d(:,:) = 0.e0 
     97       
    8798      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(:,:) 
    94111      ELSE 
    95          DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
    96             pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
    97             pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
    98          END DO 
    99          pu2d(:,:) = pu2d(:,:) * phur(:,:) 
    100          pv2d(:,:) = 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(:,:) 
    101118      ENDIF 
     119 
    102120      DO jk = 1 , jpkm1 
    103          ua(:,:,jk) = ua(:,:,jk) - pu2d(:,:) * umask(:,:,jk) 
    104          va(:,:,jk) = va(:,:,jk) - pv2d(:,:) * vmask(:,:,jk) 
     121         ua(:,:,jk) = ua(:,:,jk) - pua2d(:,:) 
     122         va(:,:,jk) = va(:,:,jk) - pva2d(:,:) 
    105123      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 
    106158 
    107159      !------------------------------------------------------- 
     
    119171 
    120172      DO jk = 1 , jpkm1 
    121          ua(:,:,jk) = ( ua(:,:,jk) + pu2d(:,:) ) * umask(:,:,jk) 
    122          va(:,:,jk) = ( va(:,:,jk) + pv2d(:,:) ) * vmask(:,:,jk) 
     173         ua(:,:,jk) = ( ua(:,:,jk) + pua2d(:,:) ) * umask(:,:,jk) 
     174         va(:,:,jk) = ( va(:,:,jk) + pva2d(:,:) ) * vmask(:,:,jk) 
    123175      END DO 
    124176 
    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)  
    126186 
    127187      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn') 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90

    r3680 r4292  
    66   !! History :  3.4  !  2011     (D. Storkey) new module as part of BDY rewrite 
    77   !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Optimization of BDY communications 
     8   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_bdy  
     
    1112   !!   'key_bdy' :                    Unstructured Open Boundary Condition 
    1213   !!---------------------------------------------------------------------- 
    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 
    1519   !!---------------------------------------------------------------------- 
    1620   USE timing          ! Timing 
     
    1822   USE dom_oce         ! ocean space and time domain 
    1923   USE bdy_oce         ! ocean open boundary conditions 
     24   USE bdylib          ! BDY library routines 
    2025   USE dynspg_oce      ! for barotropic variables 
    2126   USE phycst          ! physical constants 
     
    2631   PRIVATE 
    2732 
    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 
    2935 
    3036   !!---------------------------------------------------------------------- 
     
    4854      DO ib_bdy=1, nb_bdy 
    4955 
    50          SELECT CASE( nn_dyn2d(ib_bdy) ) 
    51          CASE(jp_none) 
     56         SELECT CASE( cn_dyn2d(ib_bdy) ) 
     57         CASE('none') 
    5258            CYCLE 
    53          CASE(jp_frs) 
     59         CASE('frs') 
    5460            CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy ) 
    55          CASE(jp_flather) 
     61         CASE('flather') 
    5662            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. ) 
    5767         CASE DEFAULT 
    5868            CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' ) 
     
    8999         ij   = idx%nbj(jb,igrd) 
    90100         zwgt = idx%nbw(jb,igrd) 
    91          pu2d(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) 
    92102      END DO 
    93103      ! 
     
    97107         ij   = idx%nbj(jb,igrd) 
    98108         zwgt = idx%nbw(jb,igrd) 
    99          pv2d(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) 
    100110      END DO  
    101       CALL lbc_bdy_lnk( pu2d, 'U', -1., ib_bdy )  
    102       CALL lbc_bdy_lnk( pv2d, 'V', -1., ib_bdy)   ! Boundary points should be updated 
     111      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )  
     112      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy)   ! Boundary points should be updated 
    103113      ! 
    104114      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_frs') 
     
    133143      INTEGER  ::   jb, igrd                         ! dummy loop indices 
    134144      INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses 
     145      REAL(wp), POINTER :: flagu, flagv              ! short cuts 
    135146      REAL(wp) ::   zcorr                            ! Flather correction 
    136147      REAL(wp) ::   zforc                            ! temporary scalar 
     148      REAL(wp) ::   zflag, z1_2                      !    "        " 
    137149      !!---------------------------------------------------------------------- 
    138150 
    139151      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_fla') 
     152 
     153      z1_2 = 0.5_wp 
    140154 
    141155      ! ---------------------------------! 
     
    160174         ii  = idx%nbi(jb,igrd) 
    161175         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  
    164179         ! 
    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)  
    168187      END DO 
    169188      ! 
     
    173192         ii  = idx%nbi(jb,igrd) 
    174193         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  
    177197         ! 
    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 )   ! 
    184208      ! 
    185209      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_fla') 
    186210      ! 
    187211   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 
    188300#else 
    189301   !!---------------------------------------------------------------------- 
     
    192304CONTAINS 
    193305   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 
    195308   END SUBROUTINE bdy_dyn2d 
     309 
    196310#endif 
    197311 
    198312   !!====================================================================== 
    199313END MODULE bdydyn2d 
     314 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90

    r3703 r4292  
    1919   USE dom_oce         ! ocean space and time domain 
    2020   USE bdy_oce         ! ocean open boundary conditions 
     21   USE bdylib          ! for orlanski library routines 
    2122   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2223   USE in_out_manager  ! 
     
    5253      DO ib_bdy=1, nb_bdy 
    5354 
    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') 
    6057            CYCLE 
    61          CASE(jp_frs) 
     58         CASE('frs') 
    6259            CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
    63          CASE(2) 
     60         CASE('specified') 
    6461            CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
    65          CASE(3) 
     62         CASE('zero') 
    6663            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. ) 
    6768         CASE DEFAULT 
    6869            CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' ) 
     
    109110         END DO 
    110111      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 )    
    112114      ! 
    113115      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
     
    204206         END DO 
    205207      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 )    
    207210      ! 
    208211      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
     
    211214 
    212215   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 
    213254 
    214255   SUBROUTINE bdy_dyn3d_dmp( kt ) 
     
    225266      REAL(wp) ::   zwgt           ! boundary weight 
    226267      INTEGER  ::  ib_bdy          ! loop index 
     268      REAL(wp), POINTER, DIMENSION(:,:) :: phur1, phvr1     ! inverse depth at u and v points 
    227269      !!---------------------------------------------------------------------- 
    228270      ! 
     
    232274      ! Remove barotropic part from before velocity 
    233275      !------------------------------------------------------- 
    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. 
    239283      DO jk = 1, jpkm1 
    240284#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) 
    243289#else 
    244          pu2d(:,:) = pu2d(:,:) + fse3u_0(:,:,jk) * ub(:,:,jk)  * umask(:,:,jk) 
    245          pv2d(:,:) = 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) 
    246292#endif 
    247293      END DO 
    248294 
    249295      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(:,:) 
    252300      ELSE 
    253          pu2d(:,:) = pv2d(:,:) * hur(:,:) 
    254          pv2d(:,:) = pu2d(:,:) * hvr(:,:) 
     301         pub2d(:,:) = pvb2d(:,:) * hur(:,:) 
     302         pvb2d(:,:) = pub2d(:,:) * hvr(:,:) 
    255303      ENDIF 
    256304 
    257305      DO ib_bdy=1, nb_bdy 
    258          IF ( ln_dyn3d_dmp(ib_bdy).and.nn_dyn3d(ib_bdy).gt.0 ) THEN 
     306         IF ( ln_dyn3d_dmp(ib_bdy) .and. cn_dyn3d(ib_bdy) /= 'none' ) THEN 
    259307            igrd = 2                      ! Relaxation of zonal velocity 
    260308            DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     
    264312               DO jk = 1, jpkm1 
    265313                  ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - & 
    266                                    ub(ii,ij,jk) + pu2d(ii,ij)) ) * umask(ii,ij,jk) 
     314                                   ub(ii,ij,jk) + pub2d(ii,ij)) ) * umask(ii,ij,jk) 
    267315               END DO 
    268316            END DO 
     
    275323               DO jk = 1, jpkm1 
    276324                  va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) -  & 
    277                                    vb(ii,ij,jk) + pv2d(ii,ij)) ) * vmask(ii,ij,jk) 
     325                                   vb(ii,ij,jk) + pvb2d(ii,ij)) ) * vmask(ii,ij,jk) 
    278326               END DO 
    279327            END DO 
     
    281329      ENDDO 
    282330      ! 
    283       CALL wrk_dealloc(jpi,jpj,pu2d,pv2d)  
     331      CALL wrk_dealloc(jpi,jpj,pub2d,pvb2d)  
    284332      ! 
    285333      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  
    2121   !!   bdy_init       : Initialization of unstructured open boundaries 
    2222   !!---------------------------------------------------------------------- 
     23   USE wrk_nemo        ! Memory Allocation 
    2324   USE timing          ! Timing 
    2425   USE oce             ! ocean dynamics and tracers variables 
     
    7980      INTEGER  ::   jpbdtau, jpbdtas                       !   -       - 
    8081      INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       - 
     82      INTEGER  ::   i_offset, j_offset                     !   -       - 
    8183      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 
    8386      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars 
    8487      INTEGER, DIMENSION (2)                  ::   kdimsz 
     
    9093      INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b  ! Flags for boundaries receiving 
    9194      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) 
    9296 
    9397      !! 
    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_lim2 
    99          &             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,                           & 
    100104#endif 
    101105         &             ln_vol, nn_volctl, nn_rimwidth 
     
    156160 
    157161        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' ) 
    163189        END SELECT 
    164         IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN 
     190        IF( cn_dyn2d(ib_bdy) /= 'none' ) THEN 
    165191           SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !  
    166192              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
     
    177203 
    178204        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' ) 
    185231        END SELECT 
    186         IF( nn_dyn3d(ib_bdy) .gt. 0 ) THEN 
     232        IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 
    187233           SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   !  
    188234              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
     
    193239 
    194240        IF ( ln_dyn3d_dmp(ib_bdy) ) THEN 
    195            IF ( nn_dyn3d(ib_bdy).EQ.0 ) THEN 
     241           IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN 
    196242              IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.' 
    197243              ln_dyn3d_dmp(ib_bdy)=.false. 
    198            ELSEIF ( nn_dyn3d(ib_bdy).EQ.1 ) THEN 
     244           ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN 
    199245              CALL ctl_stop( 'Use FRS OR relaxation' ) 
    200246           ELSE 
     
    202248              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
    203249              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. 
    204252           ENDIF 
    205253        ELSE 
     
    209257 
    210258        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' ) 
    218289        END SELECT 
    219         IF( nn_tra(ib_bdy) .gt. 0 ) THEN 
     290        IF( cn_tra(ib_bdy) /= 'none' ) THEN 
    220291           SELECT CASE( nn_tra_dta(ib_bdy) )                   !  
    221292              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
     
    226297 
    227298        IF ( ln_tra_dmp(ib_bdy) ) THEN 
    228            IF ( nn_tra(ib_bdy).EQ.0 ) THEN 
     299           IF ( cn_tra(ib_bdy) == 'none' ) THEN 
    229300              IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.' 
    230301              ln_tra_dmp(ib_bdy)=.false. 
    231            ELSEIF ( nn_tra(ib_bdy).EQ.1 ) THEN 
     302           ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN 
    232303              CALL ctl_stop( 'Use FRS OR relaxation' ) 
    233304           ELSE 
    234305              IF(lwp) WRITE(numout,*) '      + T/S relaxation zone' 
    235306              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' 
    236308              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. 
    237311           ENDIF 
    238312        ELSE 
     
    243317#if defined key_lim2 
    244318        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' ) 
    249331        END SELECT 
    250         IF( nn_ice_lim2(ib_bdy) .gt. 0 ) THEN  
    251            SELECT CASE( nn_ice_lim2_dta(ib_bdy) )                   !  
     332        IF( cn_ice_lim(ib_bdy) /= 'none' ) THEN  
     333           SELECT CASE( nn_ice_lim_dta(ib_bdy) )                   !  
    252334              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'         
    253335              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' ) 
    255360           END SELECT 
    256361        ENDIF 
     
    740845               IF(lwp) THEN         ! Since all procs read global data only need to do this check on one proc... 
    741846                  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') 
    746849                  ENDIF     
    747850               ENDIF 
     
    766869         ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) ) 
    767870         ALLOCATE( idx_bdy(ib_bdy)%nbd(ilen1,jpbgrd) ) 
     871         ALLOCATE( idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) ) 
    768872         ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) ) 
    769873         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) ) 
    772876 
    773877         ! Dispatch mapping indices and discrete distances on each processor 
     
    9371041            ENDDO 
    9381042         ENDDO  
     1043 
    9391044         ! definition of the i- and j- direction local boundaries arrays 
    9401045         ! used for sending the boudaries 
     
    9901095               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 
    9911096               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 ) &  
    9921099               & *(FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.   ! quadratic 
    9931100            END DO 
     
    10921199      ENDDO 
    10931200 
     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 
    10941211      ! Lateral boundary conditions 
     1212      CALL lbc_lnk( zfmask       , 'F', 1. ) 
    10951213      CALL lbc_lnk( fmask        , 'F', 1. )   ;   CALL lbc_lnk( bdytmask(:,:), 'T', 1. ) 
    10961214      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. ) 
     
    10981216      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components 
    10991217 
    1100          idx_bdy(ib_bdy)%flagu(:) = 0.e0 
    1101          idx_bdy(ib_bdy)%flagv(:) = 0.e0 
     1218         idx_bdy(ib_bdy)%flagu(:,:) = 0.e0 
     1219         idx_bdy(ib_bdy)%flagv(:,:) = 0.e0 
    11021220         icount = 0  
    11031221 
    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 
    11071226   
    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  
    11191261         END DO 
    11201262 
    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  
    11361302         END DO 
    11371303 
    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 
    11481305 
    11491306      ! Compute total lateral surface for volume correction: 
     
    11571314               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
    11581315               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) 
    11601317               bdysurftot = bdysurftot + hu     (nbi  , nbj)                           & 
    11611318                  &                    * e2u    (nbi  , nbj) * ABS( flagu ) & 
     
    11701327               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
    11711328               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) 
    11731330               bdysurftot = bdysurftot + hv     (nbi, nbj  )                           & 
    11741331                  &                    * e1v    (nbi, nbj  ) * ABS( flagv ) & 
     
    11861343         DEALLOCATE(nbidta, nbjdta, nbrdta) 
    11871344      ENDIF 
     1345 
     1346      CALL wrk_dealloc(jpi,jpj,zfmask)  
    11881347 
    11891348      IF( nn_timing == 1 ) CALL timing_stop('bdy_init') 
     
    15801739      itest = 0 
    15811740 
    1582       IF (nn_dyn2d(ib1)/=nn_dyn2d(ib2)) itest = itest + 1 
    1583       IF (nn_dyn3d(ib1)/=nn_dyn3d(ib2)) itest = itest + 1 
    1584       IF (nn_tra(ib1)/=nn_tra(ib2)) itest = itest + 1 
     1741      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 
    15851744      ! 
    15861745      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  
    99   !!            3.3  !  2010-09  (D.Storkey and E.O'Dea)  bug fixes 
    1010   !!            3.4  !  2012-09  (G. Reffray and J. Chanut) New inputs + mods 
     11   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes 
    1112   !!---------------------------------------------------------------------- 
    1213#if defined key_bdy 
     
    3233!   USE tide_mod       ! Useless ?? 
    3334   USE fldread, ONLY: fld_map 
     35   USE dynspg_oce, ONLY: lk_dynspg_ts 
    3436 
    3537   IMPLICIT NONE 
     
    3840   PUBLIC   bdytide_init     ! routine called in bdy_init 
    3941   PUBLIC   bdytide_update   ! routine called in bdy_dta 
     42   PUBLIC   bdy_dta_tides    ! routine called in dyn_spg_ts 
    4043 
    4144   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data 
     
    4952 
    5053   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) 
    5155 
    5256   !!---------------------------------------------------------------------- 
     
    131135            ! JC: If FRS scheme is used, we assume that tidal is needed over the whole 
    132136            ! relaxation area       
    133             IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     137            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN 
    134138               ilen0(:)=nblen(:) 
    135139            ELSE 
     
    146150            ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) ) 
    147151 
    148             td%ssh0(:,:,:) = 0.e0 
    149             td%ssh(:,:,:) = 0.e0 
    150             td%u0(:,:,:) = 0.e0 
    151             td%u(:,:,:) = 0.e0 
    152             td%v0(:,:,:) = 0.e0 
    153             td%v(:,:,:) = 0.e0 
     152            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 
    154158 
    155159            IF (ln_bdytide_2ddta) THEN 
     
    255259            ENDIF 
    256260            ! 
     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            ! 
    257271         ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2 
    258272         ! 
     
    300314      ENDIF 
    301315 
    302       IF ( nsec_day == NINT(0.5 * rdttra(1)) .AND. zflag==1 ) THEN 
     316      IF ( nsec_day == NINT(0.5_wp * rdttra(1)) .AND. zflag==1 ) THEN 
    303317        ! 
    304318        kt_tide = kt 
     
    321335          
    322336      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) ) 
    324338      ELSE                               
    325339         z_arg = ((kt-kt_tide)+time_add) * rdt 
     
    327341 
    328342      ! 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) 
    331345 
    332346      DO itide = 1, nb_harmo 
     
    354368      ! 
    355369   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 
    356498 
    357499   SUBROUTINE tide_init_elevation( idx, td ) 
     
    460602      WRITE(*,*) 'bdytide_update: You should not have seen this print! error?', kt, jit 
    461603   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 
    462610#endif 
    463611 
    464612   !!====================================================================== 
    465613END MODULE bdytides 
     614 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90

    r3777 r4292  
    2020   USE dom_oce         ! ocean space and time domain variables  
    2121   USE bdy_oce         ! ocean open boundary conditions 
     22   USE bdylib          ! for orlanski library routines 
    2223   USE bdydta, ONLY:   bf 
    2324   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     
    5152      DO ib_bdy=1, nb_bdy 
    5253 
    53          SELECT CASE( nn_tra(ib_bdy) ) 
    54          CASE(jp_none) 
     54         SELECT CASE( cn_tra(ib_bdy) ) 
     55         CASE('none') 
    5556            CYCLE 
    56          CASE(jp_frs) 
     57         CASE('frs') 
    5758            CALL bdy_tra_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
    58          CASE(2) 
     59         CASE('specified') 
    5960            CALL bdy_tra_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
    60          CASE(3) 
     61         CASE('neumann') 
    6162            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') 
    6368            CALL bdy_tra_rnf( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
    6469         CASE DEFAULT 
     
    196201      ! 
    197202   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 
    198237 
    199238   SUBROUTINE bdy_tra_rnf( idx, dta, kt ) 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90

    r3294 r4292  
    104104               ii = idx%nbi(jb,jgrd) 
    105105               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) 
    107107            END DO 
    108108         END DO 
     
    112112               ii = idx%nbi(jb,jgrd) 
    113113               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)  
    115115            END DO 
    116116         END DO 
     
    136136               ii = idx%nbi(jb,jgrd) 
    137137               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) 
    140140            END DO 
    141141         END DO 
     
    145145               ii = idx%nbi(jb,jgrd) 
    146146               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) 
    149149            END DO 
    150150         END DO 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

    r3294 r4292  
    196196      thick0(:,:) = 0._wp 
    197197      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) 
    200200      END DO 
    201201      IF( lk_mpp )   CALL mpp_sum( vol0 ) 
     
    212212               ik = mbkt(ji,jj) 
    213213               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) ) 
    215215                  sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1) 
    216216               ENDIF 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diadimg.F90

    r3294 r4292  
    112112 
    113113    CASE ( 'T') 
    114        z4dep(:)=gdept_0(:) 
     114       z4dep(:)=gdept_1d(:) 
    115115 
    116116    CASE ( 'W' ) 
    117        z4dep(:)=gdepw_0(:) 
     117       z4dep(:)=gdepw_1d(:) 
    118118 
    119119    CASE ( '2' ) 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90

    r4147 r4292  
    11MODULE diaharm  
    2  
    3 #if defined key_diaharm && defined key_tide 
    4    !!================================================================================= 
     2   !!====================================================================== 
    53   !!                       ***  MODULE  diaharm  *** 
    64   !! 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   !!---------------------------------------------------------------------- 
    913   USE oce             ! ocean dynamics and tracers variables 
    1014   USE dom_oce         ! ocean space and time domain 
    11    USE in_out_manager  ! I/O units 
    12    USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    13    USE ioipsl          ! NetCDF IPSL library 
    14    USE diadimg         ! To write dimg 
    1515   USE phycst 
    1616   USE dynspg_oce 
     
    1818   USE daymod 
    1919   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 
    2125   USE timing          ! preformance summary 
    2226   USE wrk_nemo        ! working arrays 
     
    3034   INTEGER, PARAMETER :: jpdimsparse  = jpincomax*300*24 
    3135 
    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   !!---------------------------------------------------------------------- 
    6464CONTAINS 
    6565 
     
    6767      !!---------------------------------------------------------------------- 
    6868      !!                 ***  ROUTINE dia_harm_init  *** 
    69       !!---------------------------------------------------------------------- 
    7069      !!          
    7170      !! ** Purpose :   Initialization of tidal harmonic analysis 
     
    7372      !! ** Method  :   Initialize frequency array and  nodal factor for nit000_han 
    7473      !! 
    75       !! History : 
    76       !!   9.0  O. Le Galloudec and J. Chanut (Original) 
    77       !!-------------------------------------------------------------------- 
    78       !! * Local declarations  
     74      !!-------------------------------------------------------------------- 
    7975      INTEGER :: jh, nhan, jk, ji 
    8076      INTEGER ::   ios                 ! Local integer output status for namelist read 
     
    108104      ! Basic checks on harmonic analysis time window: 
    109105      ! ---------------------------------------------- 
    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 
    127115      DO jk=1,jpmax_harmo 
    128116         DO ji=1,jpmax_harmo 
     
    157145      ! Initialize frequency array: 
    158146      ! --------------------------- 
    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 ) 
    165150 
    166151      IF(lwp) WRITE(numout,*) 'Analysed frequency  : ',nb_ana ,'Frequency ' 
     
    172157      ! Initialize temporary arrays: 
    173158      ! ---------------------------- 
    174       ALLOCATE( ana_temp(jpi,jpj,nb_ana*2,3)) 
     159      ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) ) 
    175160      ana_temp(:,:,:,:) = 0.e0 
    176161 
    177162   END SUBROUTINE dia_harm_init 
    178     
     163 
     164 
    179165   SUBROUTINE dia_harm ( kt ) 
    180166      !!---------------------------------------------------------------------- 
    181167      !!                 ***  ROUTINE dia_harm  *** 
    182       !!---------------------------------------------------------------------- 
    183168      !!          
    184169      !! ** Purpose :   Tidal harmonic analysis main routine 
     
    186171      !! ** Action  :   Sums ssh/u/v over time analysis [nit000_han,nitend_han] 
    187172      !! 
    188       !! History : 
    189       !!   9.0  O. Le Galloudec and J. Chanut (Original) 
    190       !!-------------------------------------------------------------------- 
    191       !! * Argument: 
     173      !!-------------------------------------------------------------------- 
    192174      INTEGER, INTENT( IN ) :: kt 
    193  
    194       !! * Local declarations 
     175      ! 
    195176      INTEGER  :: ji, jj, jh, jc, nhc 
    196177      REAL(wp) :: ztime, ztemp 
     
    198179      IF( nn_timing == 1 )   CALL timing_start('dia_harm') 
    199180 
    200       IF ( kt .EQ. nit000 ) CALL dia_harm_init 
     181      IF ( kt == nit000 ) CALL dia_harm_init 
    201182 
    202183      IF ( ((kt.GE.nit000_han).AND.(kt.LE.nitend_han)).AND. & 
     
    215196              DO ji = 1,jpi 
    216197                ! 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)         
    219199#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) 
    226202#endif 
    227203              END DO 
     
    233209      END IF 
    234210 
    235       IF ( kt .EQ. nitend_han ) CALL dia_harm_end 
     211      IF ( kt == nitend_han )  CALL dia_harm_end 
    236212 
    237213      IF( nn_timing == 1 )   CALL timing_stop('dia_harm') 
     
    239215   END SUBROUTINE dia_harm 
    240216 
     217 
    241218   SUBROUTINE dia_harm_end 
    242219      !!---------------------------------------------------------------------- 
    243220      !!                 ***  ROUTINE diaharm_end  *** 
    244       !!---------------------------------------------------------------------- 
    245221      !!          
    246222      !! ** Purpose :  Compute the Real and Imaginary part of tidal constituents 
     
    248224      !! ** Action  :  Decompose the signal on the harmonic constituents  
    249225      !! 
    250       !! History : 
    251       !!   9.0  O. Le Galloudec and J. Chanut (Original) 
    252       !!-------------------------------------------------------------------- 
    253  
    254       !! * Local declarations 
     226      !!-------------------------------------------------------------------- 
    255227      INTEGER :: ji, jj, jh, jc, jn, nhan, jl 
    256228      INTEGER :: ksp, kun, keq 
     
    283255               nisparse(ksp) = keq 
    284256               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)) ) 
    288259            END DO 
    289260         END DO 
    290261      END DO 
    291262 
    292       nsparse=ksp 
     263      nsparse = ksp 
    293264 
    294265      ! Elevation: 
     
    296267         DO ji = 1, jpi 
    297268            ! Fill input array 
    298             kun=0 
    299             DO jh = 1,nb_ana 
    300                DO jc = 1,2 
     269            kun = 0 
     270            DO jh = 1, nb_ana 
     271               DO jc = 1, 2 
    301272                  kun = kun + 1 
    302273                  ztmp4(kun)=ana_temp(ji,jj,kun,1) 
    303                ENDDO 
    304             ENDDO 
     274               END DO 
     275            END DO 
    305276 
    306277            CALL SUR_DETERMINE(jj) 
     
    314285      END DO 
    315286 
    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)  ) 
    319290 
    320291      DO jj = 1, jpj 
    321292         DO ji = 1, jpi 
    322293            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) 
    327298            ENDDO 
    328299         ENDDO 
     
    402373   END SUBROUTINE dia_harm_end 
    403374 
     375 
    404376   SUBROUTINE dia_wri_harm 
    405377      !!-------------------------------------------------------------------- 
    406378      !!                 ***  ROUTINE dia_wri_harm  *** 
    407       !!-------------------------------------------------------------------- 
    408379      !!          
    409380      !! ** 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      !!-------------------------------------------------------------------- 
    417382      CHARACTER(LEN=lc) :: cltext 
    418383      CHARACTER(LEN=lc) ::   & 
     
    472437#else 
    473438      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) ) 
    476441      END DO 
    477442#endif 
    478443 
    479444   END SUBROUTINE dia_wri_harm 
     445 
    480446 
    481447   SUBROUTINE SUR_DETERMINE(init) 
     
    486452   !!        
    487453   !!--------------------------------------------------------------------------------- 
    488    INTEGER, INTENT(in) :: init  
    489                 
     454   INTEGER, INTENT(in) ::   init  
     455   ! 
    490456   INTEGER                         :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd 
    491457   REAL(wp)                        :: zval1, zval2, zx1 
     
    496462   CALL wrk_alloc( jpincomax , ipos2 , ipivot        ) 
    497463             
    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      ! 
    508470      DO jk1_sd = 1, nsparse 
    509471         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) 
    514474            IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN 
    515475               ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd))  & 
    516476                                                        + valuesparse(jk1_sd)*valuesparse(jk2_sd) 
    517477            ENDIF 
    518  
    519          ENDDO 
    520       ENDDO 
     478         END DO 
     479      END DO 
    521480 
    522481      DO jj_sd = 1 ,ninco 
     
    588547   ENDDO 
    589548 
    590  
    591549   CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 ) 
    592550   CALL wrk_dealloc( jpincomax , ipos2 , ipivot        ) 
     
    594552  END SUBROUTINE SUR_DETERMINE 
    595553 
    596  
    597554#else 
    598555   !!---------------------------------------------------------------------- 
     
    601558   LOGICAL, PUBLIC, PARAMETER ::   lk_diaharm = .FALSE. 
    602559CONTAINS 
    603  
    604560   SUBROUTINE dia_harm ( kt )     ! Empty routine 
    605561      INTEGER, INTENT( IN ) :: kt   
    606562      WRITE(*,*) 'dia_harm: you should not have seen this print' 
    607563   END SUBROUTINE dia_harm 
    608  
    609  
    610 #endif 
     564#endif 
     565 
    611566   !!====================================================================== 
    612567END MODULE diaharm 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90

    r3764 r4292  
    304304      ! ----------------------------- ! 
    305305 
    306       ! find ilevel with (ilevel+1) the deepest W-level above 300m (we assume we can use e3t_0 to 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...) 
    307307      ilevel   = 0 
    308308      zthick_0 = 0._wp 
    309309      DO jk = 1, jpkm1                       
    310          zthick_0 = zthick_0 + e3t_0(jk) 
     310         zthick_0 = zthick_0 + e3t_1d(jk) 
    311311         IF( zthick_0 < 300. )   ilevel = jk 
    312312      END DO 
     
    326326            htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem) * MIN( fse3t(ji,jj,ilevel+1), zthick(ji,jj) )  & 
    327327                                                                   * 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) 
    330328         END DO 
    331329      END DO 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90

    r4148 r4292  
    259259      ! 
    260260#if defined key_mpp_mpi 
     261      ijpjjpk = jpj*jpk 
    261262      ish(1) = ijpjjpk  ;   ish2(1) = jpj   ;   ish2(2) = jpk 
    262263      zwork(1:ijpjjpk) = RESHAPE( p_fval, ish ) 
     
    314315      END DO 
    315316#if defined key_mpp_mpi 
     317      ijpjjpk = jpj*jpk 
    316318      ish(1) = jpj*jpk   ;   ish2(1) = jpj   ;   ish2(2) = jpk 
    317319      zwork(1:ijpjjpk)= RESHAPE( p_fval, ish ) 
     
    670672            CALL histbeg(clhstnam, 1, zfoo, jpj, zphi,   & 
    671673               1, 1, 1, jpj, niter, zjulian, zdt*nn_fptr, nhoridz, numptr, domain_id=nidom_ptr) 
    672             ! Vertical grids : gdept_0, gdepw_0 
     674            ! Vertical grids : gdept_1d, gdepw_1d 
    673675            CALL histvert( numptr, "deptht", "Vertical T levels",   & 
    674                &                   "m", jpk, gdept_0, ndepidzt, "down" ) 
     676               &                   "m", jpk, gdept_1d, ndepidzt, "down" ) 
    675677            CALL histvert( numptr, "depthw", "Vertical W levels",   & 
    676                &                   "m", jpk, gdepw_0, ndepidzw, "down" ) 
     678               &                   "m", jpk, gdepw_1d, ndepidzw, "down" ) 
    677679            ! 
    678680            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  
    2525   USE oce             ! ocean dynamics and tracers  
    2626   USE dom_oce         ! ocean space and time domain 
     27   USE dynadv, ONLY: ln_dynadv_vec 
    2728   USE zdf_oce         ! ocean vertical physics 
    2829   USE ldftra_oce      ! ocean active tracers: lateral physics 
     
    4445   USE diadimg         ! dimg direct access file format output 
    4546   USE diaar5, ONLY :   lk_diaar5 
     47   USE dynadv, ONLY :   ln_dynadv_vec 
    4648   USE iom 
    4749   USE ioipsl 
     
    144146      ENDIF 
    145147 
    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. 
    159176      IF( lk_zdfddm ) THEN 
    160177         CALL iom_put( "avs" , fsavs(:,:,:)                          )    ! S vert. eddy diff. coef. 
     
    252269      ! 
    253270      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 ) 
    255272      ! 
    256273      ! Output the initial state and forcings 
     
    325342            &          nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set ) 
    326343         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" ) 
    328345         !                                                            ! Index of ocean points 
    329346         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume 
     
    361378            &          nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set ) 
    362379         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" ) 
    364381         !                                                            ! Index of ocean points 
    365382         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume 
     
    374391            &          nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set ) 
    375392         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" ) 
    377394         !                                                            ! Index of ocean points 
    378395         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume 
     
    387404            &          nit000-1, zjulian, zdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set ) 
    388405         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" ) 
    390407 
    391408 
     
    397414         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn 
    398415            &          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 
    399424         !                                                                                      !!! nid_T : 2D 
    400425         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst 
     
    408433         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx 
    409434            &          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) 
    412437            &                                                                  , "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) 
    415440            &                                                                  , "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 
    418443         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr 
    419444            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     
    587612      ! --------------------- 
    588613 
    589       ! ndex(1) est utilise ssi l'avant dernier argument est diffferent de  
     614      ! ndex(1) est utilise ssi l'avant dernier argument est different de  
    590615      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument 
    591616      ! donne le nombre d'elements, et ndex la liste des indices a sortir 
     
    597622 
    598623      ! 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 
    603642      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height 
    604643      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux 
     
    606645                                                                                  ! (includes virtual salt flux beneath ice  
    607646                                                                                  ! 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 sst 
    611       zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal) 
    612       CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )             ! c/d term on sss 
    613 #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 
    614653      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux 
    615654      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux 
     
    752791      ! 
    753792      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 ) 
    755794      ! 
    756795      IF( nn_timing == 1 )   CALL timing_stop('dia_wri') 
     
    813852          1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit 
    814853      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") 
    816855 
    817856      ! Declare all the output fields as NetCDF variables 
     
    841880      CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress 
    842881         &          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 
    843886 
    844887#if defined key_lim2 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r4247 r4292  
    152152   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  gphit, gphiu   !: latitude  of t-, u-, v- and f-points (degre) 
    153153   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) 
    158158   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  e1e2t          !: surface at t-point (m2) 
    159159   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ff             !: coriolis factor (2.*omega*sin(yphi) ) (s-1) 
     
    169169   !! All coordinates 
    170170   !! --------------- 
    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-W  points (m) 
    173    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3v   , e3f     !: analytical vertical scale factors at  V--F 
    174    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e3t   , e3u     !:                                       T--U  points (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--UW points (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) 
    177177#if defined key_vvl 
    178178   LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .TRUE.    !: variable grid flag 
     
    180180   !! All coordinates 
    181181   !! --------------- 
    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) 
    190193#else 
    191194   LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .FALSE.   !: fixed grid flag 
    192195#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 
    196206 
    197207   INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1) 
     
    200210   !! z-coordinate with full steps (also used in the other cases as reference z-coordinate) 
    201211   !! =-----------------====------ 
    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 points 
     212   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 
    205215 
    206216   !! s-coordinate and hybrid z-s-coordinate 
    207217   !! =----------------======--------------- 
    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-levels 
    211  
    212    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbatv , hbatf    !: ocean depth at the vertical of  V--F 
    213    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbatt , hbatu    !:                                 T--U points (m) 
    214    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   scosrf, scobot   !: ocean surface and bottom topographies  
    215    !                                        !  (if deviating from coordinate surfaces in HYBRID) 
    216    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hifv  , hiff     !: interface depth between stretching at  V--F 
    217    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hift  , hifu     !: and quasi-uniform spacing              T--U points (m) 
    218    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rx1              !: Maximum grid stiffness ratio 
     218   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 
    219229 
    220230   !!---------------------------------------------------------------------- 
    221231   !! masks, bathymetry 
    222232   !! --------------------------------------------------------------------- 
    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 level 
    225    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbku, mbkv   !: vertical index of the bottom last U- and W- ocean level 
    226    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy        !: ocean depth (meters) 
    227    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i      !: interior domain T-point mask 
    228    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask        !: land/ocean mask of barotropic stream function 
    229  
    230    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts 
     233   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 
    231241 
    232242   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tpol, fpol          !: north fold mask (jperio= 3 or 4) 
    233243 
    234244#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      !: ??? 
    237247#endif 
    238248 
     
    316326         &      glamf(jpi,jpj) , gphif(jpi,jpj) , e1f(jpi,jpj) , e2f(jpi,jpj) , ff   (jpi,jpj) , STAT=ierr(3) )      
    317327         ! 
    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) ) 
    321331         ! 
    322332#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) ) 
    336353         ! 
    337354      ALLOCATE( hbatv (jpi,jpj) , hbatf (jpi,jpj) ,     & 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r4245 r4292  
    8787                             CALL dom_msk      ! Masks 
    8888      IF( ln_sco )           CALL dom_stiff    ! Maximum stiffness ratio/hydrostatic consistency 
    89       IF( lk_vvl         )   CALL dom_vvl      ! Vertical variable mesh 
     89      IF( lk_vvl )           CALL dom_vvl_init ! Vertical variable mesh 
    9090      ! 
    9191      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(:,:) 
    92105      ! 
    93106      hu(:,:) = 0._wp                          ! Ocean depth at U- and V-points 
    94107      hv(:,:) = 0._wp 
    95108      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) 
    98111      END DO 
    99112      !                                        ! Inverse of the local depth 
     
    407420         DO jj = 2, jpjm1 
    408421            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) ) 
    425438               zrxmax = MAXVAL(zr1(1:4)) 
    426439               rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax) 
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domstp.F90

    r2715 r4292  
    9090 
    9191            DO jk = 1, jpk 
    92                IF( gdept_0(jk) <= rdth ) rdttra(jk) = rdtmin 
    93                IF( gdept_0(jk) >  rdth ) THEN 
     92               IF( gdept_1d(jk) <= rdth ) rdttra(jk) = rdtmin 
     93               IF( gdept_1d(jk) >  rdth ) THEN 
    9494                  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. ) 
    9797               ENDIF 
    9898               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  
    66   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code 
    77   !!            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 
    1010   !!---------------------------------------------------------------------- 
    1111   !!   'key_vvl'                              variable volume 
    1212   !!---------------------------------------------------------------------- 
    13    !!   dom_vvl     : defined coefficients to distribute ssh on each layers 
    1413   !!---------------------------------------------------------------------- 
     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 
    1524   USE oce             ! ocean dynamics and tracers 
    1625   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 
    1927   USE in_out_manager  ! I/O manager 
     28   USE iom             ! I/O manager library 
     29   USE restart         ! ocean restart 
    2030   USE lib_mpp         ! distributed memory computing library 
    2131   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
     
    2636   PRIVATE 
    2737 
    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 
    3666 
    3767   !! * Substitutions 
     
    3969#  include "vectopt_loop_substitute.h90" 
    4070   !!---------------------------------------------------------------------- 
    41    !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
     71   !! NEMO/OPA 3.3 , NEMO-Consortium (2010)  
    4272   !! $Id$ 
    4373   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4474   !!---------------------------------------------------------------------- 
    45 CONTAINS        
     75 
     76CONTAINS 
    4677 
    4778   INTEGER FUNCTION dom_vvl_alloc() 
    4879      !!---------------------------------------------------------------------- 
    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(:,:,:), 'U' ) 
     153      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 
     154  &