New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 6736 – NEMO

Changeset 6736


Ignore:
Timestamp:
2016-06-24T09:50:27+02:00 (8 years ago)
Author:
jamesharle
Message:

FASTNEt code modifications

Location:
branches/NERC/dev_r3874_FASTNEt/NEMOGCM
Files:
6 added
142 edited

Legend:

Unmodified
Added
Removed

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/CONFIG/cfg.txt

    r3769 r6736  
    99ORCA2_LIM_CFC_C14b OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
    1010ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 
     11NNA_R12 OPA_SRC LIM_SRC_2  

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90

    r3625 r6736  
    1919   PUBLIC    ice_alloc_2  !  Called in iceini_2.F90 
    2020 
    21    INTEGER , PUBLIC ::   numit        !: ice iteration index 
    22    REAL(wp), PUBLIC ::   rdt_ice      !: ice time step 
     21   INTEGER , PUBLIC ::   numit     !: ice iteration index 
     22   REAL(wp), PUBLIC ::   rdt_ice   !: ice time step 
    2323 
    2424   !                                                                     !!* namicerun read in iceini  * 
    … …  
    2727   LOGICAL               , PUBLIC ::   ln_limdyn     = .TRUE.             !: flag for ice dynamics (T) or not (F) 
    2828   LOGICAL               , PUBLIC ::   ln_limdmp     = .FALSE.            !: Ice damping 
     29   LOGICAL               , PUBLIC ::   ln_vp2evp     = .FALSE.            !: restart from a vp file in an evp run  
    2930   LOGICAL               , PUBLIC ::   ln_nicep      = .TRUE.             !: flag grid points output (T) or not (F) 
    3031   REAL(wp)              , PUBLIC ::   hsndif        = 0._wp              !: snow temp. computation (0) or not (9999) 
    … …  
    9899   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qstoif        !: Energy stored in the brine pockets 
    99100   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fbif          !: Heat flux at the ice base 
    100    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdm_snw       !: Variation of snow mass over 1 time step           [Kg/m2] 
    101    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdq_snw       !: Heat content associated with rdm_snw              [J/m2] 
    102    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdm_ice       !: Variation of ice  mass over 1 time step           [Kg/m2] 
    103    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdq_ice       !: Heat content associated with rdm_ice              [J/m2] 
     101   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdmsnif       !: Variation of snow mass 
     102   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rdmicif       !: Variation of ice mass 
    104103   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qldif         !: heat balance of the lead (or of the open ocean) 
    105104   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qcmif         !: Energy needed to freeze the ocean surface layer 
    … …  
    155154 
    156155      ALLOCATE(phicif(jpi,jpj) , pfrld  (jpi,jpj) , qstoif (jpi,jpj) ,     & 
    157          &     fbif  (jpi,jpj) , rdm_snw(jpi,jpj) , rdq_snw(jpi,jpj) ,     & 
    158          &                       rdm_ice(jpi,jpj) , rdq_ice(jpi,jpj) ,     & 
     156         &     fbif  (jpi,jpj) , rdmsnif(jpi,jpj) , rdmicif(jpi,jpj) ,     & 
    159157         &     qldif (jpi,jpj) , qcmif  (jpi,jpj) , fdtcn  (jpi,jpj) ,     & 
    160158         &     qdtcn (jpi,jpj) , thcm   (jpi,jpj)                    , STAT=ierr(4) ) 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90

    r3625 r6736  
    1313   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    1414   !!---------------------------------------------------------------------- 
    15    !!   ice_init_2    : sea-ice model initialization 
    16    !!   ice_run_2     : Definition some run parameter for ice model 
     15   !!   ice_init_2       : sea-ice model initialization 
     16   !!   ice_run_2        : Definition some run parameter for ice model 
    1717   !!---------------------------------------------------------------------- 
    18    USE phycst         ! physical constants 
    19    USE dom_oce        ! ocean domain 
    20    USE sbc_oce        ! surface boundary condition: ocean 
    21    USE sbc_ice        ! LIM2 surface boundary condition 
    22    USE dom_ice_2      ! LIM2 ice domain 
    23    USE par_ice_2      ! LIM2 parameters 
    24    USE thd_ice_2      ! LIM2 thermodynamical variables 
    25    USE ice_2          ! LIM2 ice variable 
    26    USE limmsh_2       ! LIM2 mesh 
    27    USE limistate_2    ! LIM2 initial state 
    28    USE limrst_2       ! LIM2 restart 
    29    USE limsbc_2       ! LIM2 surface boundary condition 
    30    USE in_out_manager ! I/O manager 
    31    USE lib_mpp        ! MPP library 
    32    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     18   USE phycst           ! physical constants 
     19   USE dom_oce          ! ocean domain 
     20   USE sbc_oce          ! surface boundary condition: ocean 
     21   USE sbc_ice          ! LIM2 surface boundary condition 
     22   USE dom_ice_2        ! LIM2 ice domain 
     23   USE par_ice_2        ! LIM2 parameters 
     24   USE thd_ice_2        ! LIM2 thermodynamical variables 
     25   USE ice_2            ! LIM2 ice variable 
     26   USE limmsh_2         ! LIM2 mesh 
     27   USE limistate_2      ! LIM2 initial state 
     28   USE limrst_2         ! LIM2 restart 
     29   USE limsbc_2         ! LIM2 surface boundary condition 
     30   USE in_out_manager   ! I/O manager 
     31   USE lib_mpp          ! MPP library 
    3332 
    3433   IMPLICIT NONE 
    … …  
    110109      !! ** input   :   Namelist namicerun 
    111110      !!------------------------------------------------------------------- 
    112       NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, ln_limdmp, acrit, hsndif, hicdif 
     111      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, ln_limdmp, acrit, hsndif, hicdif, ln_vp2evp 
    113112      !!------------------------------------------------------------------- 
    114113      !                     
    … …  
    125124         WRITE(numout,*) '   computation of temp. in snow (=0) or not (=9999) hsndif = ', hsndif 
    126125         WRITE(numout,*) '   computation of temp. in ice  (=0) or not (=9999) hicdif = ', hicdif 
     126         WRITE(numout,*) '   Restart EVP run from VP restart file (set stresses to 0)= ', ln_vp2evp 
    127127      ENDIF 
    128128      ! 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/limadv_2.F90

    r3625 r6736  
    1414   !!   'key_lim2'                                    LIM 2.0 sea-ice model 
    1515   !!---------------------------------------------------------------------- 
    16    !!   lim_adv_x_2   : advection of sea ice on x axis 
    17    !!   lim_adv_y_2   : advection of sea ice on y axis 
     16   !!   lim_adv_x_2  : advection of sea ice on x axis 
     17   !!   lim_adv_y_2  : advection of sea ice on y axis 
    1818   !!---------------------------------------------------------------------- 
    1919   USE dom_oce 
    … …  
    2121   USE ice_2 
    2222   USE lbclnk 
    23    USE in_out_manager ! I/O manager 
    24    USE lib_mpp        ! MPP library 
    25    USE wrk_nemo       ! work arrays 
    26    USE prtctl         ! Print control 
    27    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     23   USE in_out_manager     ! I/O manager 
     24   USE lib_mpp            ! MPP library 
     25   USE wrk_nemo           ! work arrays 
     26   USE prtctl             ! Print control 
     27   USE lib_fortran        ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     28 
    2829 
    2930   IMPLICIT NONE 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/limdia_2.F90

    r3625 r6736  
    2424   USE in_out_manager  ! I/O manager 
    2525   USE lib_mpp         ! MPP library 
    26    USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2726 
    2827   IMPLICIT NONE 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/limdmp_2.F90

    r3635 r6736  
    1919   USE in_out_manager ! I/O manager 
    2020   USE lib_mpp        ! MPP library 
    21    USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2221 
    2322   IMPLICIT NONE 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90

    r3625 r6736  
    3131   USE in_out_manager   ! I/O manager 
    3232   USE prtctl           ! Print control 
    33    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3433 
    3534   IMPLICIT NONE 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90

    r3625 r6736  
    2121   USE prtctl           ! Print control 
    2222   USE in_out_manager   ! I/O manager 
    23    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2423 
    2524   IMPLICIT NONE 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/limistate_2.F90

    r3625 r6736  
    2727   USE iom 
    2828   USE in_out_manager 
    29    USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3029 
    3130   IMPLICIT NONE 
    … …  
    193192            IF(lwp) WRITE(numout,*) '                  ice state initialization with : Ice_initialization.nc' 
    194193             
     194#if defined key_lim2_initcd_alt1 
     195            CALL iom_get( inum_ice, jpdom_data, 'hicif', hicif )       
     196            CALL iom_get( inum_ice, jpdom_data, 'hsnif', hsnif )       
     197            CALL iom_get( inum_ice, jpdom_data, 'frld' , frld  )      
     198            CALL iom_get( inum_ice, jpdom_data, 'tbif1'   , tbif(:,:,1)  ) 
     199            CALL iom_get( inum_ice, jpdom_data, 'tbif2'   , tbif(:,:,2)  ) 
     200            CALL iom_get( inum_ice, jpdom_data, 'tbif3'   , tbif(:,:,3)  ) 
     201            CALL iom_get( inum_ice, jpdom_data, 'sist'   , sist  ) 
     202#elif defined key_lim2_initcd_alt2 
     203            CALL iom_get( inum_ice, jpdom_data, 'iicethic', hicif )       
     204            CALL iom_get( inum_ice, jpdom_data, 'isnowthi', hsnif )       
     205            CALL iom_get( inum_ice, jpdom_data, 'ileadfra' , frld  )      
     206            CALL iom_get( inum_ice, jpdom_data, 'isstempe'   , sist  ) 
     207            CALL iom_get( inum_ice, jpdom_unknown, 'iicetemp', tbif(1:nlci,1:nlcj,:),   & 
     208                 &        kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,jplayersp1 /) ) 
     209#else 
    195210            CALL iom_get( inum_ice, jpdom_data, 'hicif', hicif )       
    196211            CALL iom_get( inum_ice, jpdom_data, 'hsnif', hsnif )       
    … …  
    199214            CALL iom_get( inum_ice, jpdom_unknown, 'tbif', tbif(1:nlci,1:nlcj,:),   & 
    200215                 &        kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,jplayersp1 /) ) 
     216#endif          
    201217            ! put some values in the extra-halo... 
    202218            DO jj = nlcj+1, jpj   ;   tbif(1:nlci,jj,:) = tbif(1:nlci,nlej,:)   ;   END DO 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90

    r3625 r6736  
    2323   USE wrk_nemo         ! work arrays 
    2424#endif 
    25    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2625 
    2726   IMPLICIT NONE 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90

    r3680 r6736  
    3030   USE in_out_manager ! I/O manager 
    3131   USE prtctl         ! Print control 
    32    USE oce     , ONLY : snwice_mass, snwice_mass_b 
    33    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    34 #if defined key_agrif 
    35    USE agrif_lim2_interp ! nesting 
    36 #endif 
     32   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    3733 
    3834   IMPLICIT NONE 
    … …  
    8581      REAL(wp) ::   zs21_11, zs21_12, zs21_21, zs21_22 
    8682      REAL(wp) ::   zs22_11, zs22_12, zs22_21, zs22_22 
    87       REAL(wp) ::   zintb, zintn 
    8883      REAL(wp), POINTER, DIMENSION(:,:) ::   zfrld, zmass, zcorl 
    8984      REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct, za2ct, zresr 
    9085      REAL(wp), POINTER, DIMENSION(:,:) ::   zc1u, zc1v, zc2u, zc2v 
    91       REAL(wp), POINTER, DIMENSION(:,:) ::   zsang, zpice 
     86      REAL(wp), POINTER, DIMENSION(:,:) ::   zsang 
    9287      REAL(wp), POINTER, DIMENSION(:,:) ::   zu0, zv0 
    9388      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_n, zv_n 
    … …  
    9994       
    10095      CALL wrk_alloc( jpi,jpj, zfrld, zmass, zcorl, za1ct, za2ct, zresr ) 
    101       CALL wrk_alloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang, zpice ) 
     96      CALL wrk_alloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang ) 
    10297      CALL wrk_alloc( jpi,jpj+2, zu0, zv0, zu_n, zv_n, zu_a, zv_a, zviszeta, zviseta, kjstart = 0 ) 
    10398      CALL wrk_alloc( jpi,jpj+2, zzfrld, zztms, zi1, zi2, zmasst, zpresh, kjstart = 0 ) 
    … …  
    135130!i    zviszeta(:,jpj+1) = 0._wp    ;    zviseta(:,jpj+1) = 0._wp 
    136131 
    137       IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==! 
    138           ! 
    139           ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 
    140           !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1} 
    141          zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp 
    142           ! 
    143           ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 
    144           !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 
    145          zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 
    146           ! 
    147          zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0 
    148           ! 
    149          ! 
    150       ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==! 
    151          zpice(:,:) = ssh_m(:,:) 
    152       ENDIF 
    153 #if defined key_agrif 
    154       ! load the boundary value of velocity in special array zuive and zvice 
    155       CALL agrif_rhg_lim2_load 
    156 #endif 
    157132 
    158133      ! Ice mass, ice strength, and wind stress at the center            | 
    … …  
    222197 
    223198            ! Gradient of the sea surface height 
    224             zgsshx =  (   (zpice(ji  ,jj  ) - zpice(ji-1,jj  ))/e1u(ji-1,jj  )   & 
    225                &       +  (zpice(ji  ,jj-1) - zpice(ji-1,jj-1))/e1u(ji-1,jj-1)   ) * 0.5_wp 
    226             zgsshy =  (   (zpice(ji  ,jj  ) - zpice(ji  ,jj-1))/e2v(ji  ,jj-1)   & 
    227                &       +  (zpice(ji-1,jj  ) - zpice(ji-1,jj-1))/e2v(ji-1,jj-1)   ) * 0.5_wp 
     199            zgsshx =  (   (ssh_m(ji  ,jj  ) - ssh_m(ji-1,jj  ))/e1u(ji-1,jj  )   & 
     200               &       +  (ssh_m(ji  ,jj-1) - ssh_m(ji-1,jj-1))/e1u(ji-1,jj-1)   ) * 0.5_wp 
     201            zgsshy =  (   (ssh_m(ji  ,jj  ) - ssh_m(ji  ,jj-1))/e2v(ji  ,jj-1)   & 
     202               &       +  (ssh_m(ji-1,jj  ) - ssh_m(ji-1,jj-1))/e2v(ji-1,jj-1)   ) * 0.5_wp 
    228203 
    229204            ! Computation of the velocity field taking into account the ice-ice interaction.                                  
    … …  
    559534            CALL lbc_lnk( zv_n(:,1:jpj), 'I', -1. ) 
    560535 
    561 #if defined key_agrif 
    562             ! copy the boundary value from u_ice_nst and v_ice_nst to u_ice and v_ice 
    563             ! before next interations 
    564             CALL agrif_rhg_lim2(zu_n,zv_n) 
    565 #endif 
    566  
    567536            ! Test of Convergence 
    568537            DO jj = k_j1+1 , k_jpj-1 
    … …  
    607576 
    608577      CALL wrk_dealloc( jpi,jpj, zfrld, zmass, zcorl, za1ct, za2ct, zresr ) 
    609       CALL wrk_dealloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang, zpice ) 
     578      CALL wrk_dealloc( jpi,jpj, zc1u , zc1v , zc2u , zc2v , zsang ) 
    610579      CALL wrk_dealloc( jpi,jpj+2, zu0, zv0, zu_n, zv_n, zu_a, zv_a, zviszeta, zviseta, kjstart = 0 ) 
    611580      CALL wrk_dealloc( jpi,jpj+2, zzfrld, zztms, zi1, zi2, zmasst, zpresh, kjstart = 0 ) 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/limrst_2.F90

    r2528 r6736  
    225225      CALL iom_get( numrir, jpdom_autoglo, 'fsbbq'      , fsbbq  )     
    226226#if ! defined key_lim2_vp 
     227      IF ( ln_vp2evp ) THEN 
     228      stress1_i (:,:) = 0._wp                          ! EVP rheology 
     229      stress2_i (:,:) = 0._wp 
     230      stress12_i(:,:) = 0._wp 
     231      ELSE 
    227232      CALL iom_get( numrir, jpdom_autoglo, 'stress1_i'  , stress1_i  ) 
    228233      CALL iom_get( numrir, jpdom_autoglo, 'stress2_i'  , stress2_i  ) 
    229234      CALL iom_get( numrir, jpdom_autoglo, 'stress12_i' , stress12_i ) 
     235      ENDIF 
    230236#endif 
    231237      CALL iom_get( numrir, jpdom_autoglo, 'sxice'      , sxice  ) 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90

    r3625 r6736  
    99   !!            3.3  ! 2009-05 (G. Garric, C. Bricaud) addition of the lim2_evp case 
    1010   !!             -   ! 2010-11 (G. Madec) ice-ocean stress computed at each ocean time-step 
    11    !!           3.3.1 ! 2011-01 (A. R. Porter, STFC Daresbury) dynamical allocation 
    12    !!            3.5  ! 2012-11 ((G. Madec, Y. Aksenov, A. Coward) salt and heat fluxes associated with e-p 
     11   !!            4.0  ! 2011-01  (A. R. Porter, STFC Daresbury) dynamical allocation 
    1312   !!---------------------------------------------------------------------- 
    1413#if defined key_lim2 
    … …  
    2928   USE sbc_oce          ! surface boundary condition: ocean 
    3029   USE sbccpl 
    31    USE cpl_oasis3, ONLY : lk_cpl 
    32    USE oce       , ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass  
     30 
    3331   USE albedo           ! albedo parameters 
    3432   USE lbclnk           ! ocean lateral boundary condition - MPP exchanges 
    … …  
    3937   USE iom              ! I/O library 
    4038   USE prtctl           ! Print control 
    41    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     39   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
     40   USE cpl_oasis3, ONLY : lk_cpl 
    4241 
    4342   IMPLICIT NONE 
    … …  
    9089      !!              - Update the fluxes provided to the ocean 
    9190      !!      
    92       !! ** Outputs : - qsr     : sea heat flux    : solar  
    93       !!              - qns     : sea heat flux    : non solar (including heat content of the mass flux) 
    94       !!              - emp     : freshwater budget: mass flux  
    95       !!              - sfx     : freshwater budget: salt flux due to Freezing/Melting 
     91      !! ** Outputs : - qsr     : sea heat flux:    solar  
     92      !!              - qns     : sea heat flux: non solar 
     93      !!              - emp     : freshwater budget: volume flux  
     94      !!              - emps    : freshwater budget: concentration/dillution  
    9695      !!              - utau    : sea surface i-stress (ocean referential) 
    9796      !!              - vtau    : sea surface j-stress (ocean referential) 
    … …  
    109108      INTEGER  ::   ifvt, i1mfr, idfr, iflt    !   -       - 
    110109      INTEGER  ::   ial, iadv, ifral, ifrdv    !   -       - 
    111       REAL(wp) ::   zqsr,     zqns,   zfmm     ! local scalars 
    112       REAL(wp) ::   zinda,    zfsalt, zemp     !   -      - 
    113       REAL(wp) ::   zemp_snw, zqhc,   zcd      !   -      - 
    114       REAL(wp) ::   zswitch                    !   -      - 
     110      REAL(wp) ::   zqsr, zqns, zfm            ! local scalars 
     111      REAL(wp) ::   zinda, zfons, zemp         !   -      - 
    115112      REAL(wp), POINTER, DIMENSION(:,:)   ::   zqnsoce       ! 2D workspace 
    116113      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp   ! 2D/3D workspace 
    … …  
    119116      CALL wrk_alloc( jpi, jpj, zqnsoce ) 
    120117      CALL wrk_alloc( jpi, jpj, 1, zalb, zalbp ) 
    121  
    122       SELECT CASE( nn_ice_embd )                 ! levitating or embedded sea-ice option 
    123         CASE( 0    )   ;   zswitch = 1           ! (0) standard levitating sea-ice : salt exchange only 
    124         CASE( 1, 2 )   ;   zswitch = 0           ! (1) levitating sea-ice: salt and volume exchange but no pressure effect 
    125                                                  ! (2) embedded sea-ice : salt and volume fluxes and pressure 
    126       END SELECT                                 !     
    127118 
    128119      !------------------------------------------! 
    … …  
    143134            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv  
    144135 
    145 !!$            attempt to explain the tricky flags set above.... 
    146 !!$            zinda   = 1.0 - AINT( pfrld(ji,jj) )                   ! = 0. if ice-free ocean else 1. (after ice adv, but before ice thermo) 
    147 !!$            i1mfr   = 1.0 - AINT(  frld(ji,jj) )                   ! = 0. if ice-free ocean else 1. (after ice thermo) 
    148 !!$ 
    149 !!$            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      ! = zinda if previous thermodynamic step overmelted the ice??? 
    150 !!$            ELSE                             ;   ifvt = 0.         !  
     136!!$            zinda   = 1.0 - AINT( pfrld(ji,jj) )                   !   = 0. if pure ocean else 1. (at previous time) 
     137!!$ 
     138!!$            i1mfr   = 1.0 - AINT(  frld(ji,jj) )                   !   = 0. if pure ocean else 1. (at current  time) 
     139!!$ 
     140!!$            IF( phicif(ji,jj) <= 0. ) THEN   ;   ifvt = zinda      !   = 1. if (snow and no ice at previous time) else 0. ??? 
     141!!$            ELSE                             ;   ifvt = 0. 
    151142!!$            ENDIF 
    152143!!$ 
    153 !!$            IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases due to ice thermodynamics 
     144!!$            IF( frld(ji,jj) >= pfrld(ji,jj) ) THEN   ;   idfr = 0.  !   = 0. if lead fraction increases from previous to current 
    154145!!$            ELSE                                     ;   idfr = 1.    
    155146!!$            ENDIF 
    156147!!$ 
    157 !!$            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )    !   = 1. if ice (not only snow) at previous time and ice-free ocean currently 
     148!!$            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )    !   = 1. if ice (not only snow) at previous and pure ocean at current 
    158149!!$ 
    159150!!$            ial     = ifvt   * i1mfr    +    ( 1 - ifvt ) * idfr 
    160 !!$                    = i1mfr if ifvt = 1 i.e.  
    161 !!$                    = idfr  if ifvt = 0 
    162151!!$!                 snow no ice   ice         ice or nothing  lead fraction increases 
    163152!!$!                 at previous   now           at previous 
    164 !!$!                -> ice area increases  ???         -> ice area decreases ??? 
     153!!$!                -> ice aera increases  ???         -> ice aera decreases ??? 
    165154!!$ 
    166155!!$            iadv    = ( 1  - i1mfr ) * zinda 
    … …  
    186175#endif             
    187176            !  computation the non solar heat flux at ocean surface 
    188             zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr                                              &   ! part of the solar energy used in leads 
    189                &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                             & 
    190                &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice  & 
    191                &       + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdtice  
    192  
    193             fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! store residual heat flux (to put into the ocean at the next time-step) 
    194             zqhc = ( rdq_snw(ji,jj)                                     & 
    195                  & + rdq_ice(ji,jj) * ( 1.- zswitch) ) * r1_rdtice       ! heat flux due to snow ( & ice heat content,  
    196             !                                                            !           if ice/ocean mass exchange active)  
     177            zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr   &   ! part of the solar energy used in leads 
     178               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                            & 
     179               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice    & 
     180               &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) )                   * r1_rdtice  
     181 
     182            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ??? 
     183            ! 
    197184            qsr  (ji,jj) = zqsr                                          ! solar heat flux  
    198             qns  (ji,jj) = zqns - fdtcn(ji,jj) + zqhc                    ! non solar heat flux  
    199             ! 
    200             !                          !------------------------------------------! 
    201             !                          !  mass and salt flux at the ocean surface ! 
    202             !                          !------------------------------------------! 
    203             ! 
    204             ! mass flux at the ocean-atmosphere interface (open ocean fraction = leads area) 
    205 #if defined key_coupled 
    206             !                                                  ! coupled mode:  
    207             zemp = + emp_tot(ji,jj)                            &     ! net mass flux over the grid cell (ice+ocean area) 
    208                &   - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )          ! minus the mass flux intercepted by sea-ice 
    209 #else 
    210             !                                                  ! forced  mode:  
    211             zemp = + emp(ji,jj)     *         frld(ji,jj)      &     ! mass flux over open ocean fraction  
    212                &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &     ! liquid precip. over ice reaches directly the ocean 
    213                &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )          ! snow is intercepted by sea-ice (previous frld) 
    214 #endif             
    215             ! 
    216             ! mass flux at the ocean/ice interface (sea ice fraction) 
    217             zemp_snw = rdm_snw(ji,jj) * r1_rdtice                    ! snow melting = pure water that enters the ocean 
    218             zfmm     = rdm_ice(ji,jj) * r1_rdtice                    ! Freezing minus Melting (F-M) 
    219  
    220             ! salt flux at the ice/ocean interface (sea ice fraction) [PSU*kg/m2/s] 
    221             zfsalt = - sice_0(ji,jj) * zfmm                          ! F-M salt exchange 
    222             zcd    =   soce_0(ji,jj) * zfmm                          ! concentration/dilution term due to F-M 
    223             ! 
    224             ! salt flux only       : add concentration dilution term in salt flux  and no  F-M term in volume flux 
    225             ! salt and mass fluxes : non concentration dilution term in salt flux  and add F-M term in volume flux 
    226             sfx (ji,jj) = zfsalt +                  zswitch  * zcd   ! salt flux (+ C/D if no ice/ocean mass exchange) 
    227             emp (ji,jj) = zemp   + zemp_snw + ( 1.- zswitch) * zfmm  ! mass flux (+ F/M mass flux if ice/ocean mass exchange) 
    228             ! 
     185            qns  (ji,jj) = zqns - fdtcn(ji,jj)                           ! non solar heat flux 
    229186         END DO 
    230187      END DO 
    231       !                                !------------------------------------------! 
    232       !                                !    mass of snow and ice per unit area    ! 
    233       !                                !------------------------------------------! 
    234       IF( nn_ice_embd /= 0 ) THEN      ! embedded sea-ice (mass required) 
    235          snwice_mass_b(:,:) = snwice_mass(:,:)                  ! save mass from the previous ice time step 
    236          !                                                      ! new mass per unit area 
    237          snwice_mass  (:,:) = tms(:,:) * ( rhosn * hsnif(:,:) + rhoic * hicif(:,:)  ) * ( 1.0 - frld(:,:) ) 
    238          !                                                      ! time evolution of snow+ice mass 
    239          snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / rdt_ice 
    240       ENDIF 
    241188 
    242189      CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )       
    … …  
    244191      CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1.e0 - pfrld(:,:)) ) 
    245192 
     193      !------------------------------------------! 
     194      !      mass flux at the ocean surface      ! 
     195      !------------------------------------------! 
     196      DO jj = 1, jpj 
     197         DO ji = 1, jpi 
     198            ! 
     199#if defined key_coupled 
     200            ! freshwater exchanges at the ice-atmosphere / ocean interface (coupled mode) 
     201            zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  
     202               &   + rdmsnif(ji,jj) * r1_rdtice                                   !  freshwaterflux due to snow melting  
     203#else 
     204            !  computing freshwater exchanges at the ice/ocean interface 
     205            zemp = + emp(ji,jj)     *         frld(ji,jj)      &   !  e-p budget over open ocean fraction  
     206               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &   !  liquid precipitation reaches directly the ocean 
     207               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  change in ice cover within the time step 
     208               &   + rdmsnif(ji,jj) * r1_rdtice                    !  freshwater flux due to snow melting  
     209#endif             
     210            ! 
     211            !  computing salt exchanges at the ice/ocean interface 
     212            zfons = ( soce_0(ji,jj) - sice_0(ji,jj) ) * ( rdmicif(ji,jj) * r1_rdtice )  
     213            ! 
     214            !  converting the salt flux from ice to a freshwater flux from ocean 
     215            zfm  = zfons / ( sss_m(ji,jj) + epsi16 ) 
     216            ! 
     217            emps(ji,jj) = zemp + zfm      ! surface ocean concentration/dilution effect (use on SSS evolution) 
     218            emp (ji,jj) = zemp            ! surface ocean volume flux (use on sea-surface height evolution) 
     219            ! 
     220         END DO 
     221      END DO 
     222 
    246223      IF( lk_diaar5 ) THEN       ! AR5 diagnostics 
    247          CALL iom_put( 'isnwmlt_cea'  ,                 rdm_snw(:,:) * r1_rdtice ) 
    248          CALL iom_put( 'fsal_virt_cea',   soce_0(:,:) * rdm_ice(:,:) * r1_rdtice ) 
    249          CALL iom_put( 'fsal_real_cea', - sice_0(:,:) * rdm_ice(:,:) * r1_rdtice ) 
     224         CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * r1_rdtice ) 
     225         CALL iom_put( 'fsal_virt_cea',   soce_0(:,:) * rdmicif(:,:) * r1_rdtice ) 
     226         CALL iom_put( 'fsal_real_cea', - sice_0(:,:) * rdmicif(:,:) * r1_rdtice ) 
    250227      ENDIF 
    251228 
    … …  
    267244      IF(ln_ctl) THEN            ! control print 
    268245         CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ') 
    269          CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=sfx   , clinfo2=' sfx     : ') 
     246         CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps  , clinfo2=' emps    : ') 
    270247         CALL prt_ctl(tab2d_1=utau  , clinfo1=' lim_sbc: utau   : ', mask1=umask,   & 
    271248            &         tab2d_2=vtau  , clinfo2=' vtau    : '        , mask2=vmask ) 
    … …  
    463440         END WHERE 
    464441      ENDIF 
    465       !                                      ! embedded sea ice 
    466       IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass 
    467          snwice_mass  (:,:) = tms(:,:) * ( rhosn * hsnif(:,:) + rhoic * hicif(:,:)  ) * ( 1.0 - frld(:,:) ) 
    468          snwice_mass_b(:,:) = snwice_mass(:,:) 
    469       ELSE 
    470          snwice_mass  (:,:) = 0.e0           ! no mass exchanges 
    471          snwice_mass_b(:,:) = 0.e0           ! no mass exchanges 
    472       ENDIF 
    473       IF( nn_ice_embd == 2 .AND.          &  ! full embedment (case 2) & no restart :  
    474          &   .NOT.ln_rstart ) THEN           ! deplete the initial ssh below sea-ice area 
    475          sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0 
    476          sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0 
    477       ENDIF 
    478442      ! 
    479443   END SUBROUTINE lim_sbc_init_2 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90

    r3625 r6736  
    1313   !!   'key_lim2' :                                  LIM 2.0 sea-ice model 
    1414   !!---------------------------------------------------------------------- 
    15    !!   lim_thd_2       : thermodynamic of sea ice 
    16    !!   lim_thd_init_2  : initialisation of sea-ice thermodynamic 
     15   !!   lim_thd_2      : thermodynamic of sea ice 
     16   !!   lim_thd_init_2 : initialisation of sea-ice thermodynamic 
    1717   !!---------------------------------------------------------------------- 
    18    USE phycst           ! physical constants 
    19    USE dom_oce          ! ocean space and time domain variables 
     18   USE phycst          ! physical constants 
     19   USE dom_oce         ! ocean space and time domain variables 
    2020   USE domvvl 
    2121   USE lbclnk 
    22    USE in_out_manager   ! I/O manager 
     22   USE in_out_manager  ! I/O manager 
    2323   USE lib_mpp 
    24    USE wrk_nemo         ! work arrays 
    25    USE iom              ! IOM library 
    26    USE ice_2            ! LIM sea-ice variables 
    27    USE sbc_oce          !  
    28    USE sbc_ice          !  
    29    USE thd_ice_2        ! LIM thermodynamic sea-ice variables 
    30    USE dom_ice_2        ! LIM sea-ice domain 
     24   USE wrk_nemo        ! work arrays 
     25   USE iom             ! IOM library 
     26   USE ice_2           ! LIM sea-ice variables 
     27   USE sbc_oce         !  
     28   USE sbc_ice         !  
     29   USE thd_ice_2       ! LIM thermodynamic sea-ice variables 
     30   USE dom_ice_2       ! LIM sea-ice domain 
    3131   USE limthd_zdf_2 
    3232   USE limthd_lac_2 
    3333   USE limtab_2 
    34    USE prtctl           ! Print control 
    35    USE cpl_oasis3, ONLY :   lk_cpl 
    36    USE diaar5    , ONLY :   lk_diaar5 
    37    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    38     
     34   USE prtctl          ! Print control 
     35   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
     36   USE cpl_oasis3, ONLY : lk_cpl 
     37   USE diaar5, ONLY :   lk_diaar5 
     38       
    3939   IMPLICIT NONE 
    4040   PRIVATE 
    … …  
    5656   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    5757   !!---------------------------------------------------------------------- 
     58 
    5859CONTAINS 
    5960 
    … …  
    8990      REAL(wp) ::   za , zh, zthsnice    ! 
    9091      REAL(wp) ::   zfric_u              ! friction velocity  
     92      REAL(wp) ::   zfnsol               ! total non solar heat 
     93      REAL(wp) ::   zfontn               ! heat flux from snow thickness 
    9194      REAL(wp) ::   zfntlat, zpareff     ! test. the val. of lead heat budget 
    9295 
    … …  
    127130      zdvolif(:,:) = 0.e0   ! total variation of ice volume 
    128131      zdvonif(:,:) = 0.e0   ! transformation of snow to sea-ice volume 
     132!      zdvonif(:,:) = 0.e0   ! lateral variation of ice volume 
    129133      zlicegr(:,:) = 0.e0   ! lateral variation of ice volume 
    130134      zdvomif(:,:) = 0.e0   ! variation of ice volume at bottom due to melting only 
    … …  
    134138      ffltbif(:,:) = 0.e0   ! linked with fstric 
    135139      qfvbq  (:,:) = 0.e0   ! linked with fstric 
    136       rdm_snw(:,:) = 0.e0   ! variation of snow mass over 1 time step 
    137       rdq_snw(:,:) = 0.e0   ! heat content associated with rdm_snw 
    138       rdm_ice(:,:) = 0.e0   ! variation of ice mass over 1 time step 
    139       rdq_ice(:,:) = 0.e0   ! heat content associated with rdm_ice 
     140      rdmsnif(:,:) = 0.e0   ! variation of snow mass per unit area 
     141      rdmicif(:,:) = 0.e0   ! variation of ice mass per unit area 
    140142      zmsk (:,:,:) = 0.e0 
    141143 
    … …  
    198200      !-------------------------------------------------------------------------- 
    199201 
    200       !CDIR NOVERRCHK 
    201       DO jj = 1, jpj 
    202          !CDIR NOVERRCHK 
     202      sst_m(:,:) = sst_m(:,:) + rt0 
     203 
     204!CDIR NOVERRCHK 
     205      DO jj = 1, jpj 
     206!CDIR NOVERRCHK 
    203207         DO ji = 1, jpi 
    204208            zthsnice       = hsnif(ji,jj) + hicif(ji,jj) 
    … …  
    214218            !  temperature and turbulent mixing (McPhee, 1992) 
    215219            zfric_u        = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin )  ! friction velocity 
    216             fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_m(ji,jj) + rt0 - tfu(ji,jj) )  
     220            fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_m(ji,jj) - tfu(ji,jj) )  
    217221            qdtcn(ji,jj)  = zindb * fdtcn(ji,jj) * frld(ji,jj) * rdt_ice 
    218222                         
    219223            !  partial computation of the lead energy budget (qldif) 
    220224#if defined key_coupled  
    221             qldif(ji,jj)   = tms(ji,jj) * rdt_ice                                                  & 
     225            qldif(ji,jj)   = tms(ji,jj) * rdt_ice                                             & 
    222226               &    * (   ( qsr_tot(ji,jj) - qsr_ice(ji,jj,1) * zfricp ) * ( 1.0 - thcm(ji,jj) )   & 
    223227               &        + ( qns_tot(ji,jj) - qns_ice(ji,jj,1) * zfricp )                           & 
    224228               &        + frld(ji,jj) * ( fdtcn(ji,jj) + ( 1.0 - zindb ) * fsbbq(ji,jj) )   ) 
    225229#else 
    226             qldif(ji,jj)   = tms(ji,jj) * rdt_ice * frld(ji,jj)                    & 
    227                &                        * (  qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )    & 
    228                &                           + qns(ji,jj)  +  fdtcn(ji,jj)           & 
    229                &                           + ( 1.0 - zindb ) * fsbbq(ji,jj)      ) 
     230            zfontn         = ( sprecip(ji,jj) / rhosn ) * xlsn  !   energy for melting solid precipitation 
     231            zfnsol         = qns(ji,jj)                         !  total non solar flux over the ocean 
     232            qldif(ji,jj)   = tms(ji,jj) * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )   & 
     233               &                               + zfnsol + fdtcn(ji,jj) - zfontn     & 
     234               &                               + ( 1.0 - zindb ) * fsbbq(ji,jj) )   & 
     235               &                        * frld(ji,jj) * rdt_ice     
     236!!$            qldif(ji,jj)   = tms(ji,jj) * rdt_ice * frld(ji,jj)  
     237!!$               &           * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )      & 
     238!!$               &             + qns(ji,jj)  + fdtcn(ji,jj) - zfontn     & 
     239!!$               &             + ( 1.0 - zindb ) * fsbbq(ji,jj)      )   & 
    230240#endif 
    231241            !  parlat : percentage of energy used for lateral ablation (0.0)  
    … …  
    237247             
    238248            !  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 ) 
     249            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1)   & 
     250                &          * ( tfu(ji,jj) - sst_m(ji,jj) ) * ( 1 - zinda ) 
    240251             
    241252            !  calculate oceanic heat flux. 
    … …  
    247258      END DO 
    248259       
     260      sst_m(:,:) = sst_m(:,:) - rt0 
     261                
    249262      !         Select icy points and fulfill arrays for the vectorial grid. 
    250263      !---------------------------------------------------------------------- 
    … …  
    300313         CALL tab_2d_1d_2( nbpb, qldif_1d   (1:nbpb)     , qldif      , jpi, jpj, npb(1:nbpb) ) 
    301314         CALL tab_2d_1d_2( nbpb, qstbif_1d  (1:nbpb)     , qstoif     , jpi, jpj, npb(1:nbpb) ) 
    302          CALL tab_2d_1d_2( nbpb, rdm_ice_1d (1:nbpb)     , rdm_ice    , jpi, jpj, npb(1:nbpb) ) 
    303          CALL tab_2d_1d_2( nbpb, rdq_ice_1d (1:nbpb)     , rdq_ice    , jpi, jpj, npb(1:nbpb) ) 
     315         CALL tab_2d_1d_2( nbpb, rdmicif_1d (1:nbpb)     , rdmicif    , jpi, jpj, npb(1:nbpb) ) 
    304316         CALL tab_2d_1d_2( nbpb, dmgwi_1d   (1:nbpb)     , dmgwi      , jpi, jpj, npb(1:nbpb) ) 
    305          CALL tab_2d_1d_2( nbpb, rdm_snw_1d (1:nbpb)     , rdm_snw    , jpi, jpj, npb(1:nbpb) ) 
    306          CALL tab_2d_1d_2( nbpb, rdq_snw_1d (1:nbpb)     , rdq_snw    , jpi, jpj, npb(1:nbpb) ) 
    307317         CALL tab_2d_1d_2( nbpb, qlbbq_1d   (1:nbpb)     , zqlbsbq    , jpi, jpj, npb(1:nbpb) ) 
    308318         ! 
    … …  
    323333         CALL tab_1d_2d_2( nbpb, qfvbq      , npb, qfvbq_1d  (1:nbpb)     , jpi, jpj ) 
    324334         CALL tab_1d_2d_2( nbpb, qstoif     , npb, qstbif_1d (1:nbpb)     , jpi, jpj ) 
    325          CALL tab_1d_2d_2( nbpb, rdm_ice    , npb, rdm_ice_1d(1:nbpb)     , jpi, jpj ) 
    326          CALL tab_1d_2d_2( nbpb, rdq_ice    , npb, rdq_ice_1d(1:nbpb)     , jpi, jpj ) 
     335         CALL tab_1d_2d_2( nbpb, rdmicif    , npb, rdmicif_1d(1:nbpb)     , jpi, jpj ) 
    327336         CALL tab_1d_2d_2( nbpb, dmgwi      , npb, dmgwi_1d  (1:nbpb)     , jpi, jpj ) 
    328          CALL tab_1d_2d_2( nbpb, rdm_snw    , npb, rdm_snw_1d(1:nbpb)     , jpi, jpj ) 
    329          CALL tab_1d_2d_2( nbpb, rdq_snw    , npb, rdq_snw_1d(1:nbpb)     , jpi, jpj ) 
     337         CALL tab_1d_2d_2( nbpb, rdmsnif    , npb, rdmsnif_1d(1:nbpb)     , jpi, jpj ) 
    330338         CALL tab_1d_2d_2( nbpb, zdvosif    , npb, dvsbq_1d  (1:nbpb)     , jpi, jpj ) 
    331339         CALL tab_1d_2d_2( nbpb, zdvobif    , npb, dvbbq_1d  (1:nbpb)     , jpi, jpj ) 
    … …  
    386394      IF( nbpac > 0 ) THEN 
    387395         ! 
    388          zlicegr(:,:) = rdm_ice(:,:)      ! to output the lateral sea-ice growth  
     396         zlicegr(:,:) = rdmicif(:,:)      ! to output the lateral sea-ice growth  
    389397         !...Put the variable in a 1-D array for lateral accretion 
    390398         CALL tab_2d_1d_2( nbpac, frld_1d   (1:nbpac)     , frld       , jpi, jpj, npac(1:nbpac) ) 
    … …  
    397405         CALL tab_2d_1d_2( nbpac, qcmif_1d  (1:nbpac)     , qcmif      , jpi, jpj, npac(1:nbpac) ) 
    398406         CALL tab_2d_1d_2( nbpac, qstbif_1d (1:nbpac)     , qstoif     , jpi, jpj, npac(1:nbpac) ) 
    399          CALL tab_2d_1d_2( nbpac, rdm_ice_1d(1:nbpac)     , rdm_ice    , jpi, jpj, npac(1:nbpac) ) 
    400          CALL tab_2d_1d_2( nbpac, rdq_ice_1d(1:nbpac)     , rdq_ice    , jpi, jpj, npac(1:nbpac) ) 
     407         CALL tab_2d_1d_2( nbpac, rdmicif_1d(1:nbpac)     , rdmicif    , jpi, jpj, npac(1:nbpac) ) 
    401408         CALL tab_2d_1d_2( nbpac, dvlbq_1d  (1:nbpac)     , zdvolif    , jpi, jpj, npac(1:nbpac) ) 
    402409         CALL tab_2d_1d_2( nbpac, tfu_1d    (1:nbpac)     , tfu        , jpi, jpj, npac(1:nbpac) ) 
    … …  
    412419         CALL tab_1d_2d_2( nbpac, tbif(:,:,3), npac(1:nbpac), tbif_1d   (1:nbpac , 3 ), jpi, jpj ) 
    413420         CALL tab_1d_2d_2( nbpac, qstoif     , npac(1:nbpac), qstbif_1d (1:nbpac)     , jpi, jpj ) 
    414          CALL tab_1d_2d_2( nbpac, rdm_ice    , npac(1:nbpac), rdm_ice_1d(1:nbpac)     , jpi, jpj ) 
    415          CALL tab_1d_2d_2( nbpac, rdq_ice    , npac(1:nbpac), rdq_ice_1d(1:nbpac)     , jpi, jpj ) 
     421         CALL tab_1d_2d_2( nbpac, rdmicif    , npac(1:nbpac), rdmicif_1d(1:nbpac)     , jpi, jpj ) 
    416422         CALL tab_1d_2d_2( nbpac, zdvolif    , npac(1:nbpac), dvlbq_1d  (1:nbpac)     , jpi, jpj ) 
    417423         ! 
    … …  
    444450      CALL iom_put( 'iceprod_cea' , hicifp (:,:) * zztmp     )   ! Ice produced               [m/s] 
    445451      IF( lk_diaar5 ) THEN 
    446          CALL iom_put( 'snowmel_cea' , rdm_snw(:,:) * zztmp     )   ! Snow melt                  [kg/m2/s] 
     452         CALL iom_put( 'snowmel_cea' , rdmsnif(:,:) * zztmp     )   ! Snow melt                  [kg/m2/s] 
    447453         zztmp = rhoic / rdt_ice 
    448454         CALL iom_put( 'sntoice_cea' , zdvonif(:,:) * zztmp     )   ! Snow to Ice transformation [kg/m2/s] 
    449455         CALL iom_put( 'ticemel_cea' , zdvosif(:,:) * zztmp     )   ! Melt at Sea Ice top        [kg/m2/s] 
    450456         CALL iom_put( 'bicemel_cea' , zdvomif(:,:) * zztmp     )   ! Melt at Sea Ice bottom     [kg/m2/s] 
    451          zlicegr(:,:) = MAX( 0.e0, rdm_ice(:,:)-zlicegr(:,:) ) 
    452          CALL iom_put( 'licepro_cea' , zlicegr(:,:) * zztmp     )   ! Lateral sea ice growth     [kg/m2/s] 
     457         zlicegr(:,:) = MAX( 0.e0, rdmicif(:,:)-zlicegr(:,:) ) 
     458         CALL iom_put( 'licepro_cea' , zlicegr(:,:) * zztmp     )   ! Latereal sea ice growth    [kg/m2/s] 
    453459      ENDIF 
    454460      ! 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/limthd_lac_2.F90

    r3625 r6736  
    77 
    88   !!---------------------------------------------------------------------- 
    9    !!   lim_lat_acr_2 : lateral accretion of ice 
    10    !!---------------------------------------------------------------------- 
    11    USE par_oce        ! ocean parameters 
     9   !!   lim_lat_acr_2   : lateral accretion of ice 
     10   !!---------------------------------------------------------------------- 
     11   USE par_oce          ! ocean parameters 
    1212   USE phycst 
    1313   USE thd_ice_2 
    1414   USE ice_2 
    1515   USE limistate_2  
    16    USE lib_mpp        ! MPP library 
    17    USE wrk_nemo       ! work arrays 
    18    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     16   USE lib_mpp          ! MPP library 
     17   USE wrk_nemo         ! work arrays 
     18   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    1919 
    2020   IMPLICIT NONE 
    … …  
    146146         frld_1d   (ji) = MAX( zfrlnew , zfrlmin(ji) ) 
    147147         !--computation of the remaining part of ice thickness which has been already used 
    148          zdhicbot(ji) =  ( frld_1d(ji) - zfrlnew ) * zhice0(ji) / ( 1.0 - zfrlmin(ji) )   &  
    149             &         -  (  ( 1.0 - zfrrate ) / ( 1.0 - frld_1d(ji) ) )  * ( zqbgow(ji) / xlic )  
     148         zdhicbot(ji) =  ( frld_1d(ji) - zfrlnew ) * zhice0(ji) / ( 1.0 - zfrlmin(ji) ) &  
     149                      -  (  ( 1.0 - zfrrate ) / ( 1.0 - frld_1d(ji) ) )  * ( zqbgow(ji) / xlic )  
    150150      END DO 
    151151  
    … …  
    197197            &          ) / zah 
    198198          
    199          tbif_1d(ji,3) =     ( iiceform * ( zhnews2 - zdh3 )                                           * zta1  & 
     199         tbif_1d(ji,3) =     (  iiceform * ( zhnews2 - zdh3 )                                          * zta1  & 
    200200            &              + ( iiceform * zdh3 + ( 1 - iiceform ) * zdh1 )                             * zta2  & 
    201201            &              + ( iiceform * ( zhnews2 - zdh5 ) + ( 1 - iiceform ) * ( zhnews2 - zdh1 ) ) * zta3  &  
    … …  
    218218      DO ji = kideb , kiut 
    219219         dvlbq_1d  (ji) = ( 1. - frld_1d(ji) ) * h_ice_1d(ji) - ( 1. - zfrl_old(ji) ) * zhice_old(ji) 
    220          rdm_ice_1d(ji) = rdm_ice_1d(ji) + rhoic * dvlbq_1d(ji) 
    221          rdq_ice_1d(ji) = rdq_ice_1d(ji) + rcpic * dvlbq_1d(ji) * ( tfu_1d(ji) - rt0 )      ! heat content relative to rt0 
     220         rdmicif_1d(ji) = rdmicif_1d(ji) + rhoic * dvlbq_1d(ji) 
    222221      END DO 
    223222       

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90

    r3625 r6736  
    1818   USE ice_2 
    1919   USE limistate_2 
    20    USE cpl_oasis3, ONLY : lk_cpl 
    2120   USE in_out_manager 
    2221   USE lib_mpp          ! MPP library 
    2322   USE wrk_nemo         ! work arrays 
    24    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    25      
     23   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
     24   USE cpl_oasis3, ONLY : lk_cpl 
     25       
    2626   IMPLICIT NONE 
    2727   PRIVATE 
    … …  
    8787      REAL(wp), POINTER, DIMENSION(:) ::   zrcpdt         ! h_su*rho_su*cp_su/dt(h_su being the thick. of surf. layer) 
    8888      REAL(wp), POINTER, DIMENSION(:) ::   zts_old        ! previous surface temperature 
    89       REAL(wp), POINTER, DIMENSION(:) ::   zidsn , z1midsn , zidsnic ! temporary variables 
     89      REAL(wp), POINTER, DIMENSION(:) ::   zidsn , z1midsn , zidsnic ! tempory variables 
    9090      REAL(wp), POINTER, DIMENSION(:) ::   zfnet          ! net heat flux at the top surface( incl. conductive heat flux) 
    9191      REAL(wp), POINTER, DIMENSION(:) ::   zsprecip       ! snow accumulation 
    … …  
    9999      REAL(wp), POINTER, DIMENSION(:) ::   zep            ! internal temperature of the 2nd layer of the snow/ice system 
    100100      REAL(wp), DIMENSION(3) :: &  
    101             zplediag  &    ! principle diagonal, subdiag. and supdiag. of the  
     101          zplediag  &    ! principle diagonal, subdiag. and supdiag. of the  
    102102          , zsubdiag  &    ! tri-diagonal matrix coming from the computation 
    103103          , zsupdiag  &    ! of the temperatures inside the snow-ice system 
    104104          , zsmbr          ! second member 
    105        REAL(wp) ::    & 
    106             zhsu      &    ! thickness of surface layer 
    107           , zhe       &    ! effective thickness for compu. of equ. thermal conductivity 
    108           , zheshth   &    ! = zhe / thth 
    109           , zghe      &    ! correction factor of the thermal conductivity 
    110           , zumsb     &    ! parameter for numerical method to solve heat-diffusion eq. 
    111           , zkhsn     &    ! conductivity at the snow layer 
    112           , zkhic     &    ! conductivity at the ice layers 
    113           , zkint     &    ! equivalent conductivity at the snow-ice interface 
    114           , zkhsnint  &    ! = zkint*dt / (hsn*rhosn*cpsn)   
    115           , zkhicint  &    ! = 2*zkint*dt / (hic*rhoic*cpic) 
    116           , zpiv1, zpiv2 & ! temporary scalars used to solve the tri-diagonal system 
    117           , zb2, zd2  &    ! temporary scalars used to solve the tri-diagonal system 
    118           , zb3, zd3  &    ! temporary scalars used to solve the tri-diagonal system 
     105       REAL(wp) :: &  
     106          zhsu     &     ! thickness of surface layer 
     107          , zhe      &     ! effective thickness for compu. of equ. thermal conductivity 
     108          , zheshth  &     ! = zhe / thth 
     109          , zghe     &     ! correction factor of the thermal conductivity 
     110          , zumsb    &     ! parameter for numerical method to solve heat-diffusion eq. 
     111          , zkhsn    &     ! conductivity at the snow layer 
     112          , zkhic    &     ! conductivity at the ice layers 
     113          , zkint    &     ! equivalent conductivity at the snow-ice interface 
     114          , zkhsnint &     ! = zkint*dt / (hsn*rhosn*cpsn)   
     115          , zkhicint &     ! = 2*zkint*dt / (hic*rhoic*cpic) 
     116          , zpiv1 , zpiv2  &       ! tempory scalars used to solve the tri-diagonal system 
     117          , zb2 , zd2 , zb3 , zd3 & 
    119118          , ztint          ! equivalent temperature at the snow-ice interface 
    120        REAL(wp) ::    &  
    121             zexp      &    ! exponential function of the ice thickness 
    122           , zfsab     &    ! part of solar radiation stored in brine pockets 
    123           , zfts      &    ! value of energy balance function when the temp. equal surf. temp. 
    124           , zdfts     &    ! value of derivative of ztfs when the temp. equal surf. temp. 
    125           , zdts      &    ! surface temperature increment 
    126           , zqsnw_mlt &    ! energy needed to melt snow 
    127           , zdhsmlt   &    ! change in snow thickness due to melt 
    128           , zhsn      &    ! snow thickness (previous+accumulation-melt) 
    129           , zqsn_mlt_rem & ! remaining heat coming from snow melting 
    130           , zqice_top_mlt &! energy used to melt ice at top surface 
    131           , zdhssub   &    ! change in snow thick. due to sublimation or evaporation 
    132           , zdhisub   &    ! change in ice thick. due to sublimation or evaporation     
    133           , zdhsn     &    ! snow ice thickness increment 
    134           , zdtsn     &    ! snow internal temp. increment 
    135           , zdtic     &    ! ice internal temp. increment 
     119       REAL(wp) :: &  
     120          zexp      &     ! exponential function of the ice thickness 
     121          , zfsab     &     ! part of solar radiation stored in brine pockets 
     122          , zfts      &     ! value of energy balance function when the temp. equal surf. temp. 
     123          , zdfts     &     ! value of derivative of ztfs when the temp. equal surf. temp. 
     124          , zdts      &     ! surface temperature increment 
     125          , zqsnw_mlt &     ! energy needed to melt snow 
     126          , zdhsmlt   &     ! change in snow thickness due to melt 
     127          , zhsn      &     ! snow thickness (previous+accumulation-melt) 
     128          , zqsn_mlt_rem &  ! remaining heat coming from snow melting 
     129          , zqice_top_mlt & ! energy used to melt ice at top surface 
     130          , zdhssub      &  ! change in snow thick. due to sublimation or evaporation 
     131          , zdhisub      &  ! change in ice thick. due to sublimation or evaporation     
     132          , zdhsn        &  ! snow ice thickness increment 
     133          , zdtsn        &  ! snow internal temp. increment 
     134          , zdtic        &  ! ice internal temp. increment 
    136135          , zqnes          ! conductive energy due to ice melting in the first ice layer 
    137        REAL(wp) ::    &  
    138             ztbot     &    ! temperature at the bottom surface 
    139           , zfcbot    &    ! conductive heat flux at bottom surface 
    140           , zqice_bot &    ! energy used for bottom melting/growing 
    141           , zqice_bot_mlt &! energy used for bottom melting 
    142           , zqstbif_bot  & ! part of energy stored in brine pockets used for bottom melting 
    143           , zqstbif_old  & ! temporary var. for zqstbif_bot 
    144           , zdhicmlt  &    ! change in ice thickness due to bottom melting 
    145           , zdhicm    &    ! change in ice thickness var.  
    146           , zdhsnm    &    ! change in snow thickness var.  
    147           , zhsnfi    &    ! snow thickness var.  
    148           , zc1, zpc1 &    ! temporary variables 
    149           , zc2, zpc2 &    ! temporary variables 
    150           , zp1, zp2  &    ! temporary variables 
    151           , ztb2, ztb3     ! temporary variables 
    152        REAL(wp) ::    &  
    153             zdrmh     &    ! change in snow/ice thick. after snow-ice formation 
    154           , zhicnew   &    ! new ice thickness 
    155           , zhsnnew   &    ! new snow thickness 
    156           , zquot     & 
    157           , ztneq     &    ! temporary temp. variables 
    158           , zqice     & 
    159           , zqicetot  &    ! total heat inside the snow/ice system 
    160           , zdfrl     &    ! change in ice concentration 
    161           , zdvsnvol  &    ! change in snow volume 
    162           , zdrfrl1, zdrfrl2, zihsn, zidhb, zihic &  ! temporary scalars 
    163           , zihe, zihq, ziexp, ziqf, zihnf        &  ! temporary scalars 
    164           , zibmlt, ziqr, zihgnew, zind, ztmp        ! temporary scalars 
     136       REAL(wp) :: &  
     137          ztbot     &      ! temperature at the bottom surface 
     138          , zfcbot    &      ! conductive heat flux at bottom surface 
     139          , zqice_bot &      ! energy used for bottom melting/growing 
     140          , zqice_bot_mlt &  ! energy used for bottom melting 
     141          , zqstbif_bot  &  ! part of energy stored in brine pockets used for bottom melting 
     142          , zqstbif_old  &  ! tempory var. for zqstbif_bot 
     143          , zdhicmlt      &  ! change in ice thickness due to bottom melting 
     144          , zdhicm        &  ! change in ice thickness var.  
     145          , zdhsnm        &  ! change in snow thickness var.  
     146          , zhsnfi        &  ! snow thickness var.  
     147          , zc1, zpc1, zc2, zpc2, zp1, zp2 & ! tempory variables 
     148          , ztb2, ztb3 
     149       REAL(wp) :: &  
     150          zdrmh         &   ! change in snow/ice thick. after snow-ice formation 
     151          , zhicnew       &   ! new ice thickness 
     152          , zhsnnew       &   ! new snow thickness 
     153          , zquot , ztneq &   ! tempory temp. variables 
     154          , zqice, zqicetot & ! total heat inside the snow/ice system 
     155          , zdfrl         &   ! change in ice concentration 
     156          , zdvsnvol      &   ! change in snow volume 
     157          , zdrfrl1, zdrfrl2 &  ! tempory scalars 
     158          , zihsn, zidhb, zihic, zihe, zihq, ziexp, ziqf, zihnf, zibmlt, ziqr, zihgnew, zind 
    165159       !!---------------------------------------------------------------------- 
    166160       CALL wrk_alloc( jpij, ztsmlt, ztbif  , zksn    , zkic    , zksndh , zfcsu  , zfcsudt , zi0      , z1mi0   , zqmax    ) 
    … …  
    176170        
    177171       DO ji = kideb , kiut 
    178           ! do nothing if the snow (ice) thickness falls below its minimum thickness 
    179172          zihsn = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) ) 
    180173          zihic = MAX( zzero , SIGN( zone , hicdif - h_ice_1d(ji) ) ) 
    181           !--energy required to bring snow to its melting point (rt0_snow) 
    182           zqcmlts(ji) = ( MAX ( zzero , rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn ) 
    183           !--energy required to bring ice to its melting point (rt0_ice) 
    184           zqcmltb(ji) = ( MAX( zzero , rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) )  & 
    185              &          + MAX( zzero , rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) )  & 
    186              &          ) * ( 1.0 - zihic  ) 
    187           !--limitation of snow/ice system internal temperature 
     174          !--computation of energy due to surface melting 
     175          zqcmlts(ji) = ( MAX ( zzero ,  & 
     176             &                   rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn ) 
     177          !--computation of energy due to bottom melting 
     178          zqcmltb(ji) = ( MAX( zzero , & 
     179             &                  rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) & 
     180             &           + MAX( zzero , & 
     181             &                  rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) & 
     182             &           ) * ( 1.0 - zihic  ) 
     183          !--limitation of  snow/ice system internal temperature 
    188184          tbif_1d(ji,1)   = MIN( rt0_snow, tbif_1d(ji,1) ) 
    189185          tbif_1d(ji,2)   = MIN( rt0_ice , tbif_1d(ji,2) ) 
    … …  
    485481          dvsbq_1d(ji) =  ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnw_old(ji) - zsprecip(ji) ) 
    486482          dvsbq_1d(ji) =  MIN( zzero , dvsbq_1d(ji) ) 
    487           ztmp = rhosn * dvsbq_1d(ji) 
    488           rdm_snw_1d(ji) =  ztmp 
    489           !--heat content of the water provided to the ocean (referenced to rt0) 
    490           rdq_snw_1d(ji) =  cpic * ztmp * ( rt0_snow - rt0 ) 
     483          rdmsnif_1d(ji) =  rhosn * dvsbq_1d(ji) 
    491484          !-- If the snow is completely melted the remaining heat is used to melt ice 
    492485          zqsn_mlt_rem  = MAX( zzero , -zhsn ) * xlsn 
    … …  
    631624          !---updating new ice thickness and computing the newly formed ice mass 
    632625          zhicnew   =  zihgnew * zhicnew 
    633           ztmp    =  ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) * rhoic 
    634           rdm_ice_1d(ji) =  rdm_ice_1d(ji) + ztmp 
    635           !---heat content of the water provided to the ocean (referenced to rt0) 
    636           !   use of rt0_ice is OK for melting ice; in the case of freezing, tfu_1d should be used.  
    637           !   This is done in 9.5 section (see below) 
    638           rdq_ice_1d(ji) =  cpic * ztmp * ( rt0_ice - rt0 ) 
     626          rdmicif_1d(ji) =  rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) ) * rhoic 
    639627          !---updating new snow thickness and computing the newly formed snow mass 
    640628          zhsnfi   = zhsn + zdhsnm 
    641629          h_snow_1d(ji) = MAX( zzero , zhsnfi ) 
    642           ztmp = ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsn ) * rhosn 
    643           rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp 
    644           !---updating the heat content of the water provided to the ocean (referenced to rt0) 
    645           rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( rt0_snow - rt0 ) 
     630          rdmsnif_1d(ji) =  rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) ) * ( h_snow_1d(ji) - zhsn ) * rhosn 
    646631          !--remaining energy in case of total ablation 
    647632          zqocea(ji) = - ( zihsn * xlic * zdhicm + xlsn * ( zhsnfi - h_snow_1d(ji) ) ) * ( 1.0 - frld_1d(ji) ) 
    … …  
    675660          tbif_1d(ji,3) =  zihgnew * ztb3 + ( 1.0 - zihgnew ) * tfu_1d(ji) 
    676661          h_ice_1d(ji)  =  zhicnew 
    677           ! update the ice heat content given to the ocean in freezing case  
    678           ! (part due to difference between rt0_ice and tfu_1d) 
    679           ztmp = ( 1. - zidhb ) * rhoic * dvbbq_1d(ji) 
    680           rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * ( tfu_1d(ji) - rt0_ice ) 
    681662       END DO 
    682663 
    … …  
    720701          dmgwi_1d(ji) = dmgwi_1d(ji) + ( 1.0 -frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnnew ) * rhosn 
    721702          !---  volume change of ice and snow (used for ocean-ice freshwater flux computation) 
    722           ztmp = ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d (ji) ) * rhoic 
    723           rdm_ice_1d(ji) = rdm_ice_1d(ji) + ztmp 
    724           rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * ( tfu_1d(ji) - rt0 ) 
    725           !!gm BUG ??   snow ==>  only needed for nn_ice_embd == 0  (standard levitating sea-ice) 
    726           ztmp = ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn 
    727           rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp 
    728           rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( rt0_snow - rt0 ) 
     703          rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhicnew - h_ice_1d (ji) ) * rhoic 
     704          rdmsnif_1d(ji) = rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn 
    729705 
    730706          !---  Actualize new snow and ice thickness. 
    … …  
    773749          !--variation of ice volume and ice mass  
    774750          dvlbq_1d(ji)   = zihic * ( zfrl_old(ji) - frld_1d(ji) ) * h_ice_1d(ji) 
    775           ztmp = dvlbq_1d(ji) * rhoic 
    776           rdm_ice_1d(ji) = rdm_ice_1d(ji) + ztmp 
    777 !!gm 
    778 !!gm   This should be split in two parts: 
    779 !!gm         1-  heat required to bring sea-ice to tfu  : this part should be added to the heat flux taken from the ocean 
    780 !!gm                 cpic * ztmp * 0.5 * ( tbif_1d(ji,2) + tbif_1d(ji,3) - 2.* rt0_ice ) 
    781 !!gm         2-  heat content of lateral ablation referenced to rt0 : this part only put in rdq_ice_1d 
    782 !!gm                 cpic * ztmp * ( rt0_ice - rt0 ) 
    783 !!gm   Currently we put all the heat in rdq_ice_1d 
    784           rdq_ice_1d(ji) = rdq_ice_1d(ji) + cpic * ztmp * 0.5 * ( tbif_1d(ji,2) + tbif_1d(ji,3) - 2.* rt0 ) 
    785           ! 
     751          rdmicif_1d(ji) = rdmicif_1d(ji) + dvlbq_1d(ji) * rhoic 
    786752          !--variation of snow volume and snow mass  
    787           zdvsnvol = zihsn * ( zfrl_old(ji) - frld_1d(ji) ) * h_snow_1d(ji) 
    788           ztmp     = zdvsnvol * rhosn 
    789           rdm_snw_1d(ji) = rdm_snw_1d(ji) + ztmp 
    790 !!gm 
    791 !!gm   This should be split in two parts: 
    792 !!gm         1-  heat required to bring snow to tfu  : this part should be added to the heat flux taken from the ocean 
    793 !!gm                 cpic * ztmp * ( tbif_1d(ji,1) - rt0_snow ) 
    794 !!gm         2-  heat content of lateral ablation referenced to rt0 : this part only put in rdq_snw_1d 
    795 !!gm                 cpic * ztmp * ( rt0_snow - rt0 ) 
    796 !!gm   Currently we put all the heat in rdq_snw_1d 
    797           rdq_snw_1d(ji) = rdq_snw_1d(ji) + cpic * ztmp * ( tbif_1d(ji,1) - rt0 ) 
    798  
     753          zdvsnvol    = zihsn * ( zfrl_old(ji) - frld_1d(ji) ) * h_snow_1d(ji) 
     754          rdmsnif_1d(ji) = rdmsnif_1d(ji) + zdvsnvol * rhosn 
    799755          h_snow_1d(ji)  = ziqf * h_snow_1d(ji) 
    800756 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/limtrp_2.F90

    r3764 r6736  
    2828   USE lib_mpp         ! MPP library 
    2929   USE wrk_nemo        ! work arrays 
    30 # if defined key_agrif 
    31    USE agrif_lim2_interp ! nesting 
    32 # endif 
    3330   USE lib_fortran     ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    3431 
    … …  
    8481 
    8582      IF( kt == nit000  )   CALL lim_trp_init_2      ! Initialization (first time-step only) 
    86  
    87 # if defined key_agrif 
    88       CALL agrif_trp_lim2_load      ! First interpolation 
    89 # endif 
    9083 
    9184      zsm(:,:) = area(:,:) 
    … …  
    277270      ENDIF 
    278271      ! 
    279 # if defined key_agrif 
    280       CALL agrif_trp_lim2      ! Fill boundaries of the fine grid 
    281 # endif 
    282       !  
    283272      CALL wrk_dealloc( jpi, jpj, zui_u , zvi_v , zsm, zs0ice, zs0sn , zs0a, zs0c0 , zs0c1 , zs0c2 , zs0st ) 
    284273      ! 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/limwri_2.F90

    r3764 r6736  
    1313   !!---------------------------------------------------------------------- 
    1414   !!---------------------------------------------------------------------- 
    15    !!   lim_wri_2       : write of the diagnostics variables in ouput file  
    16    !!   lim_wri_init_2  : initialization and namelist read 
     15   !!   lim_wri_2      : write of the diagnostics variables in ouput file  
     16   !!   lim_wri_init_2 : initialization and namelist read 
    1717   !!   lim_wri_state_2 : write for initial state or/and abandon: 
    1818   !!                     > output.init.nc (if ninist = 1 in namelist) 
    … …  
    2626   USE ice_2 
    2727 
    28    USE dianam           ! build name of file (routine) 
     28   USE dianam          ! build name of file (routine) 
    2929   USE lbclnk 
    3030   USE in_out_manager 
    31    USE lib_mpp          ! MPP library 
    32    USE wrk_nemo         ! work arrays 
     31   USE lib_mpp         ! MPP library 
     32   USE wrk_nemo        ! work arrays 
    3333   USE iom 
    3434   USE ioipsl 
    35    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3635 
    3736   IMPLICIT NONE 
    … …  
    185184            zcmo(ji,jj,13) = qns(ji,jj) 
    186185            ! See thersf for the coefficient 
    187             zcmo(ji,jj,14) = - sfx (ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce    !!gm ??? 
     186            zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce    !!gm ??? 
    188187            zcmo(ji,jj,15) = utau_ice(ji,jj) 
    189188            zcmo(ji,jj,16) = vtau_ice(ji,jj) 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/limwri_dimg_2.h90

    r3764 r6736  
    125125          zcmo(ji,jj,13) = qns(ji,jj) 
    126126          ! See thersf for the coefficient 
    127           zcmo(ji,jj,14) = - sfx (ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 
     127          zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 
    128128          zcmo(ji,jj,15) = utau_ice(ji,jj) 
    129129          zcmo(ji,jj,16) = vtau_ice(ji,jj) 
    … …  
    173173                rcmoy(ji,jj,13) = qns(ji,jj) 
    174174                ! See thersf for the coefficient 
    175                 rcmoy(ji,jj,14) = - sfx (ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 
     175                rcmoy(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce 
    176176                rcmoy(ji,jj,15) = utau_ice(ji,jj) 
    177177                rcmoy(ji,jj,16) = vtau_ice(ji,jj) 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_2/thd_ice_2.F90

    r3625 r6736  
    6868      qstbif_1d   ,     &  !:    "                  "      qstoif 
    6969      fbif_1d     ,     &  !:    "                  "      fbif 
    70       rdm_ice_1d  ,     &  !:    "                  "      rdm_ice 
    71       rdq_ice_1d  ,     &  !:    "                  "      rdq_ice 
    72       rdm_snw_1d  ,     &  !:    "                  "      rdm_snw 
    73       rdq_snw_1d  ,     &  !:    "                  "      rdq_snw 
     70      rdmicif_1d  ,     &  !:    "                  "      rdmicif 
     71      rdmsnif_1d  ,     &  !:    "                  "      rdmsnif 
    7472      qlbbq_1d    ,     &  !:    "                  "      qlbsbq 
    7573      dmgwi_1d    ,     &  !:    "                  "      dmgwi 
    … …  
    110108         &      qstbif_1d(jpij),  fbif_1d(jpij),  Stat=ierr(2)) 
    111109         ! 
    112       ALLOCATE( rdm_ice_1d(jpij), rdq_ice_1d(jpij)                  , & 
    113          &      rdm_snw_1d(jpij), rdq_snw_1d(jpij), qlbbq_1d(jpij)  , & 
     110      ALLOCATE( rdmicif_1d(jpij), rdmsnif_1d(jpij), qlbbq_1d(jpij),   & 
    114111         &      dmgwi_1d(jpij)  , dvsbq_1d(jpij)  , rdvomif_1d(jpij), & 
    115112         &      dvbbq_1d(jpij)  , dvlbq_1d(jpij)  , dvnbq_1d(jpij)  , & 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r3791 r6736  
    88   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy  
    99   !!            3.3  !  2009-05  (G.Garric) addition of the lim2_evp cas 
    10    !!            3.4  !  2011-01  (A. Porter)  dynamical allocation  
    11    !!            3.5  !  2012-08  (R. Benshila)  AGRIF  
     10   !!            3.4  !  2011-01  (A Porter)  dynamical allocation  
    1211   !!---------------------------------------------------------------------- 
    1312#if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp ) 
    … …  
    1615   !!   'key_lim2' AND NOT 'key_lim2_vp'            EVP LIM-2 sea-ice model 
    1716   !!---------------------------------------------------------------------- 
    18    !!   lim_rhg       : computes ice velocities 
     17   !!   lim_rhg   : computes ice velocities 
    1918   !!---------------------------------------------------------------------- 
    20    USE phycst         ! Physical constant 
    21    USE oce     , ONLY :  snwice_mass, snwice_mass_b 
    22    USE par_oce        ! Ocean parameters 
    23    USE dom_oce        ! Ocean domain 
    24    USE sbc_oce        ! Surface boundary condition: ocean fields 
    25    USE sbc_ice        ! Surface boundary condition: ice fields 
     19   USE phycst           ! Physical constant 
     20   USE par_oce          ! Ocean parameters 
     21   USE dom_oce          ! Ocean domain 
     22   USE sbc_oce          ! Surface boundary condition: ocean fields 
     23   USE sbc_ice          ! Surface boundary condition: ice fields 
     24   USE lbclnk           ! Lateral Boundary Condition / MPP link 
     25   USE lib_mpp          ! MPP library 
     26   USE wrk_nemo         ! work arrays 
     27   USE in_out_manager   ! I/O manager 
     28   USE prtctl           ! Print control 
     29   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
    2630#if defined key_lim3 
    27    USE ice            ! LIM-3: ice variables 
    28    USE dom_ice        ! LIM-3: ice domain 
    29    USE limitd_me      ! LIM-3:  
     31   USE ice              ! LIM-3: ice variables 
     32   USE dom_ice          ! LIM-3: ice domain 
     33   USE limitd_me        ! LIM-3:  
    3034#else 
    31    USE ice_2          ! LIM-2: ice variables 
    32    USE dom_ice_2      ! LIM-2: ice domain 
    33 #endif 
    34    USE lbclnk         ! Lateral Boundary Condition / MPP link 
    35    USE lib_mpp        ! MPP library 
    36    USE wrk_nemo       ! work arrays 
    37    USE in_out_manager ! I/O manager 
    38    USE prtctl         ! Print control 
    39    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    40 #if defined key_agrif && defined key_lim2 
    41    USE agrif_lim2_interp 
     35   USE ice_2            ! LIM2: ice variables 
     36   USE dom_ice_2        ! LIM2: ice domain 
    4237#endif 
    4338 
    … …  
    130125      REAL(wp) ::   zindb         ! ice (1) or not (0)       
    131126      REAL(wp) ::   zdummy        ! dummy argument 
    132       REAL(wp) ::   zintb, zintn  ! dummy argument 
    133127 
    134128      REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh           ! temporary array for ice strength 
    … …  
    152146      REAL(wp), POINTER, DIMENSION(:,:) ::   zs12             ! Non-diagonal stress tensor component zs12 
    153147      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr   ! Local error on velocity 
    154       REAL(wp), POINTER, DIMENSION(:,:) ::   zpice            ! array used for the calculation of ice surface slope: 
    155                                                               !   ocean surface (ssh_m) if ice is not embedded 
    156                                                               !   ice top surface if ice is embedded    
    157148      !!------------------------------------------------------------------- 
    158149 
    … …  
    160151      CALL wrk_alloc( jpi,jpj, zc1   , u_oce1, u_oce2, u_ice2, zusw  , v_oce1 , v_oce2, v_ice1                ) 
    161152      CALL wrk_alloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds  , zdst  ) 
    162       CALL wrk_alloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 ) 
     153      CALL wrk_alloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr                         ) 
    163154 
    164155#if  defined key_lim2 && ! defined key_lim2_vp 
    … …  
    171162# endif 
    172163     at_i(:,:) = 1. - frld(:,:) 
    173 #endif 
    174 #if defined key_agrif && defined key_lim2  
    175     CALL agrif_rhg_lim2_load      ! First interpolation of coarse values 
    176164#endif 
    177165      ! 
    … …  
    244232      !  v_oce2: ocean v component on v points                         
    245233 
    246       IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==! 
    247           !                                             
    248           ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 
    249           !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1} 
    250          zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp      
    251           ! 
    252           ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 
    253           !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 
    254          zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 
    255           ! 
    256          zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0 
    257           ! 
    258       ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==! 
    259          zpice(:,:) = ssh_m(:,:) 
    260       ENDIF 
    261  
    262234      DO jj = k_j1+1, k_jpj-1 
    263235         DO ji = fs_2, fs_jpim1 
    … …  
    302274            ! include it later 
    303275 
    304             zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj) 
    305             zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj) 
     276            zdsshx =  ( ssh_m(ji+1,jj) - ssh_m(ji,jj) ) / e1u(ji,jj) 
     277            zdsshy =  ( ssh_m(ji,jj+1) - ssh_m(ji,jj) ) / e2v(ji,jj) 
    306278 
    307279            za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx 
    … …  
    520492 
    521493            CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
    522 #if defined key_agrif 
    523             CALL agrif_rhg_lim2( jter, nevp, 'U' ) 
    524 #endif 
    525494 
    526495!CDIR NOVERRCHK 
    … …  
    548517 
    549518            CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
    550 #if defined key_agrif 
    551             CALL agrif_rhg_lim2( jter, nevp, 'V' ) 
    552 #endif 
    553519 
    554520         ELSE  
    … …  
    577543 
    578544            CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
    579 #if defined key_agrif 
    580             CALL agrif_rhg_lim2( jter, nevp , 'V' ) 
    581 #endif 
    582545 
    583546!CDIR NOVERRCHK 
    … …  
    608571 
    609572            CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
    610 #if defined key_agrif 
    611             CALL agrif_rhg_lim2( jter, nevp, 'U' ) 
    612 #endif 
    613573 
    614574         ENDIF 
    … …  
    651611      CALL lbc_lnk( u_ice(:,:), 'U', -1. )  
    652612      CALL lbc_lnk( v_ice(:,:), 'V', -1. )  
    653 #if defined key_agrif 
    654       CALL agrif_rhg_lim2( nevp , nevp, 'U' ) 
    655       CALL agrif_rhg_lim2( nevp , nevp, 'V' ) 
    656 #endif 
    657613 
    658614      DO jj = k_j1+1, k_jpj-1  
    … …  
    790746      CALL wrk_dealloc( jpi,jpj, zc1   , u_oce1, u_oce2, u_ice2, zusw  , v_oce1 , v_oce2, v_ice1                ) 
    791747      CALL wrk_dealloc( jpi,jpj, zf1   , deltat, zu_ice, zf2   , deltac, zv_ice , zdd   , zdt    , zds  , zdst  ) 
    792       CALL wrk_dealloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 ) 
     748      CALL wrk_dealloc( jpi,jpj, zdd   , zdt   , zds   , zs1   , zs2   , zs12   , zresr                         ) 
    793749 
    794750   END SUBROUTINE lim_rhg 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r3785 r6736  
    1414 
    1515   !!---------------------------------------------------------------------- 
    16    !!   'key_asminc'   : Switch on the assimilation increment interface 
     16   !!   'key_asminc' : Switch on the assimilation increment interface 
    1717   !!---------------------------------------------------------------------- 
    18    !!   asm_inc_init   : Initialize the increment arrays and IAU weights 
    19    !!   calc_date      : Compute the calendar date YYYYMMDD on a given step 
    20    !!   tra_asm_inc    : Apply the tracer (T and S) increments 
    21    !!   dyn_asm_inc    : Apply the dynamic (u and v) increments 
    22    !!   ssh_asm_inc    : Apply the SSH increment 
    23    !!   seaice_asm_inc : Apply the seaice increment 
     18   !!   asm_inc_init : Initialize the increment arrays and IAU weights 
     19   !!   calc_date    : Compute the calendar date YYYYMMDD on a given step 
     20   !!   tra_asm_inc  : Apply the tracer (T and S) increments 
     21   !!   dyn_asm_inc  : Apply the dynamic (u and v) increments 
     22   !!   ssh_asm_inc  : Apply the SSH increment 
     23   !!   seaice_asm_inc  : Apply the seaice increment 
    2424   !!---------------------------------------------------------------------- 
    2525   USE wrk_nemo         ! Memory Allocation 
    2626   USE par_oce          ! Ocean space and time domain variables 
    2727   USE dom_oce          ! Ocean space and time domain 
    28    USE domvvl           ! domain: variable volume level 
    2928   USE oce              ! Dynamics and active tracers defined in memory 
    3029   USE ldfdyn_oce       ! ocean dynamics: lateral physics 
    … …  
    4039#endif 
    4140   USE sbc_oce          ! Surface boundary condition variables. 
     41   USE domvvl 
    4242 
    4343   IMPLICIT NONE 
    … …  
    9292#  include "ldfdyn_substitute.h90" 
    9393#  include "vectopt_loop_substitute.h90" 
     94 
    9495   !!---------------------------------------------------------------------- 
    9596   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    … …  
    109110      !! ** Action  :  
    110111      !!---------------------------------------------------------------------- 
    111       INTEGER :: ji, jj, jk 
     112      !! 
     113      !! 
     114      INTEGER :: ji,jj,jk 
    112115      INTEGER :: jt 
    113116      INTEGER :: imid 
    … …  
    939942            ! Update before fields 
    940943            sshb(:,:) = sshn(:,:)          
    941  
     944#if ! defined key_jth_fix 
    942945            IF( lk_vvl ) THEN 
    943946               DO jk = 1, jpk 
    … …  
    945948               END DO 
    946949            ENDIF 
    947  
     950#endif 
    948951            DEALLOCATE( ssh_bkg    ) 
    949952            DEALLOCATE( ssh_bkginc ) 
    … …  
    955958   END SUBROUTINE ssh_asm_inc 
    956959 
    957  
    958960   SUBROUTINE seaice_asm_inc( kt, kindic ) 
    959961      !!---------------------------------------------------------------------- 
    … …  
    966968      !! ** Action  :  
    967969      !! 
    968       !!---------------------------------------------------------------------- 
     970      !! History : 
     971      !!        !  07-2011  (D. Lea)  Initial version based on ssh_asm_inc 
     972      !!---------------------------------------------------------------------- 
     973 
    969974      IMPLICIT NONE 
    970       ! 
    971       INTEGER, INTENT(in)           ::   kt   ! Current time step 
    972       INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation 
    973       ! 
    974       INTEGER  ::   it 
    975       REAL(wp) ::   zincwgt   ! IAU weight for current time step 
     975 
     976      !! * Arguments 
     977      INTEGER, INTENT(IN) :: kt   ! Current time step 
     978      INTEGER, OPTIONAL, INTENT(IN) :: kindic ! flag for disabling the deallocation 
     979 
     980      !! * Local declarations 
     981      INTEGER :: it 
     982      REAL(wp) :: zincwgt  ! IAU weight for current time step 
     983 
    976984#if defined key_lim2 
    977       REAL(wp), DIMENSION(jpi,jpj) ::   zofrld, zohicif, zseaicendg, zhicifinc  ! LIM 
    978       REAL(wp) ::   zhicifmin = 0.5_wp      ! ice minimum depth in metres 
    979 #endif 
    980       !!---------------------------------------------------------------------- 
     985      REAL(wp), DIMENSION(jpi,jpj) :: zofrld, zohicif, zseaicendg, zhicifinc  ! LIM 
     986      REAL(wp) :: zhicifmin=0.5_wp      ! ice minimum depth in metres 
     987 
     988#endif 
     989 
    981990 
    982991      IF ( ln_asmiau ) THEN 
    … …  
    9991008            ENDIF 
    10001009 
    1001             ! Sea-ice : LIM-3 case (to add) 
    1002  
    10031010#if defined key_lim2 
    1004             ! Sea-ice : LIM-2 case 
    1005             zofrld (:,:) = frld(:,:) 
    1006             zohicif(:,:) = hicif(:,:) 
    1007             ! 
    1008             frld  = MIN( MAX( frld (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 
     1011 
     1012            zofrld(:,:)=frld(:,:) 
     1013            zohicif(:,:)=hicif(:,:) 
     1014 
     1015            frld = MIN( MAX( frld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 
    10091016            pfrld = MIN( MAX( pfrld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 
    10101017            fr_i(:,:) = 1.0_wp - frld(:,:)        ! adjust ice fraction 
    1011             ! 
    1012             zseaicendg(:,:) = zofrld(:,:) - frld(:,:)   ! find out actual sea ice nudge applied 
    1013             ! 
     1018 
     1019            zseaicendg(:,:)=zofrld(:,:) - frld(:,:)         ! find out actual sea ice nudge applied 
     1020 
    10141021            ! Nudge sea ice depth to bring it up to a required minimum depth 
     1022 
    10151023            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )  
    10161024               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt     
    … …  
    10181026               zhicifinc(:,:) = 0.0_wp 
    10191027            END WHERE 
    1020             ! 
    1021             ! nudge ice depth 
    1022             hicif (:,:) = hicif (:,:) + zhicifinc(:,:) 
    1023             phicif(:,:) = phicif(:,:) + zhicifinc(:,:)        
    1024             ! 
    1025             ! seaice salinity balancing (to add) 
     1028 
     1029! nudge ice depth 
     1030            hicif(:,:)=hicif(:,:) + zhicifinc(:,:) 
     1031            phicif(:,:)=phicif(:,:) + zhicifinc(:,:)        
     1032 
     1033! seaice salinity balancing (to add) 
     1034 
    10261035#endif 
    10271036 
    10281037#if defined key_cice 
    1029             ! Sea-ice : CICE case. Pass ice increment tendency into CICE 
     1038 
     1039! Pass ice increment tendency into CICE 
    10301040            ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt 
     1041 
    10311042#endif 
    10321043 
    … …  
    10381049 
    10391050#if defined key_cice 
    1040             ! Sea-ice : CICE case. Zero ice increment tendency into CICE 
     1051 
     1052! Zero ice increment tendency into CICE 
    10411053            ndaice_da(:,:) = 0.0_wp 
     1054 
    10421055#endif 
    10431056 
    … …  
    10541067            neuler = 0                    ! Force Euler forward step 
    10551068 
    1056             ! Sea-ice : LIM-3 case (to add) 
    1057  
    10581069#if defined key_lim2 
    1059             ! Sea-ice : LIM-2 case. 
     1070 
    10601071            zofrld(:,:)=frld(:,:) 
    10611072            zohicif(:,:)=hicif(:,:) 
    1062             !  
     1073  
    10631074            ! Initialize the now fields the background + increment 
    1064             frld (:,:) = MIN( MAX( frld(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp) 
     1075 
     1076            frld(:,:) = MIN( MAX( frld(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp) 
    10651077            pfrld(:,:) = frld(:,:)  
    1066             fr_i (:,:) = 1.0_wp - frld(:,:)                ! adjust ice fraction 
    1067             zseaicendg(:,:) = zofrld(:,:) - frld(:,:)      ! find out actual sea ice nudge applied 
    1068             ! 
     1078            fr_i(:,:) = 1.0_wp - frld(:,:)        ! adjust ice fraction 
     1079 
     1080            zseaicendg(:,:)=zofrld(:,:) - frld(:,:)         ! find out actual sea ice nudge applied 
     1081 
    10691082            ! Nudge sea ice depth to bring it up to a required minimum depth 
     1083 
    10701084            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )  
    10711085               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt     
    … …  
    10731087               zhicifinc(:,:) = 0.0_wp 
    10741088            END WHERE 
    1075             ! 
    1076             ! nudge ice depth 
    1077             hicif (:,:) = hicif (:,:) + zhicifinc(:,:) 
    1078             phicif(:,:) = phicif(:,:)        
    1079             ! 
    1080             ! seaice salinity balancing (to add) 
     1089 
     1090! nudge ice depth 
     1091            hicif(:,:)=hicif(:,:) + zhicifinc(:,:) 
     1092            phicif(:,:)=phicif(:,:)        
     1093 
     1094! seaice salinity balancing (to add) 
     1095   
    10811096#endif 
    10821097  
    10831098#if defined key_cice 
    1084             ! Sea-ice : CICE case. Pass ice increment tendency into CICE - is this correct? 
     1099 
     1100! Pass ice increment tendency into CICE - is this correct? 
    10851101           ndaice_da(:,:) = seaice_bkginc(:,:) / rdt 
     1102 
    10861103#endif 
    10871104           IF ( .NOT. PRESENT(kindic) ) THEN 
    … …  
    10921109 
    10931110#if defined key_cice 
    1094             ! Sea-ice : CICE case. Zero ice increment tendency into CICE  
     1111 
     1112! Zero ice increment tendency into CICE  
    10951113            ndaice_da(:,:) = 0.0_wp 
     1114 
    10961115#endif 
    10971116          
    10981117         ENDIF 
    10991118 
    1100 !#if defined defined key_lim2 || defined key_cice 
     1119!#if defined key_lim2 || defined key_cice 
    11011120! 
    11021121!            IF (ln_seaicebal ) THEN        
    … …  
    11731192!#endif 
    11741193 
     1194 
    11751195      ENDIF 
    11761196 
    11771197   END SUBROUTINE seaice_asm_inc 
    1178     
    11791198   !!====================================================================== 
    11801199END MODULE asminc 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r3651 r6736  
    2828      INTEGER, POINTER, DIMENSION(:,:)   ::  nbmap 
    2929      REAL   , POINTER, DIMENSION(:,:)   ::  nbw 
    30       REAL   , POINTER, DIMENSION(:,:)   ::  nbd 
     30      REAL   , POINTER, DIMENSION(:,:,:)   ::  nbz 
    3131      REAL   , POINTER, DIMENSION(:)     ::  flagu 
    3232      REAL   , POINTER, DIMENSION(:)     ::  flagv 
    … …  
    4646      REAL, POINTER, DIMENSION(:)     ::  hsnif 
    4747#endif 
    48    END TYPE OBC_DATA 
     48   END TYPE OBC_DATA  
    4949 
    5050   !!---------------------------------------------------------------------- 
    … …  
    7474   INTEGER, DIMENSION(jp_bdy) ::   nn_tra_dta             !: = 0 use the initial state as bdy dta ;  
    7575                                                            !: = 1 read it in a NetCDF file 
    76    LOGICAL, DIMENSION(jp_bdy) ::   ln_tra_dmp               !: =T Tracer damping 
    77    LOGICAL, DIMENSION(jp_bdy) ::   ln_dyn3d_dmp             !: =T Baroclinic velocity damping 
    78    REAL,    DIMENSION(jp_bdy) ::   rn_time_dmp              !: Damping time scale in days 
    79  
     76   INTEGER                    ::   nb_jpk              ! Number of levels in the bdy data (set < 0 if consistent with planned run) 
    8077#if defined key_lim2 
    8178   INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim2              ! Choice of boundary condition for sea ice variables  
    … …  
    106103   INTEGER,  DIMENSION(jp_bdy)                     ::   nn_dta            !: =0 => *all* data is set to initial conditions 
    107104                                                                          !: =1 => some data to be read in from data files 
    108    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays (unstr.  bdy) 
    109    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2       !: workspace for reading in global data arrays (struct. bdy) 
     105   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays 
     106   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global_1      !: workspace for reading in global data arrays 
     107   REAL(wp), ALLOCATABLE, DIMENSION(:,:  ), TARGET ::   dta_global_2      !: workspace for reading in global data arrays 
    110108   TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process) 
    111109   TYPE(OBC_DATA) , DIMENSION(jp_bdy)              ::   dta_bdy           !: bdy external data (local process) 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r3851 r6736  
    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.4  !  2013-04  (J. Harle) add in option to read bdy data with 
     14   !!                                        different vertical coordinates 
    1315   !!---------------------------------------------------------------------- 
    1416#if defined key_bdy 
    … …  
    3234   USE ice_2 
    3335#endif 
    34    USE sbcapr 
    3536 
    3637   IMPLICIT NONE 
    … …  
    109110 
    110111            IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 0 ) THEN  
    111                ilen1(:) = nblen(:) 
     112               IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     113                  ilen1(:) = nblen(:) 
     114               ELSE 
     115                  ilen1(:) = nblenrim(:) 
     116               ENDIF 
    112117               igrd = 1 
    113118               DO ib = 1, ilen1(igrd) 
    … …  
    131136 
    132137            IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN  
    133                ilen1(:) = nblen(:) 
     138               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 
     139                  ilen1(:) = nblen(:) 
     140               ELSE 
     141                  ilen1(:) = nblenrim(:) 
     142               ENDIF 
    134143               igrd = 2  
    135144               DO ib = 1, ilen1(igrd) 
    … …  
    151160 
    152161            IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 0 ) THEN  
    153                ilen1(:) = nblen(:) 
     162               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 
     163                  ilen1(:) = nblen(:) 
     164               ELSE 
     165                  ilen1(:) = nblenrim(:) 
     166               ENDIF 
    154167               igrd = 1                       ! Everything is at T-points here 
    155168               DO ib = 1, ilen1(igrd) 
    … …  
    165178#if defined key_lim2 
    166179            IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN  
    167                ilen1(:) = nblen(:) 
     180               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
     181                  ilen1(:) = nblen(:) 
     182               ELSE 
     183                  ilen1(:) = nblenrim(:) 
     184               ENDIF 
    168185               igrd = 1                       ! Everything is at T-points here 
    169186               DO ib = 1, ilen1(igrd) 
    … …  
    192209            IF( PRESENT(jit) ) THEN 
    193210               ! Update barotropic boundary conditions only 
    194                ! jit is optional argument for fld_read and bdytide_update 
     211               ! jit is optional argument for fld_read and tide_update 
    195212               IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN 
    196213                  IF( nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
    … …  
    199216                     dta_bdy(ib_bdy)%v2d(:) = 0.0 
    200217                  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 
    208                         CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  & 
    209                                      & 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 
    213                         IF( ln_full_vel_array(ib_bdy) .AND.                                             & 
    214                           &    ( nn_dyn2d_dta(ib_bdy) .EQ. 1 .OR. nn_dyn2d_dta(ib_bdy) .EQ. 3 .OR.  & 
    215                           &      nn_dyn3d_dta(ib_bdy) .EQ. 1 ) )THEN 
    216  
    217                            igrd = 2                      ! zonal velocity 
    218                            dta_bdy(ib_bdy)%u2d(:) = 0.0 
    219                            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    220                               ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    221                               ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    222                               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) 
    225                               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 
    230                            END DO 
    231                            igrd = 3                      ! meridional velocity 
    232                            dta_bdy(ib_bdy)%v2d(:) = 0.0 
    233                            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    234                               ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    235                               ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    236                               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) 
    239                               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 
    244                            END DO 
    245                         ENDIF                     
    246                      ENDIF 
    247                      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),   &  
    249                           &                 jit=jit, time_offset=time_offset ) 
    250                      ENDIF 
     218                  IF( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN ! update external data 
     219                     jend = jstart + 2 
     220                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),   & 
     221                     &              jit=jit, time_offset=time_offset ) 
    251222                  ENDIF 
     223                  IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
     224                     CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy),   &  
     225                     &                 jit=jit, time_offset=time_offset ) 
     226                  ENDIF 
    252227               ENDIF 
    253228            ELSE 
    254                IF (nn_tra(ib_bdy).eq.4) then      ! runoff condition 
    255                   jend = nb_bdy_fld(ib_bdy) 
    256                   CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend),  & 
    257                                & map=nbmap_ptr(jstart:jend), kt_offset=time_offset ) 
    258                   ! 
    259                   igrd = 2                      ! zonal velocity 
    260                   DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    261                      ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    262                      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) ) 
     229               IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
     230                  dta_bdy(ib_bdy)%ssh(:) = 0.0 
     231                  dta_bdy(ib_bdy)%u2d(:) = 0.0 
     232                  dta_bdy(ib_bdy)%v2d(:) = 0.0 
     233               ENDIF 
     234               IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data 
     235                  jend = jstart + nb_bdy_fld(ib_bdy) - 1 
     236                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset,& 
     237                &                                                                                          jpk_1=nb_jpk ) 
     238               ENDIF 
     239               IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
     240                  CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), time_offset=time_offset ) 
     241               ENDIF 
     242            ENDIF 
     243            jstart = jend+1 
     244 
     245            ! If full velocities in boundary data then split into barotropic and baroclinic data 
     246            ! (Note that we have already made sure that you can't use ln_full_vel = .true. at the same 
     247            ! time as the dynspg_ts option).  
     248 
     249            IF( ln_full_vel_array(ib_bdy) .and.                                             &  
     250           &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN  
     251 
     252               igrd = 2                      ! zonal velocity 
     253               dta_bdy(ib_bdy)%u2d(:) = 0.0 
     254               DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     255                  ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     256                  ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     257                  DO ik = 1, jpkm1 
     258                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) & 
     259              &                                + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik) 
    264260                  END DO 
    265                   ! 
    266                   igrd = 3                      ! meridional velocity 
    267                   DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    268                      ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    269                      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) ) 
     261                  dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 
     262                  DO ik = 1, jpkm1 
     263                     dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib)  
    271264                  END DO 
    272                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 
    277                   ENDIF 
    278                   IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data 
    279                      jend = nb_bdy_fld(ib_bdy) 
    280                      CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), & 
    281                                   & map=nbmap_ptr(jstart:jend), kt_offset=time_offset ) 
    282                   ENDIF 
    283                   ! If full velocities in boundary data then split into barotropic and baroclinic data 
    284                   IF( ln_full_vel_array(ib_bdy) .and.                                             & 
    285                     & ( nn_dyn2d_dta(ib_bdy) .EQ. 1 .OR. nn_dyn2d_dta(ib_bdy) .EQ. 3 .OR. & 
    286                     &   nn_dyn3d_dta(ib_bdy) .EQ. 1 ) ) THEN 
    287                      igrd = 2                      ! zonal velocity 
    288                      dta_bdy(ib_bdy)%u2d(:) = 0.0 
    289                      DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    290                         ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    291                         ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    292                         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) 
    295                         END DO 
    296                         dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 
    297                         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) 
    299                         END DO 
    300                      END DO 
    301                      igrd = 3                      ! meridional velocity 
    302                      dta_bdy(ib_bdy)%v2d(:) = 0.0 
    303                      DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    304                         ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    305                         ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    306                         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) 
    309                         END DO 
    310                         dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 
    311                         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) 
    313                         END DO 
    314                      END DO 
    315                   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 
     265               END DO 
     266 
     267               igrd = 3                      ! meridional velocity 
     268               dta_bdy(ib_bdy)%v2d(:) = 0.0 
     269               DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
     270                  ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     271                  ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     272                  DO ik = 1, jpkm1 
     273                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) & 
     274              &                                + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik) 
     275                  END DO 
     276                  dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 
     277                  DO ik = 1, jpkm1 
     278                     dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib)  
     279                  END DO 
     280               END DO 
     281     
     282            ENDIF 
     283 
    323284         END IF ! nn_dta(ib_bdy) = 1 
    324285      END DO  ! ib_bdy 
    325  
    326       IF ( ln_apr_obc ) THEN 
    327          DO ib_bdy = 1, nb_bdy 
    328             IF (nn_tra(ib_bdy).NE.4)THEN 
    329                igrd = 1                      ! meridional velocity 
    330                DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    331                   ii   = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    332                   ij   = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    333                   dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + ssh_ib(ii,ij) 
    334                ENDDO 
    335             ENDIF 
    336          ENDDO 
    337       ENDIF 
    338286 
    339287      IF( nn_timing == 1 ) CALL timing_stop('bdy_dta') 
    … …  
    381329      IF( nn_timing == 1 ) CALL timing_start('bdy_dta_init') 
    382330 
    383       IF(lwp) WRITE(numout,*) 
    384       IF(lwp) WRITE(numout,*) 'bdy_dta_ini : initialization of data at the open boundaries' 
    385       IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
    386       IF(lwp) WRITE(numout,*) '' 
    387  
    388331      ! Set nn_dta 
    389332      DO ib_bdy = 1, nb_bdy 
    … …  
    417360         ENDIF 
    418361#endif                
    419          IF(lwp) WRITE(numout,*) 'Maximum number of files to open =',nb_bdy_fld(ib_bdy) 
    420362      ENDDO             
    421363 
    … …  
    469411            ln_full_vel_array(ib_bdy) = ln_full_vel 
    470412 
     413            IF( ln_full_vel_array(ib_bdy) .and. lk_dynspg_ts )  THEN 
     414               CALL ctl_stop( 'bdy_dta_init: ERROR, cannot specify full velocities in boundary data',& 
     415            &                  'with dynspg_ts option' )   ;   RETURN   
     416            ENDIF              
     417 
    471418            nblen => idx_bdy(ib_bdy)%nblen 
    472419            nblenrim => idx_bdy(ib_bdy)%nblenrim 
    … …  
    476423            IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN  
    477424 
    478                IF( nn_tra(ib_bdy) .ne. 4 ) THEN ! runoff condition : no ssh reading 
     425               IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN 
    479426                  jfld = jfld + 1 
    480427                  blf_i(jfld) = bn_ssh 
    481428                  ibdy(jfld) = ib_bdy 
    482429                  igrid(jfld) = 1 
    483                   ilen1(jfld) = nblen(igrid(jfld)) 
     430                  ilen1(jfld) = nblenrim(igrid(jfld)) 
    484431                  ilen3(jfld) = 1 
    485432               ENDIF 
    486433 
    487434               IF( .not. ln_full_vel_array(ib_bdy) ) THEN 
     435 
    488436                  jfld = jfld + 1 
    489437                  blf_i(jfld) = bn_u2d 
    490438                  ibdy(jfld) = ib_bdy 
    491439                  igrid(jfld) = 2 
    492                   ilen1(jfld) = nblen(igrid(jfld)) 
     440                  IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     441                     ilen1(jfld) = nblen(igrid(jfld)) 
     442                  ELSE 
     443                     ilen1(jfld) = nblenrim(igrid(jfld)) 
     444                  ENDIF 
    493445                  ilen3(jfld) = 1 
    494446 
    … …  
    497449                  ibdy(jfld) = ib_bdy 
    498450                  igrid(jfld) = 3 
    499                   ilen1(jfld) = nblen(igrid(jfld)) 
     451                  IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     452                     ilen1(jfld) = nblen(igrid(jfld)) 
     453                  ELSE 
     454                     ilen1(jfld) = nblenrim(igrid(jfld)) 
     455                  ENDIF 
    500456                  ilen3(jfld) = 1 
     457 
    501458               ENDIF 
    502459 
    … …  
    512469               ibdy(jfld) = ib_bdy 
    513470               igrid(jfld) = 2 
    514                ilen1(jfld) = nblen(igrid(jfld)) 
     471               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 
     472                  ilen1(jfld) = nblen(igrid(jfld)) 
     473               ELSE 
     474                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     475               ENDIF 
    515476               ilen3(jfld) = jpk 
    516477 
    … …  
    519480               ibdy(jfld) = ib_bdy 
    520481               igrid(jfld) = 3 
    521                ilen1(jfld) = nblen(igrid(jfld)) 
     482               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 
     483                  ilen1(jfld) = nblen(igrid(jfld)) 
     484               ELSE 
     485                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     486               ENDIF 
    522487               ilen3(jfld) = jpk 
    523488 
    … …  
    531496               ibdy(jfld) = ib_bdy 
    532497               igrid(jfld) = 1 
    533                ilen1(jfld) = nblen(igrid(jfld)) 
     498               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 
     499                  ilen1(jfld) = nblen(igrid(jfld)) 
     500               ELSE 
     501                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     502               ENDIF 
    534503               ilen3(jfld) = jpk 
    535504 
    … …  
    538507               ibdy(jfld) = ib_bdy 
    539508               igrid(jfld) = 1 
    540                ilen1(jfld) = nblen(igrid(jfld)) 
     509               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 
     510                  ilen1(jfld) = nblen(igrid(jfld)) 
     511               ELSE 
     512                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     513               ENDIF 
    541514               ilen3(jfld) = jpk 
    542515 
    … …  
    551524               ibdy(jfld) = ib_bdy 
    552525               igrid(jfld) = 1 
    553                ilen1(jfld) = nblen(igrid(jfld)) 
     526               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
     527                  ilen1(jfld) = nblen(igrid(jfld)) 
     528               ELSE 
     529                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     530               ENDIF 
    554531               ilen3(jfld) = 1 
    555532 
    … …  
    558535               ibdy(jfld) = ib_bdy 
    559536               igrid(jfld) = 1 
    560                ilen1(jfld) = nblen(igrid(jfld)) 
     537               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
     538                  ilen1(jfld) = nblen(igrid(jfld)) 
     539               ELSE 
     540                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     541               ENDIF 
    561542               ilen3(jfld) = 1 
    562543 
    … …  
    565546               ibdy(jfld) = ib_bdy 
    566547               igrid(jfld) = 1 
    567                ilen1(jfld) = nblen(igrid(jfld)) 
     548               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
     549                  ilen1(jfld) = nblen(igrid(jfld)) 
     550               ELSE 
     551                  ilen1(jfld) = nblenrim(igrid(jfld)) 
     552               ENDIF 
    568553               ilen3(jfld) = 1 
    569554 
    … …  
    584569      ENDDO ! ib_bdy 
    585570 
     571 
    586572      DO jfld = 1, nb_bdy_fld_sum 
    587573         ALLOCATE( bf(jfld)%fnow(ilen1(jfld),1,ilen3(jfld)) ) 
    … …  
    594580      jstart = 1 
    595581      DO ib_bdy = 1, nb_bdy 
    596          jend = nb_bdy_fld(ib_bdy)  
     582         jend = jstart + nb_bdy_fld(ib_bdy) - 1 
    597583         CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(ib_bdy), 'bdy_dta',   & 
    598584         &              'open boundary conditions', 'nambdy_dta' ) 
    … …  
    613599         IF (nn_dyn2d(ib_bdy) .gt. 0) THEN 
    614600            IF( nn_dyn2d_dta(ib_bdy) .eq. 0 .or. nn_dyn2d_dta(ib_bdy) .eq. 2 .or. ln_full_vel_array(ib_bdy) ) THEN 
    615                ilen0(1:3) = nblen(1:3) 
     601               IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     602                  ilen0(1:3) = nblen(1:3) 
     603               ELSE 
     604                  ilen0(1:3) = nblenrim(1:3) 
     605               ENDIF 
     606               ALLOCATE( dta_bdy(ib_bdy)%ssh(ilen0(1)) ) 
    616607               ALLOCATE( dta_bdy(ib_bdy)%u2d(ilen0(2)) ) 
    617608               ALLOCATE( dta_bdy(ib_bdy)%v2d(ilen0(3)) ) 
    618                IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3) THEN 
    619                   jfld = jfld + 1 
    620                   dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1) 
    621                ELSE 
    622                   ALLOCATE( dta_bdy(ib_bdy)%ssh(nblen(1)) ) 
    623                ENDIF 
    624609            ELSE 
    625610               IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN 
    … …  
    635620 
    636621         IF ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 
    637             ilen0(1:3) = nblen(1:3) 
     622            IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 
     623               ilen0(1:3) = nblen(1:3) 
     624            ELSE 
     625               ilen0(1:3) = nblenrim(1:3) 
     626            ENDIF 
    638627            ALLOCATE( dta_bdy(ib_bdy)%u3d(ilen0(2),jpk) ) 
    639628            ALLOCATE( dta_bdy(ib_bdy)%v3d(ilen0(3),jpk) ) 
    … …  
    650639         IF (nn_tra(ib_bdy) .gt. 0) THEN 
    651640            IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN 
    652                ilen0(1:3) = nblen(1:3) 
     641               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 
     642                  ilen0(1:3) = nblen(1:3) 
     643               ELSE 
     644                  ilen0(1:3) = nblenrim(1:3) 
     645               ENDIF 
    653646               ALLOCATE( dta_bdy(ib_bdy)%tem(ilen0(1),jpk) ) 
    654647               ALLOCATE( dta_bdy(ib_bdy)%sal(ilen0(1),jpk) ) 
    … …  
    664657         IF (nn_ice_lim2(ib_bdy) .gt. 0) THEN 
    665658            IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 
    666                ilen0(1:3) = nblen(1:3) 
     659               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 
     660                  ilen0(1:3) = nblen(1:3) 
     661               ELSE 
     662                  ilen0(1:3) = nblenrim(1:3) 
     663               ENDIF 
    667664               ALLOCATE( dta_bdy(ib_bdy)%frld(ilen0(1)) ) 
    668665               ALLOCATE( dta_bdy(ib_bdy)%hicif(ilen0(1)) ) 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90

    r3680 r6736  
    55   !!====================================================================== 
    66   !! History :  3.4  !  2011     (D. Storkey) new module as part of BDY rewrite 
    7    !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Optimization of BDY communications 
    87   !!---------------------------------------------------------------------- 
    98#if defined key_bdy  
    … …  
    5251            CYCLE 
    5352         CASE(jp_frs) 
    54             CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy ) 
     53            CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy) ) 
    5554         CASE(jp_flather) 
    56             CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy ) 
     55            CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy) ) 
    5756         CASE DEFAULT 
    5857            CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' ) 
    … …  
    6261   END SUBROUTINE bdy_dyn2d 
    6362 
    64    SUBROUTINE bdy_dyn2d_frs( idx, dta, ib_bdy ) 
     63   SUBROUTINE bdy_dyn2d_frs( idx, dta ) 
    6564      !!---------------------------------------------------------------------- 
    6665      !!                  ***  SUBROUTINE bdy_dyn2d_frs  *** 
    … …  
    7574      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
    7675      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
    77       INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index 
    7876      !! 
    7977      INTEGER  ::   jb, jk         ! dummy loop indices 
    … …  
    9997         pv2d(ii,ij) = ( pv2d(ii,ij) + zwgt * ( dta%v2d(jb) - pv2d(ii,ij) ) ) * vmask(ii,ij,1) 
    10098      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 
     99      CALL lbc_lnk( pu2d, 'U', -1. )  
     100      CALL lbc_lnk( pv2d, 'V', -1. )   ! Boundary points should be updated 
    103101      ! 
    104102      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_frs') 
    … …  
    108106 
    109107 
    110    SUBROUTINE bdy_dyn2d_fla( idx, dta, ib_bdy ) 
     108   SUBROUTINE bdy_dyn2d_fla( idx, dta ) 
    111109      !!---------------------------------------------------------------------- 
    112110      !!                 ***  SUBROUTINE bdy_dyn2d_fla  *** 
    … …  
    129127      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices 
    130128      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data 
    131       INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index 
    132129 
    133130      INTEGER  ::   jb, igrd                         ! dummy loop indices 
    … …  
    180177         pv2d(ii,ij) = zforc + zcorr * vmask(ii,ij,1) 
    181178      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 )   ! 
     179      CALL lbc_lnk( pu2d, 'U', -1. )   ! Boundary points should be updated 
     180      CALL lbc_lnk( pv2d, 'V', -1. )   ! 
    184181      ! 
    185182      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_fla') 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90

    r3703 r6736  
    55   !!====================================================================== 
    66   !! History :  3.4  !  2011     (D. Storkey) new module as part of BDY rewrite  
    7    !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Optimization of BDY communications 
    87   !!---------------------------------------------------------------------- 
    98#if defined key_bdy  
    … …  
    1514   !!---------------------------------------------------------------------- 
    1615   USE timing          ! Timing 
    17    USE wrk_nemo        ! Memory Allocation 
    1816   USE oce             ! ocean dynamics and tracers  
    1917   USE dom_oce         ! ocean space and time domain 
    … …  
    2119   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2220   USE in_out_manager  ! 
    23    Use phycst 
    2421 
    2522   IMPLICIT NONE 
    … …  
    2724 
    2825   PUBLIC   bdy_dyn3d     ! routine called by bdy_dyn 
    29    PUBLIC   bdy_dyn3d_dmp ! routine called by step 
    3026 
    31    !! * Substitutions 
    32 #  include "domzgr_substitute.h90" 
    3327   !!---------------------------------------------------------------------- 
    3428   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    … …  
    6054            CYCLE 
    6155         CASE(jp_frs) 
    62             CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
    63          CASE(2) 
    64             CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
    65          CASE(3) 
    66             CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
     56            CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 
    6757         CASE DEFAULT 
    6858            CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' ) 
    … …  
    7262   END SUBROUTINE bdy_dyn3d 
    7363 
    74    SUBROUTINE bdy_dyn3d_spe( idx, dta, kt , ib_bdy ) 
    75       !!---------------------------------------------------------------------- 
    76       !!                  ***  SUBROUTINE bdy_dyn3d_spe  *** 
    77       !! 
    78       !! ** Purpose : - Apply a specified value for baroclinic velocities 
    79       !!                at open boundaries. 
    80       !! 
    81       !!---------------------------------------------------------------------- 
    82       INTEGER                     ::   kt 
    83       TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
    84       TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
    85       INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index 
    86       !! 
    87       INTEGER  ::   jb, jk         ! dummy loop indices 
    88       INTEGER  ::   ii, ij, igrd   ! local integers 
    89       REAL(wp) ::   zwgt           ! boundary weight 
    90       !!---------------------------------------------------------------------- 
    91       ! 
    92       IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_spe') 
    93       ! 
    94       igrd = 2                      ! Relaxation of zonal velocity 
    95       DO jb = 1, idx%nblenrim(igrd) 
    96          DO jk = 1, jpkm1 
    97             ii   = idx%nbi(jb,igrd) 
    98             ij   = idx%nbj(jb,igrd) 
    99             ua(ii,ij,jk) = dta%u3d(jb,jk) * umask(ii,ij,jk) 
    100          END DO 
    101       END DO 
    102       ! 
    103       igrd = 3                      ! Relaxation of meridional velocity 
    104       DO jb = 1, idx%nblenrim(igrd) 
    105          DO jk = 1, jpkm1 
    106             ii   = idx%nbi(jb,igrd) 
    107             ij   = idx%nbj(jb,igrd) 
    108             va(ii,ij,jk) = dta%v3d(jb,jk) * vmask(ii,ij,jk) 
    109          END DO 
    110       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       ! 
    113       IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
    114  
    115       IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_spe') 
    116  
    117    END SUBROUTINE bdy_dyn3d_spe 
    118  
    119    SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy ) 
    120       !!---------------------------------------------------------------------- 
    121       !!                  ***  SUBROUTINE bdy_dyn3d_zro  *** 
    122       !! 
    123       !! ** Purpose : - baroclinic velocities = 0. at open boundaries. 
    124       !! 
    125       !!---------------------------------------------------------------------- 
    126       INTEGER                     ::   kt 
    127       TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
    128       TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
    129       INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index 
    130       !! 
    131       INTEGER  ::   ib, ik         ! dummy loop indices 
    132       INTEGER  ::   ii, ij, igrd, zcoef   ! local integers 
    133       REAL(wp) ::   zwgt           ! boundary weight 
    134       !!---------------------------------------------------------------------- 
    135       ! 
    136       IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zro') 
    137       ! 
    138       igrd = 2                       ! Everything is at T-points here 
    139       DO ib = 1, idx%nblenrim(igrd) 
    140          ii = idx%nbi(ib,igrd) 
    141          ij = idx%nbj(ib,igrd) 
    142          DO ik = 1, jpkm1 
    143             ua(ii,ij,ik) = 0._wp 
    144          END DO 
    145       END DO 
    146  
    147       igrd = 3                       ! Everything is at T-points here 
    148       DO ib = 1, idx%nblenrim(igrd) 
    149          ii = idx%nbi(ib,igrd) 
    150          ij = idx%nbj(ib,igrd) 
    151          DO ik = 1, jpkm1 
    152             va(ii,ij,ik) = 0._wp 
    153          END DO 
    154       END DO 
    155       ! 
    156       CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ;   CALL lbc_bdy_lnk( va, 'V', -1.,ib_bdy )   ! Boundary points should be updated 
    157       ! 
    158       IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
    159  
    160       IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zro') 
    161  
    162    END SUBROUTINE bdy_dyn3d_zro 
    163  
    164    SUBROUTINE bdy_dyn3d_frs( idx, dta, kt, ib_bdy ) 
     64   SUBROUTINE bdy_dyn3d_frs( idx, dta, kt ) 
    16565      !!---------------------------------------------------------------------- 
    16666      !!                  ***  SUBROUTINE bdy_dyn3d_frs  *** 
    … …  
    17676      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
    17777      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
    178       INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index 
    17978      !! 
    18079      INTEGER  ::   jb, jk         ! dummy loop indices 
    … …  
    204103         END DO 
    205104      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 
     105      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated 
    207106      ! 
    208107      IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 
    … …  
    212111   END SUBROUTINE bdy_dyn3d_frs 
    213112 
    214    SUBROUTINE bdy_dyn3d_dmp( kt ) 
    215       !!---------------------------------------------------------------------- 
    216       !!                  ***  SUBROUTINE bdy_dyn3d_dmp  *** 
    217       !! 
    218       !! ** Purpose : Apply damping for baroclinic velocities at open boundaries. 
    219       !! 
    220       !!---------------------------------------------------------------------- 
    221       INTEGER                     ::   kt 
    222       !! 
    223       INTEGER  ::   jb, jk         ! dummy loop indices 
    224       INTEGER  ::   ii, ij, igrd   ! local integers 
    225       REAL(wp) ::   zwgt           ! boundary weight 
    226       INTEGER  ::  ib_bdy          ! loop index 
    227       !!---------------------------------------------------------------------- 
    228       ! 
    229       IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_dmp') 
    230       ! 
    231       !------------------------------------------------------- 
    232       ! Remove barotropic part from before velocity 
    233       !------------------------------------------------------- 
    234       CALL wrk_alloc(jpi,jpj,pu2d,pv2d)  
    235  
    236       pu2d(:,:) = 0.e0 
    237       pv2d(:,:) = 0.e0 
    238  
    239       DO jk = 1, jpkm1 
    240 #if defined key_vvl 
    241          pu2d(:,:) = pu2d(:,:) + fse3u_b(:,:,jk)* ub(:,:,jk)   *umask(:,:,jk)  
    242          pv2d(:,:) = pv2d(:,:) + fse3v_b(:,:,jk)* vb(:,:,jk)   *vmask(:,:,jk) 
    243 #else 
    244          pu2d(:,:) = pu2d(:,:) + fse3u_0(:,:,jk) * ub(:,:,jk)  * umask(:,:,jk) 
    245          pv2d(:,:) = pv2d(:,:) + fse3v_0(:,:,jk) * vb(:,:,jk)  * vmask(:,:,jk) 
    246 #endif 
    247       END DO 
    248  
    249       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) ) 
    252       ELSE 
    253          pu2d(:,:) = pv2d(:,:) * hur(:,:) 
    254          pv2d(:,:) = pu2d(:,:) * hvr(:,:) 
    255       ENDIF 
    256  
    257       DO ib_bdy=1, nb_bdy 
    258          IF ( ln_dyn3d_dmp(ib_bdy).and.nn_dyn3d(ib_bdy).gt.0 ) THEN 
    259             igrd = 2                      ! Relaxation of zonal velocity 
    260             DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    261                ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
    262                ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
    263                zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd) 
    264                DO jk = 1, jpkm1 
    265                   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) 
    267                END DO 
    268             END DO 
    269             ! 
    270             igrd = 3                      ! Relaxation of meridional velocity 
    271             DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    272                ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
    273                ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
    274                zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd) 
    275                DO jk = 1, jpkm1 
    276                   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) 
    278                END DO 
    279             END DO 
    280          ENDIF 
    281       ENDDO 
    282       ! 
    283       CALL wrk_dealloc(jpi,jpj,pu2d,pv2d)  
    284       ! 
    285       CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated 
    286       ! 
    287       IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_dmp') 
    288  
    289    END SUBROUTINE bdy_dyn3d_dmp 
    290113 
    291114#else 
    … …  
    295118CONTAINS 
    296119   SUBROUTINE bdy_dyn3d( kt )      ! Empty routine 
    297       WRITE(*,*) 'bdy_dyn3d: You should not have seen this print! error?', kt 
     120      WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt 
    298121   END SUBROUTINE bdy_dyn3d 
    299  
    300    SUBROUTINE bdy_dyn3d_dmp( kt )      ! Empty routine 
    301       WRITE(*,*) 'bdy_dyn3d_dmp: You should not have seen this print! error?', kt 
    302    END SUBROUTINE bdy_dyn3d_dmp 
    303  
    304122#endif 
    305123 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim2.F90

    r3680 r6736  
    66   !!  History :  3.3  !  2010-09 (D. Storkey)  Original code 
    77   !!             3.4  !  2011    (D. Storkey) rewrite in preparation for OBC-BDY merge 
    8    !!             3.5  !  2012    (S. Mocavero, I. Epicoco) Optimization of BDY communications 
    98   !!---------------------------------------------------------------------- 
    109#if defined   key_bdy   &&   defined key_lim2 
    … …  
    5453            CYCLE 
    5554         CASE(jp_frs) 
    56             CALL bdy_ice_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy ) 
     55            CALL bdy_ice_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy) ) 
    5756         CASE DEFAULT 
    5857            CALL ctl_stop( 'bdy_ice_lim_2 : unrecognised option for open boundaries for ice fields' ) 
    … …  
    6261   END SUBROUTINE bdy_ice_lim_2 
    6362 
    64    SUBROUTINE bdy_ice_frs( idx, dta, ib_bdy ) 
     63   SUBROUTINE bdy_ice_frs( idx, dta ) 
    6564      !!------------------------------------------------------------------------------ 
    6665      !!                 ***  SUBROUTINE bdy_ice_frs  *** 
    … …  
    7473      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices 
    7574      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
    76       INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index 
    7775      !! 
    78       INTEGER  ::   jb, jk, jgrd   ! dummy loop indices 
     76      INTEGER  ::   jb, jgrd   ! dummy loop indices 
    7977      INTEGER  ::   ii, ij         ! local scalar 
    8078      REAL(wp) ::   zwgt, zwgt1    ! local scalar 
    … …  
    8684      ! 
    8785      DO jb = 1, idx%nblen(jgrd) 
    88          DO jk = 1, jpkm1 
    8986            ii    = idx%nbi(jb,jgrd) 
    9087            ij    = idx%nbj(jb,jgrd) 
    9188            zwgt  = idx%nbw(jb,jgrd) 
    9289            zwgt1 = 1.e0 - idx%nbw(jb,jgrd) 
     90#if defined key_lim2_iceconc 
     91            frld (ii,ij) = ( frld (ii,ij) * zwgt1 + ( 1._wp - dta%frld (jb) ) * zwgt ) * tmask(ii,ij,1)     ! Leads fraction from ice fraction 
     92#else 
    9393            frld (ii,ij) = ( frld (ii,ij) * zwgt1 + dta%frld (jb) * zwgt ) * tmask(ii,ij,1)     ! Leads fraction  
     94#endif 
    9495            hicif(ii,ij) = ( hicif(ii,ij) * zwgt1 + dta%hicif(jb) * zwgt ) * tmask(ii,ij,1)     ! Ice depth  
    9596            hsnif(ii,ij) = ( hsnif(ii,ij) * zwgt1 + dta%hsnif(jb) * zwgt ) * tmask(ii,ij,1)     ! Snow depth 
    96          END DO 
    9797      END DO  
    98       CALL lbc_bdy_lnk( frld, 'T', 1., ib_bdy )                                         ! lateral boundary conditions 
    99       CALL lbc_bdy_lnk( hicif, 'T', 1., ib_bdy )   ;   CALL lbc_bdy_lnk( hsnif, 'T', 1., ib_bdy ) 
     98      CALL lbc_lnk( frld, 'T', 1. )                                         ! lateral boundary conditions 
     99      CALL lbc_lnk( hicif, 'T', 1. )   ;   CALL lbc_lnk( hsnif, 'T', 1. ) 
    100100      !       
    101101      IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_frs') 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r3703 r6736  
    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.4  !  2012     (J. Chanut) straight open boundary case update 
    14    !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Updates for the  
    15    !!                             optimization of BDY communications 
    1613   !!---------------------------------------------------------------------- 
    1714#if defined key_bdy 
    … …  
    2926   USE lib_mpp         ! for mpp_sum   
    3027   USE iom             ! I/O 
    31    USE sbctide, ONLY: lk_tide ! Tidal forcing or not 
    32    USE phycst, ONLY: rday 
    33  
    34    IMPLICIT NONE 
     28 
     29   IMPLICIT NONE  
    3530   PRIVATE 
    3631 
    3732   PUBLIC   bdy_init   ! routine called in nemo_init 
    3833 
    39    INTEGER, PARAMETER          :: jp_nseg = 100 
    40    INTEGER, PARAMETER          :: nrimmax = 20 ! maximum rimwidth in structured 
    41                                                ! open boundary data files 
    42    ! Straight open boundary segment parameters: 
    43    INTEGER  :: nbdysege, nbdysegw, nbdysegn, nbdysegs  
    44    INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft, npckge 
    45    INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft, npckgw 
    46    INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft, npckgn 
    47    INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft, npckgs 
    4834   !!---------------------------------------------------------------------- 
    4935   !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
    … …  
    6753      ! namelist variables 
    6854      !------------------- 
    69       CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile 
    70       CHARACTER(LEN=1)   ::   ctypebdy 
    71       INTEGER :: nbdyind, nbdybeg, nbdyend 
     55      INTEGER, PARAMETER          :: jp_nseg = 100 
     56      INTEGER                     :: nbdysege, nbdysegw, nbdysegn, nbdysegs  
     57      INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft 
     58      INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft 
     59      INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft 
     60      INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft 
    7261 
    7362      ! local variables 
    … …  
    7766      INTEGER  ::   iw, ie, is, in, inum, id_dummy         !   -       - 
    7867      INTEGER  ::   igrd_start, igrd_end, jpbdta           !   -       - 
    79       INTEGER  ::   jpbdtau, jpbdtas                       !   -       - 
    80       INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       - 
    8168      INTEGER, POINTER  ::  nbi, nbj, nbr                  ! short cuts 
    8269      REAL   , POINTER  ::  flagu, flagv                   !    -   - 
    8370      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars 
    84       INTEGER, DIMENSION (2)                  ::   kdimsz 
     71      INTEGER, DIMENSION (2)                ::   kdimsz 
    8572      INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays  
    8673      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta 
    8774      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points 
    88       CHARACTER(LEN=1),DIMENSION(jpbgrd)      ::   cgrid 
    89       INTEGER :: com_east, com_west, com_south, com_north          ! Flags for boundaries sending 
    90       INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b  ! Flags for boundaries receiving 
    91       INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4)                ! Arrays for neighbours coordinates 
    92  
     75      CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile 
     76      CHARACTER(LEN=1),DIMENSION(jpbgrd)   ::   cgrid 
    9377      !! 
    9478      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,             & 
    9579         &             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,              & 
     80         &             nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta, nb_jpk, &   
    9881#if defined key_lim2 
    9982         &             nn_ice_lim2, nn_ice_lim2_dta,                       & 
    … …  
    10184         &             ln_vol, nn_volctl, nn_rimwidth 
    10285      !! 
    103       NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend 
     86      NAMELIST/nambdy_index/ nbdysege, jpieob, jpjedt, jpjeft,             & 
     87                             nbdysegw, jpiwob, jpjwdt, jpjwft,             & 
     88                             nbdysegn, jpjnob, jpindt, jpinft,             & 
     89                             nbdysegs, jpjsob, jpisdt, jpisft 
    10490 
    10591      !!---------------------------------------------------------------------- 
    … …  
    118104 
    119105      cgrid= (/'t','u','v'/) 
    120        
     106 
    121107      ! ----------------------------------------- 
    122108      ! Initialise and read namelist parameters 
    … …  
    134120      nn_tra(:)         = 0 
    135121      nn_tra_dta(:)     = -1  ! uninitialised flag 
    136       ln_tra_dmp(:)     = .false. 
    137       ln_dyn3d_dmp(:)   = .false. 
    138       rn_time_dmp(:)    = 1. 
     122      nb_jpk            = -1 
    139123#if defined key_lim2 
    140124      nn_ice_lim2(:)    = 0 
    … …  
    151135      ! Check and write out namelist parameters 
    152136      ! ----------------------------------------- 
     137 
    153138      !                                   ! control prints 
    154139      IF(lwp) WRITE(numout,*) '         nambdy' 
    … …  
    173158        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  ' 
    174159        SELECT CASE( nn_dyn2d(ib_bdy) )                   
    175           CASE(jp_none)         ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    176           CASE(jp_frs)          ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    177           CASE(jp_flather)      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition' 
     160          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     161          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     162          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition' 
    178163          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' ) 
    179164        END SELECT 
    … …  
    186171              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' ) 
    187172           END SELECT 
    188            IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(.NOT.lk_tide)) THEN 
    189              CALL ctl_stop( 'You must activate key_tide to add tidal forcing at open boundaries' ) 
    190            ENDIF 
    191173        ENDIF 
    192174        IF(lwp) WRITE(numout,*) 
    … …  
    194176        IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  ' 
    195177        SELECT CASE( nn_dyn3d(ib_bdy) )                   
    196           CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    197           CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    198           CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value' 
    199           CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
     178          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     179          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    200180          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn3d' ) 
    201181        END SELECT 
    … …  
    207187           END SELECT 
    208188        ENDIF 
    209  
    210         IF ( ln_dyn3d_dmp(ib_bdy) ) THEN 
    211            IF ( nn_dyn3d(ib_bdy).EQ.0 ) THEN 
    212               IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.' 
    213               ln_dyn3d_dmp(ib_bdy)=.false. 
    214            ELSEIF ( nn_dyn3d(ib_bdy).EQ.1 ) THEN 
    215               CALL ctl_stop( 'Use FRS OR relaxation' ) 
    216            ELSE 
    217               IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone' 
    218               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
    219               IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 
    220            ENDIF 
    221         ELSE 
    222            IF(lwp) WRITE(numout,*) '      NO relaxation on baroclinic velocities' 
    223         ENDIF 
    224189        IF(lwp) WRITE(numout,*) 
    225190 
    226191        IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  ' 
    227192        SELECT CASE( nn_tra(ib_bdy) )                   
    228           CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    229           CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    230           CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value' 
    231           CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Neumann conditions' 
    232           CASE( 4 )      ;   IF(lwp) WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity' 
     193          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
     194          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
    233195          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' ) 
    234196        END SELECT 
    … …  
    239201              CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' ) 
    240202           END SELECT 
    241         ENDIF 
    242  
    243         IF ( ln_tra_dmp(ib_bdy) ) THEN 
    244            IF ( nn_tra(ib_bdy).EQ.0 ) THEN 
    245               IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.' 
    246               ln_tra_dmp(ib_bdy)=.false. 
    247            ELSEIF ( nn_tra(ib_bdy).EQ.1 ) THEN 
    248               CALL ctl_stop( 'Use FRS OR relaxation' ) 
    249            ELSE 
    250               IF(lwp) WRITE(numout,*) '      + T/S relaxation zone' 
    251               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
    252               IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' ) 
    253            ENDIF 
    254         ELSE 
    255            IF(lwp) WRITE(numout,*) '      NO T/S relaxation' 
    256203        ENDIF 
    257204        IF(lwp) WRITE(numout,*) 
    … …  
    274221#endif 
    275222 
    276         IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy) 
     223        IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS scheme = ', nn_rimwidth(ib_bdy) 
    277224        IF(lwp) WRITE(numout,*) 
    278225 
    279226      ENDDO 
    280227 
    281      IF (nb_bdy .gt. 0) THEN 
    282         IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value) 
    283           IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 
    284           IF(lwp) WRITE(numout,*) 
    285           SELECT CASE ( nn_volctl ) 
    286             CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant' 
    287             CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux' 
    288             CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 
    289           END SELECT 
    290           IF(lwp) WRITE(numout,*) 
    291         ELSE 
    292           IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 
    293           IF(lwp) WRITE(numout,*) 
    294         ENDIF 
     228     IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value) 
     229       IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 
     230       IF(lwp) WRITE(numout,*) 
     231       SELECT CASE ( nn_volctl ) 
     232         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant' 
     233         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux' 
     234         CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 
     235       END SELECT 
     236       IF(lwp) WRITE(numout,*) 
     237     ELSE 
     238       IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 
     239       IF(lwp) WRITE(numout,*) 
    295240     ENDIF 
    296241 
    … …  
    302247      ! --------------------------------------------- 
    303248      REWIND( numnam )                     
    304                 
    305       nblendta(:,:) = 0 
    306       nbdysege = 0 
    307       nbdysegw = 0 
    308       nbdysegn = 0 
    309       nbdysegs = 0 
    310       icount   = 0 ! count user defined segments 
    311       ! Dimensions below are used to allocate arrays to read external data 
    312       jpbdtas = 1 ! Maximum size of boundary data (structured case) 
    313       jpbdtau = 1 ! Maximum size of boundary data (unstructured case) 
    314  
    315249      DO ib_bdy = 1, nb_bdy 
    316250 
     251         jpbdta = 1 
    317252         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters 
    318253  
    319             icount = icount + 1 
    320254            ! No REWIND here because may need to read more than one nambdy_index namelist. 
    321255            READ  ( numnam, nambdy_index ) 
    322256 
    323             SELECT CASE ( TRIM(ctypebdy) ) 
    324               CASE( 'N' ) 
    325                  IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
    326                     nbdyind  = jpjglo - 2  ! set boundary to whole side of model domain. 
    327                     nbdybeg  = 2 
    328                     nbdyend  = jpiglo - 1 
    329                  ENDIF 
    330                  nbdysegn = nbdysegn + 1 
    331                  npckgn(nbdysegn) = ib_bdy ! Save bdy package number 
    332                  jpjnob(nbdysegn) = nbdyind 
    333                  jpindt(nbdysegn) = nbdybeg 
    334                  jpinft(nbdysegn) = nbdyend 
    335                  ! 
    336               CASE( 'S' ) 
    337                  IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
    338                     nbdyind  = 2           ! set boundary to whole side of model domain. 
    339                     nbdybeg  = 2 
    340                     nbdyend  = jpiglo - 1 
    341                  ENDIF 
    342                  nbdysegs = nbdysegs + 1 
    343                  npckgs(nbdysegs) = ib_bdy ! Save bdy package number 
    344                  jpjsob(nbdysegs) = nbdyind 
    345                  jpisdt(nbdysegs) = nbdybeg 
    346                  jpisft(nbdysegs) = nbdyend 
    347                  ! 
    348               CASE( 'E' ) 
    349                  IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
    350                     nbdyind  = jpiglo - 2  ! set boundary to whole side of model domain. 
    351                     nbdybeg  = 2 
    352                     nbdyend  = jpjglo - 1 
    353                  ENDIF 
    354                  nbdysege = nbdysege + 1  
    355                  npckge(nbdysege) = ib_bdy ! Save bdy package number 
    356                  jpieob(nbdysege) = nbdyind 
    357                  jpjedt(nbdysege) = nbdybeg 
    358                  jpjeft(nbdysege) = nbdyend 
    359                  ! 
    360               CASE( 'W' ) 
    361                  IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1 
    362                     nbdyind  = 2           ! set boundary to whole side of model domain. 
    363                     nbdybeg  = 2 
    364                     nbdyend  = jpjglo - 1 
    365                  ENDIF 
    366                  nbdysegw = nbdysegw + 1 
    367                  npckgw(nbdysegw) = ib_bdy ! Save bdy package number 
    368                  jpiwob(nbdysegw) = nbdyind 
    369                  jpjwdt(nbdysegw) = nbdybeg 
    370                  jpjwft(nbdysegw) = nbdyend 
    371                  ! 
    372               CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' ) 
    373             END SELECT 
    374  
    375             ! For simplicity we assume that in case of straight bdy, arrays have the same length 
    376             ! (even if it is true that last tangential velocity points 
    377             ! are useless). This simplifies a little bit boundary data format (and agrees with format 
    378             ! used so far in obc package) 
    379  
    380             nblendta(1:jpbgrd,ib_bdy) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(ib_bdy) 
    381             jpbdtas = MAX(jpbdtas, (nbdyend - nbdybeg + 1)) 
    382             IF (lwp.and.(nn_rimwidth(ib_bdy)>nrimmax)) & 
    383             & CALL ctl_stop( 'rimwidth must be lower than nrimmax' ) 
     257            ! Automatic boundary definition: if nbdysegX = -1 
     258            ! set boundary to whole side of model domain. 
     259            IF( nbdysege == -1 ) THEN 
     260               nbdysege = 1 
     261               jpieob(1) = jpiglo - 1 
     262               jpjedt(1) = 2 
     263               jpjeft(1) = jpjglo - 1 
     264            ENDIF 
     265            IF( nbdysegw == -1 ) THEN 
     266               nbdysegw = 1 
     267               jpiwob(1) = 2 
     268               jpjwdt(1) = 2 
     269               jpjwft(1) = jpjglo - 1 
     270            ENDIF 
     271            IF( nbdysegn == -1 ) THEN 
     272               nbdysegn = 1 
     273               jpjnob(1) = jpjglo - 1 
     274               jpindt(1) = 2 
     275               jpinft(1) = jpiglo - 1 
     276            ENDIF 
     277            IF( nbdysegs == -1 ) THEN 
     278               nbdysegs = 1 
     279               jpjsob(1) = 2 
     280               jpisdt(1) = 2 
     281               jpisft(1) = jpiglo - 1 
     282            ENDIF 
     283 
     284            nblendta(:,ib_bdy) = 0 
     285            DO iseg = 1, nbdysege 
     286               igrd = 1 
     287               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1                
     288               igrd = 2 
     289               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1                
     290               igrd = 3 
     291               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg)                
     292            ENDDO 
     293            DO iseg = 1, nbdysegw 
     294               igrd = 1 
     295               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1                
     296               igrd = 2 
     297               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1                
     298               igrd = 3 
     299               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg)                
     300            ENDDO 
     301            DO iseg = 1, nbdysegn 
     302               igrd = 1 
     303               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1                
     304               igrd = 2 
     305               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg)                
     306               igrd = 3 
     307               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1 
     308            ENDDO 
     309            DO iseg = 1, nbdysegs 
     310               igrd = 1 
     311               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1                
     312               igrd = 2 
     313               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) 
     314               igrd = 3 
     315               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1                
     316            ENDDO 
     317 
     318            nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy) 
     319            jpbdta = MAXVAL(nblendta(:,ib_bdy))                
     320 
    384321 
    385322         ELSE            ! Read size of arrays in boundary coordinates file. 
     323 
     324 
    386325            CALL iom_open( cn_coords_file(ib_bdy), inum ) 
     326            jpbdta = 1 
    387327            DO igrd = 1, jpbgrd 
    388328               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )   
    389329               nblendta(igrd,ib_bdy) = kdimsz(1) 
    390                jpbdtau = MAX(jpbdtau, kdimsz(1)) 
    391             ENDDO 
    392             CALL iom_close( inum ) 
     330               jpbdta = MAX(jpbdta, kdimsz(1)) 
     331            ENDDO 
    393332 
    394333         ENDIF  
    … …  
    396335      ENDDO ! ib_bdy 
    397336 
    398       IF (nb_bdy>0) THEN 
    399          jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy)) 
    400  
    401          ! Allocate arrays 
    402          !--------------- 
    403          ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    & 
    404             &      nbrdta(jpbdta, jpbgrd, nb_bdy) ) 
    405  
    406          ALLOCATE( dta_global(jpbdtau, 1, jpk) ) 
    407          IF ( icount>0 ) ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) ) 
    408          !  
    409       ENDIF 
    410  
    411       ! Now look for crossings in user (namelist) defined open boundary segments: 
    412       !-------------------------------------------------------------------------- 
    413       IF ( icount>0 ) CALL bdy_ctl_seg 
     337      ! Allocate arrays 
     338      !--------------- 
     339      ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    & 
     340         &      nbrdta(jpbdta, jpbgrd, nb_bdy) ) 
     341 
     342      ALLOCATE( dta_global(jpbdta, 1, jpk) ) 
     343      ALLOCATE( dta_global_1(jpbdta, 1, jpk) ) 
     344      ALLOCATE( dta_global_2(jpbdta, jpk) ) 
    414345 
    415346      ! Calculate global boundary index arrays or read in from file 
    416       !------------------------------------------------------------                
    417       ! 1. Read global index arrays from boundary coordinates file. 
     347      !------------------------------------------------------------ 
     348      REWIND( numnam )                     
    418349      DO ib_bdy = 1, nb_bdy 
    419350 
    420          IF( ln_coords_file(ib_bdy) ) THEN 
    421  
    422             CALL iom_open( cn_coords_file(ib_bdy), inum ) 
     351         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Calculate global index arrays from namelist parameters 
     352 
     353            ! No REWIND here because may need to read more than one nambdy_index namelist. 
     354            READ  ( numnam, nambdy_index ) 
     355 
     356            ! Automatic boundary definition: if nbdysegX = -1 
     357            ! set boundary to whole side of model domain. 
     358            IF( nbdysege == -1 ) THEN 
     359               nbdysege = 1 
     360               jpieob(1) = jpiglo - 1 
     361               jpjedt(1) = 2 
     362               jpjeft(1) = jpjglo - 1 
     363            ENDIF 
     364            IF( nbdysegw == -1 ) THEN 
     365               nbdysegw = 1 
     366               jpiwob(1) = 2 
     367               jpjwdt(1) = 2 
     368               jpjwft(1) = jpjglo - 1 
     369            ENDIF 
     370            IF( nbdysegn == -1 ) THEN 
     371               nbdysegn = 1 
     372               jpjnob(1) = jpjglo - 1 
     373               jpindt(1) = 2 
     374               jpinft(1) = jpiglo - 1 
     375            ENDIF 
     376            IF( nbdysegs == -1 ) THEN 
     377               nbdysegs = 1 
     378               jpjsob(1) = 2 
     379               jpisdt(1) = 2 
     380               jpisft(1) = jpiglo - 1 
     381            ENDIF 
     382 
     383            ! ------------ T points ------------- 
     384            igrd = 1   
     385            icount = 0 
     386            DO ir = 1, nn_rimwidth(ib_bdy) 
     387               ! east 
     388               DO iseg = 1, nbdysege 
     389                  DO ij = jpjedt(iseg), jpjeft(iseg) 
     390                     icount = icount + 1 
     391                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1 
     392                     nbjdta(icount, igrd, ib_bdy) = ij 
     393                     nbrdta(icount, igrd, ib_bdy) = ir 
     394                  ENDDO 
     395               ENDDO 
     396               ! west 
     397               DO iseg = 1, nbdysegw 
     398                  DO ij = jpjwdt(iseg), jpjwft(iseg) 
     399                     icount = icount + 1 
     400                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     401                     nbjdta(icount, igrd, ib_bdy) = ij 
     402                     nbrdta(icount, igrd, ib_bdy) = ir 
     403                  ENDDO 
     404               ENDDO 
     405               ! north 
     406               DO iseg = 1, nbdysegn 
     407                  DO ii = jpindt(iseg), jpinft(iseg) 
     408                     icount = icount + 1 
     409                     nbidta(icount, igrd, ib_bdy) = ii 
     410                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1 
     411                     nbrdta(icount, igrd, ib_bdy) = ir 
     412                  ENDDO 
     413               ENDDO 
     414               ! south 
     415               DO iseg = 1, nbdysegs 
     416                  DO ii = jpisdt(iseg), jpisft(iseg) 
     417                     icount = icount + 1 
     418                     nbidta(icount, igrd, ib_bdy) = ii 
     419                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
     420                     nbrdta(icount, igrd, ib_bdy) = ir 
     421                  ENDDO 
     422               ENDDO 
     423            ENDDO 
     424 
     425            ! ------------ U points ------------- 
     426            igrd = 2   
     427            icount = 0 
     428            DO ir = 1, nn_rimwidth(ib_bdy) 
     429               ! east 
     430               DO iseg = 1, nbdysege 
     431                  DO ij = jpjedt(iseg), jpjeft(iseg) 
     432                     icount = icount + 1 
     433                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir 
     434                     nbjdta(icount, igrd, ib_bdy) = ij 
     435                     nbrdta(icount, igrd, ib_bdy) = ir 
     436                  ENDDO 
     437               ENDDO 
     438               ! west 
     439               DO iseg = 1, nbdysegw 
     440                  DO ij = jpjwdt(iseg), jpjwft(iseg) 
     441                     icount = icount + 1 
     442                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     443                     nbjdta(icount, igrd, ib_bdy) = ij 
     444                     nbrdta(icount, igrd, ib_bdy) = ir 
     445                  ENDDO 
     446               ENDDO 
     447               ! north 
     448               DO iseg = 1, nbdysegn 
     449                  DO ii = jpindt(iseg), jpinft(iseg) - 1 
     450                     icount = icount + 1 
     451                     nbidta(icount, igrd, ib_bdy) = ii 
     452                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1 
     453                     nbrdta(icount, igrd, ib_bdy) = ir 
     454                  ENDDO 
     455               ENDDO 
     456               ! south 
     457               DO iseg = 1, nbdysegs 
     458                  DO ii = jpisdt(iseg), jpisft(iseg) - 1 
     459                     icount = icount + 1 
     460                     nbidta(icount, igrd, ib_bdy) = ii 
     461                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
     462                     nbrdta(icount, igrd, ib_bdy) = ir 
     463                  ENDDO 
     464               ENDDO 
     465            ENDDO 
     466 
     467            ! ------------ V points ------------- 
     468            igrd = 3   
     469            icount = 0 
     470            DO ir = 1, nn_rimwidth(ib_bdy) 
     471               ! east 
     472               DO iseg = 1, nbdysege 
     473                  DO ij = jpjedt(iseg), jpjeft(iseg) - 1 
     474                     icount = icount + 1 
     475                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1 
     476                     nbjdta(icount, igrd, ib_bdy) = ij 
     477                     nbrdta(icount, igrd, ib_bdy) = ir 
     478                  ENDDO 
     479               ENDDO 
     480               ! west 
     481               DO iseg = 1, nbdysegw 
     482                  DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 
     483                     icount = icount + 1 
     484                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
     485                     nbjdta(icount, igrd, ib_bdy) = ij 
     486                     nbrdta(icount, igrd, ib_bdy) = ir 
     487                  ENDDO 
     488               ENDDO 
     489               ! north 
     490               DO iseg = 1, nbdysegn 
     491                  DO ii = jpindt(iseg), jpinft(iseg) 
     492                     icount = icount + 1 
     493                     nbidta(icount, igrd, ib_bdy) = ii 
     494                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir 
     495                     nbrdta(icount, igrd, ib_bdy) = ir 
     496                  ENDDO 
     497               ENDDO 
     498               ! south 
     499               DO iseg = 1, nbdysegs 
     500                  DO ii = jpisdt(iseg), jpisft(iseg) 
     501                     icount = icount + 1 
     502                     nbidta(icount, igrd, ib_bdy) = ii 
     503                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
     504                     nbrdta(icount, igrd, ib_bdy) = ir 
     505                  ENDDO 
     506               ENDDO 
     507            ENDDO 
     508 
     509         ELSE            ! Read global index arrays from boundary coordinates file. 
     510 
    423511            DO igrd = 1, jpbgrd 
    424512               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) 
    … …  
    441529               IF (ibr_max < nn_rimwidth(ib_bdy))   & 
    442530                     CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) ) 
     531 
    443532            END DO 
    444533            CALL iom_close( inum ) 
    … …  
    446535         ENDIF  
    447536 
    448       ENDDO       
    449      
    450       ! 2. Now fill indices corresponding to straight open boundary arrays: 
    451       ! East 
    452       !----- 
    453       DO iseg = 1, nbdysege 
    454          ib_bdy = npckge(iseg) 
    455          ! 
    456          ! ------------ T points ------------- 
    457          igrd=1 
    458          icount=0 
    459          DO ir = 1, nn_rimwidth(ib_bdy) 
    460             DO ij = jpjedt(iseg), jpjeft(iseg) 
    461                icount = icount + 1 
    462                nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 
    463                nbjdta(icount, igrd, ib_bdy) = ij 
    464                nbrdta(icount, igrd, ib_bdy) = ir 
    465             ENDDO 
    466          ENDDO 
    467          ! 
    468          ! ------------ U points ------------- 
    469          igrd=2 
    470          icount=0 
    471          DO ir = 1, nn_rimwidth(ib_bdy) 
    472             DO ij = jpjedt(iseg), jpjeft(iseg) 
    473                icount = icount + 1 
    474                nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir 
    475                nbjdta(icount, igrd, ib_bdy) = ij 
    476                nbrdta(icount, igrd, ib_bdy) = ir 
    477             ENDDO 
    478          ENDDO 
    479          ! 
    480          ! ------------ V points ------------- 
    481          igrd=3 
    482          icount=0 
    483          DO ir = 1, nn_rimwidth(ib_bdy) 
    484 !            DO ij = jpjedt(iseg), jpjeft(iseg) - 1 
    485             DO ij = jpjedt(iseg), jpjeft(iseg) 
    486                icount = icount + 1 
    487                nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir 
    488                nbjdta(icount, igrd, ib_bdy) = ij 
    489                nbrdta(icount, igrd, ib_bdy) = ir 
    490             ENDDO 
    491             nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    492             nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    493          ENDDO 
    494       ENDDO 
    495       ! 
    496       ! West 
    497       !----- 
    498       DO iseg = 1, nbdysegw 
    499          ib_bdy = npckgw(iseg) 
    500          ! 
    501          ! ------------ T points ------------- 
    502          igrd=1 
    503          icount=0 
    504          DO ir = 1, nn_rimwidth(ib_bdy) 
    505             DO ij = jpjwdt(iseg), jpjwft(iseg) 
    506                icount = icount + 1 
    507                nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    508                nbjdta(icount, igrd, ib_bdy) = ij 
    509                nbrdta(icount, igrd, ib_bdy) = ir 
    510             ENDDO 
    511          ENDDO 
    512          ! 
    513          ! ------------ U points ------------- 
    514          igrd=2 
    515          icount=0 
    516          DO ir = 1, nn_rimwidth(ib_bdy) 
    517             DO ij = jpjwdt(iseg), jpjwft(iseg) 
    518                icount = icount + 1 
    519                nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    520                nbjdta(icount, igrd, ib_bdy) = ij 
    521                nbrdta(icount, igrd, ib_bdy) = ir 
    522             ENDDO 
    523          ENDDO 
    524          ! 
    525          ! ------------ V points ------------- 
    526          igrd=3 
    527          icount=0 
    528          DO ir = 1, nn_rimwidth(ib_bdy) 
    529 !            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1 
    530             DO ij = jpjwdt(iseg), jpjwft(iseg) 
    531                icount = icount + 1 
    532                nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1 
    533                nbjdta(icount, igrd, ib_bdy) = ij 
    534                nbrdta(icount, igrd, ib_bdy) = ir 
    535             ENDDO 
    536             nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    537             nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    538          ENDDO 
    539       ENDDO 
    540       ! 
    541       ! North 
    542       !----- 
    543       DO iseg = 1, nbdysegn 
    544          ib_bdy = npckgn(iseg) 
    545          ! 
    546          ! ------------ T points ------------- 
    547          igrd=1 
    548          icount=0 
    549          DO ir = 1, nn_rimwidth(ib_bdy) 
    550             DO ii = jpindt(iseg), jpinft(iseg) 
    551                icount = icount + 1 
    552                nbidta(icount, igrd, ib_bdy) = ii 
    553                nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir  
    554                nbrdta(icount, igrd, ib_bdy) = ir 
    555             ENDDO 
    556          ENDDO 
    557          ! 
    558          ! ------------ U points ------------- 
    559          igrd=2 
    560          icount=0 
    561          DO ir = 1, nn_rimwidth(ib_bdy) 
    562 !            DO ii = jpindt(iseg), jpinft(iseg) - 1 
    563             DO ii = jpindt(iseg), jpinft(iseg) 
    564                icount = icount + 1 
    565                nbidta(icount, igrd, ib_bdy) = ii 
    566                nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 
    567                nbrdta(icount, igrd, ib_bdy) = ir 
    568             ENDDO 
    569             nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    570             nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    571          ENDDO 
    572          ! 
    573          ! ------------ V points ------------- 
    574          igrd=3 
    575          icount=0 
    576          DO ir = 1, nn_rimwidth(ib_bdy) 
    577             DO ii = jpindt(iseg), jpinft(iseg) 
    578                icount = icount + 1 
    579                nbidta(icount, igrd, ib_bdy) = ii 
    580                nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir 
    581                nbrdta(icount, igrd, ib_bdy) = ir 
    582             ENDDO 
    583          ENDDO 
    584       ENDDO 
    585       ! 
    586       ! South 
    587       !----- 
    588       DO iseg = 1, nbdysegs 
    589          ib_bdy = npckgs(iseg) 
    590          ! 
    591          ! ------------ T points ------------- 
    592          igrd=1 
    593          icount=0 
    594          DO ir = 1, nn_rimwidth(ib_bdy) 
    595             DO ii = jpisdt(iseg), jpisft(iseg) 
    596                icount = icount + 1 
    597                nbidta(icount, igrd, ib_bdy) = ii 
    598                nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    599                nbrdta(icount, igrd, ib_bdy) = ir 
    600             ENDDO 
    601          ENDDO 
    602          ! 
    603          ! ------------ U points ------------- 
    604          igrd=2 
    605          icount=0 
    606          DO ir = 1, nn_rimwidth(ib_bdy) 
    607 !            DO ii = jpisdt(iseg), jpisft(iseg) - 1 
    608             DO ii = jpisdt(iseg), jpisft(iseg) 
    609                icount = icount + 1 
    610                nbidta(icount, igrd, ib_bdy) = ii 
    611                nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    612                nbrdta(icount, igrd, ib_bdy) = ir 
    613             ENDDO 
    614             nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    615             nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point 
    616          ENDDO 
    617          ! 
    618          ! ------------ V points ------------- 
    619          igrd=3 
    620          icount=0 
    621          DO ir = 1, nn_rimwidth(ib_bdy) 
    622             DO ii = jpisdt(iseg), jpisft(iseg) 
    623                icount = icount + 1 
    624                nbidta(icount, igrd, ib_bdy) = ii 
    625                nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1 
    626                nbrdta(icount, igrd, ib_bdy) = ir 
    627             ENDDO 
    628          ENDDO 
    629       ENDDO 
    630  
    631       !  Deal with duplicated points 
    632       !----------------------------- 
    633       ! We assign negative indices to duplicated points (to remove them from bdy points to be updated) 
    634       ! if their distance to the bdy is greater than the other 
    635       ! If their distance are the same, just keep only one to avoid updating a point twice 
    636       DO igrd = 1, jpbgrd 
    637          DO ib_bdy1 = 1, nb_bdy 
    638             DO ib_bdy2 = 1, nb_bdy 
    639                IF (ib_bdy1/=ib_bdy2) THEN 
    640                   DO ib1 = 1, nblendta(igrd,ib_bdy1) 
    641                      DO ib2 = 1, nblendta(igrd,ib_bdy2) 
    642                         IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. & 
    643                         &   (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN 
    644 !                           IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &  
    645 !                                                       &              nbidta(ib1, igrd, ib_bdy1),      &  
    646 !                                                       &              nbjdta(ib2, igrd, ib_bdy2) 
    647                            ! keep only points with the lowest distance to boundary: 
    648                            IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN 
    649                              nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2 
    650                              nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2 
    651                            ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN 
    652                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 
    653                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 
    654                            ! Arbitrary choice if distances are the same: 
    655                            ELSE 
    656                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1 
    657                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1 
    658                            ENDIF 
    659                         END IF 
    660                      END DO 
    661                   END DO 
    662                ENDIF 
    663             END DO 
    664          END DO 
    665       END DO 
     537      ENDDO  
    666538 
    667539      ! Work out dimensions of boundary data on each processor 
    668540      ! ------------------------------------------------------ 
    669  
    670       ! Rather assume that boundary data indices are given on global domain 
    671       ! TO BE DISCUSSED ? 
    672 !      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2 
    673 !      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1 
    674 !      is = mjg(1) + 1            ! if monotasking and no zoom, is=2 
    675 !      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1       
    676       iw = mig(1) - jpizoom + 2         ! if monotasking and no zoom, iw=2 
    677       ie = mig(1) + nlci - jpizoom - 1  ! if monotasking and no zoom, ie=jpim1 
    678       is = mjg(1) - jpjzoom + 2         ! if monotasking and no zoom, is=2 
    679       in = mjg(1) + nlcj - jpjzoom - 1  ! if monotasking and no zoom, in=jpjm1 
    680  
    681       ALLOCATE( nbondi_bdy(nb_bdy)) 
    682       ALLOCATE( nbondj_bdy(nb_bdy)) 
    683       nbondi_bdy(:)=2 
    684       nbondj_bdy(:)=2 
    685       ALLOCATE( nbondi_bdy_b(nb_bdy)) 
    686       ALLOCATE( nbondj_bdy_b(nb_bdy)) 
    687       nbondi_bdy_b(:)=2 
    688       nbondj_bdy_b(:)=2 
    689  
    690       ! Work out dimensions of boundary data on each neighbour process 
    691       IF(nbondi .eq. 0) THEN 
    692          iw_b(1) = jpizoom + nimppt(nowe+1) 
    693          ie_b(1) = jpizoom + nimppt(nowe+1)+nlcit(nowe+1)-3 
    694          is_b(1) = jpjzoom + njmppt(nowe+1) 
    695          in_b(1) = jpjzoom + njmppt(nowe+1)+nlcjt(nowe+1)-3 
    696  
    697          iw_b(2) = jpizoom + nimppt(noea+1) 
    698          ie_b(2) = jpizoom + nimppt(noea+1)+nlcit(noea+1)-3 
    699          is_b(2) = jpjzoom + njmppt(noea+1) 
    700          in_b(2) = jpjzoom + njmppt(noea+1)+nlcjt(noea+1)-3 
    701       ELSEIF(nbondi .eq. 1) THEN 
    702          iw_b(1) = jpizoom + nimppt(nowe+1) 
    703          ie_b(1) = jpizoom + nimppt(nowe+1)+nlcit(nowe+1)-3 
    704          is_b(1) = jpjzoom + njmppt(nowe+1) 
    705          in_b(1) = jpjzoom + njmppt(nowe+1)+nlcjt(nowe+1)-3 
    706       ELSEIF(nbondi .eq. -1) THEN 
    707          iw_b(2) = jpizoom + nimppt(noea+1) 
    708          ie_b(2) = jpizoom + nimppt(noea+1)+nlcit(noea+1)-3 
    709          is_b(2) = jpjzoom + njmppt(noea+1) 
    710          in_b(2) = jpjzoom + njmppt(noea+1)+nlcjt(noea+1)-3 
    711       ENDIF 
    712  
    713       IF(nbondj .eq. 0) THEN 
    714          iw_b(3) = jpizoom + nimppt(noso+1) 
    715          ie_b(3) = jpizoom + nimppt(noso+1)+nlcit(noso+1)-3 
    716          is_b(3) = jpjzoom + njmppt(noso+1) 
    717          in_b(3) = jpjzoom + njmppt(noso+1)+nlcjt(noso+1)-3 
    718  
    719          iw_b(4) = jpizoom + nimppt(nono+1) 
    720          ie_b(4) = jpizoom + nimppt(nono+1)+nlcit(nono+1)-3 
    721          is_b(4) = jpjzoom + njmppt(nono+1) 
    722          in_b(4) = jpjzoom + njmppt(nono+1)+nlcjt(nono+1)-3 
    723       ELSEIF(nbondj .eq. 1) THEN 
    724          iw_b(3) = jpizoom + nimppt(noso+1) 
    725          ie_b(3) = jpizoom + nimppt(noso+1)+nlcit(noso+1)-3 
    726          is_b(3) = jpjzoom + njmppt(noso+1) 
    727          in_b(3) = jpjzoom + njmppt(noso+1)+nlcjt(noso+1)-3 
    728       ELSEIF(nbondj .eq. -1) THEN 
    729          iw_b(4) = jpizoom + nimppt(nono+1) 
    730          ie_b(4) = jpizoom + nimppt(nono+1)+nlcit(nono+1)-3 
    731          is_b(4) = jpjzoom + njmppt(nono+1) 
    732          in_b(4) = jpjzoom + njmppt(nono+1)+nlcjt(nono+1)-3 
    733       ENDIF 
     541      
     542      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2 
     543      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1 
     544      is = mjg(1) + 1            ! if monotasking and no zoom, is=2 
     545      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1 
    734546 
    735547      DO ib_bdy = 1, nb_bdy 
    … …  
    744556               IF(lwp) THEN         ! Since all procs read global data only need to do this check on one proc... 
    745557                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN 
    746                      CALL ctl_stop('bdy_init : ERROR : boundary data in file  & 
    747                                     must be defined in order of distance from edge nbr.', & 
    748                                    'A utility for re-ordering boundary coordinates and data & 
    749                                     files exists in the TOOLS/OBC directory') 
     558                     CALL ctl_stop('bdy_init : ERROR : boundary data in file must be defined in order of distance from edge nbr.', & 
     559                    'A utility for re-ordering boundary coordinates and data files exists in the TOOLS/OBC directory') 
    750560                  ENDIF     
    751561               ENDIF 
    … …  
    769579         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) ) 
    770580         ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) ) 
    771          ALLOCATE( idx_bdy(ib_bdy)%nbd(ilen1,jpbgrd) ) 
     581         ALLOCATE( idx_bdy(ib_bdy)%nbz(ilen1,jpbgrd,jpk) )     ! jdha addition TODO use this instead of calculating in fldread?  
    772582         ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) ) 
    773583         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) ) 
    … …  
    778588         ! ----------------------------------------------------------------- 
    779589 
    780          com_east = 0 
    781          com_west = 0 
    782          com_south = 0 
    783          com_north = 0 
    784  
    785          com_east_b = 0 
    786          com_west_b = 0 
    787          com_south_b = 0 
    788          com_north_b = 0 
    789590         DO igrd = 1, jpbgrd 
    790591            icount  = 0 
    … …  
    798599                     ! 
    799600                     icount = icount  + 1 
    800  
    801                      ! Rather assume that boundary data indices are given on global domain 
    802                      ! TO BE DISCUSSED ? 
    803 !                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1 
    804 !                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 
    805                      idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+jpizoom 
    806                      idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+jpjzoom 
    807                      ! check if point has to be sent 
    808                      ii = idx_bdy(ib_bdy)%nbi(icount,igrd) 
    809                      ij = idx_bdy(ib_bdy)%nbj(icount,igrd) 
    810                      if((com_east .ne. 1) .and. (ii .eq. (nlci-1)) .and. (nbondi .le. 0)) then 
    811                         com_east = 1 
    812                      elseif((com_west .ne. 1) .and. (ii .eq. 2) .and. (nbondi .ge. 0) .and. (nbondi .ne. 2)) then 
    813                         com_west = 1 
    814                      endif  
    815                      if((com_south .ne. 1) .and. (ij .eq. 2) .and. (nbondj .ge. 0) .and. (nbondj .ne. 2)) then 
    816                         com_south = 1 
    817                      elseif((com_north .ne. 1) .and. (ij .eq. (nlcj-1)) .and. (nbondj .le. 0)) then 
    818                         com_north = 1 
    819                      endif  
     601                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1 
     602                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1 
    820603                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy) 
     604                     DO ik = 1,jpk 
     605                        idx_bdy(ib_bdy)%nbz(icount,igrd,ik) =                                  & 
     606                     &  gdept_1(idx_bdy(ib_bdy)%nbi(icount,igrd),idx_bdy(ib_bdy)%nbj(icount,igrd),ik)   ! if using in step could use fsdept? 
     607                     ENDDO 
    821608                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib 
    822609                  ENDIF 
    823                   ! check if point has to be received from a neighbour 
    824                   IF(nbondi .eq. 0) THEN 
    825                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND.   & 
    826                        & nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND.   & 
    827                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    828                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2 
    829                        if((com_west_b .ne. 1) .and. (ii .eq. (nlcit(nowe+1)-1))) then 
    830                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2 
    831                           if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then 
    832                             com_south = 1 
    833                           elseif((ij .eq. nlcjt(nowe+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then 
    834                             com_north = 1 
    835                           endif 
    836                           com_west_b = 1 
    837                        endif  
    838                      ENDIF 
    839                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND.   & 
    840                        & nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND.   & 
    841                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    842                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2 
    843                        if((com_east_b .ne. 1) .and. (ii .eq. 2)) then 
    844                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2 
    845                           if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then 
    846                             com_south = 1 
    847                           elseif((ij .eq. nlcjt(noea+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then 
    848                             com_north = 1 
    849                           endif 
    850                           com_east_b = 1 
    851                        endif  
    852                      ENDIF 
    853                   ELSEIF(nbondi .eq. 1) THEN 
    854                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND.   & 
    855                        & nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND.   & 
    856                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    857                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2 
    858                        if((com_west_b .ne. 1) .and. (ii .eq. (nlcit(nowe+1)-1))) then 
    859                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2 
    860                           if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then 
    861                             com_south = 1 
    862                           elseif((ij .eq. nlcjt(nowe+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then 
    863                             com_north = 1 
    864                           endif 
    865                           com_west_b = 1 
    866                        endif  
    867                      ENDIF 
    868                   ELSEIF(nbondi .eq. -1) THEN 
    869                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND.   & 
    870                        & nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND.   & 
    871                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    872                        ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2 
    873                        if((com_east_b .ne. 1) .and. (ii .eq. 2)) then 
    874                           ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2 
    875                           if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then 
    876                             com_south = 1 
    877                           elseif((ij .eq. nlcjt(noea+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then 
    878                             com_north = 1 
    879                           endif 
    880                           com_east_b = 1 
    881                        endif  
    882                      ENDIF 
    883                   ENDIF 
    884                   IF(nbondj .eq. 0) THEN 
    885                      IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1  & 
    886                        & .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. & 
    887                        & nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN 
    888                        com_north_b = 1  
    889                      ENDIF 
    890                      IF(com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1  & 
    891                        &.OR. nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. & 
    892                        & nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN 
    893                        com_south_b = 1  
    894                      ENDIF 
    895                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND.   & 
    896                        & nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND.   & 
    897                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    898                        ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2 
    899                        if((com_south_b .ne. 1) .and. (ij .eq. (nlcjt(noso+1)-1))) then 
    900                           com_south_b = 1 
    901                        endif  
    902                      ENDIF 
    903                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND.   & 
    904                        & nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND.   & 
    905                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    906                        ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2 
    907                        if((com_north_b .ne. 1) .and. (ij .eq. 2)) then 
    908                           com_north_b = 1 
    909                        endif  
    910                      ENDIF 
    911                   ELSEIF(nbondj .eq. 1) THEN 
    912                      IF( com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 .OR. & 
    913                        & nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. & 
    914                        & nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN 
    915                        com_south_b = 1  
    916                      ENDIF 
    917                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND.   & 
    918                        & nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND.   & 
    919                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    920                        ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2 
    921                        if((com_south_b .ne. 1) .and. (ij .eq. (nlcjt(noso+1)-1))) then 
    922                           com_south_b = 1 
    923                        endif  
    924                      ENDIF 
    925                   ELSEIF(nbondj .eq. -1) THEN 
    926                      IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1  & 
    927                        & .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. & 
    928                        & nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN 
    929                        com_north_b = 1  
    930                      ENDIF 
    931                      IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND.   & 
    932                        & nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND.   & 
    933                        & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN 
    934                        ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2 
    935                        if((com_north_b .ne. 1) .and. (ij .eq. 2)) then 
    936                           com_north_b = 1 
    937                        endif  
    938                      ENDIF 
    939                   ENDIF 
    940610               ENDDO 
    941611            ENDDO 
    942612         ENDDO  
    943          ! definition of the i- and j- direction local boundaries arrays 
    944          ! used for sending the boudaries 
    945          IF((com_east .eq. 1) .and. (com_west .eq. 1)) THEN 
    946             nbondi_bdy(ib_bdy) = 0 
    947          ELSEIF ((com_east .eq. 1) .and. (com_west .eq. 0)) THEN 
    948             nbondi_bdy(ib_bdy) = -1 
    949          ELSEIF ((com_east .eq. 0) .and. (com_west .eq. 1)) THEN 
    950             nbondi_bdy(ib_bdy) = 1 
    951          ENDIF 
    952  
    953          IF((com_north .eq. 1) .and. (com_south .eq. 1)) THEN 
    954             nbondj_bdy(ib_bdy) = 0 
    955          ELSEIF ((com_north .eq. 1) .and. (com_south .eq. 0)) THEN 
    956             nbondj_bdy(ib_bdy) = -1 
    957          ELSEIF ((com_north .eq. 0) .and. (com_south .eq. 1)) THEN 
    958             nbondj_bdy(ib_bdy) = 1 
    959          ENDIF 
    960  
    961          ! definition of the i- and j- direction local boundaries arrays 
    962          ! used for receiving the boudaries 
    963          IF((com_east_b .eq. 1) .and. (com_west_b .eq. 1)) THEN 
    964             nbondi_bdy_b(ib_bdy) = 0 
    965          ELSEIF ((com_east_b .eq. 1) .and. (com_west_b .eq. 0)) THEN 
    966             nbondi_bdy_b(ib_bdy) = -1 
    967          ELSEIF ((com_east_b .eq. 0) .and. (com_west_b .eq. 1)) THEN 
    968             nbondi_bdy_b(ib_bdy) = 1 
    969          ENDIF 
    970  
    971          IF((com_north_b .eq. 1) .and. (com_south_b .eq. 1)) THEN 
    972             nbondj_bdy_b(ib_bdy) = 0 
    973          ELSEIF ((com_north_b .eq. 1) .and. (com_south_b .eq. 0)) THEN 
    974             nbondj_bdy_b(ib_bdy) = -1 
    975          ELSEIF ((com_north_b .eq. 0) .and. (com_south_b .eq. 1)) THEN 
    976             nbondj_bdy_b(ib_bdy) = 1 
    977          ENDIF 
    978613 
    979614         ! Compute rim weights for FRS scheme 
    … …  
    983618               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 
    984619               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation 
    985 !               idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.  ! quadratic 
    986 !               idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy))       ! linear 
    987             END DO 
    988          END DO  
    989  
    990          ! Compute damping coefficients 
    991          ! ---------------------------- 
    992          DO igrd = 1, jpbgrd 
    993             DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 
    994                nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 
    995                idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) &  
    996                & *(FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.   ! quadratic 
     620!              idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2      ! quadratic 
     621!              idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth)          ! linear 
    997622            END DO 
    998623         END DO  
    … …  
    1014639         CALL iom_close( inum ) 
    1015640 
     641               IF(lwp) WRITE(numout,*) 'get bdytmask', bdytmask   
    1016642         ! Derive mask on U and V grid from mask on T grid 
    1017643         bdyumask(:,:) = 0.e0 
    … …  
    1053679       
    1054680      bdytmask(:,:) = tmask(:,:,1) 
     681      IF( .not. ln_mask_file ) THEN 
     682         ! If .not. ln_mask_file then we need to derive mask on U and V grid  
     683         ! from mask on T grid here. 
     684         bdyumask(:,:) = 0.e0 
     685         bdyvmask(:,:) = 0.e0 
     686         DO ij=1, jpjm1 
     687            DO ii=1, jpim1 
     688               bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij ) 
     689               bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1)   
     690            END DO 
     691         END DO 
     692         CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond. 
     693      ENDIF 
    1055694 
    1056695      ! bdy masks and bmask are now set to zero on boundary points: 
    … …  
    1126765            END IF 
    1127766         END DO 
    1128  
     767  
    1129768         IF( icount /= 0 ) THEN 
    1130769            IF(lwp) WRITE(numout,*) 
    … …  
    1140779      ! Compute total lateral surface for volume correction: 
    1141780      ! ---------------------------------------------------- 
    1142       ! JC: this must be done at each time step with key_vvl 
    1143781      bdysurftot = 0.e0  
    1144782      IF( ln_vol ) THEN   
    … …  
    1174812      ! Tidy up 
    1175813      !-------- 
    1176       IF (nb_bdy>0) THEN 
    1177          DEALLOCATE(nbidta, nbjdta, nbrdta) 
    1178       ENDIF 
     814      DEALLOCATE(nbidta, nbjdta, nbrdta) 
    1179815 
    1180816      IF( nn_timing == 1 ) CALL timing_stop('bdy_init') 
    1181817 
    1182818   END SUBROUTINE bdy_init 
    1183  
    1184    SUBROUTINE bdy_ctl_seg 
    1185       !!---------------------------------------------------------------------- 
    1186       !!                 ***  ROUTINE bdy_ctl_seg  *** 
    1187       !! 
    1188       !! ** Purpose :   Check straight open boundary segments location 
    1189       !! 
    1190       !! ** Method  :   - Look for open boundary corners 
    1191       !!                - Check that segments start or end on land  
    1192       !!---------------------------------------------------------------------- 
    1193       INTEGER  ::   ib, ib1, ib2, ji ,jj, itest   
    1194       INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns   
    1195       REAL(wp), DIMENSION(2) ::   ztestmask 
    1196       !!---------------------------------------------------------------------- 
    1197       ! 
    1198       IF (lwp) WRITE(numout,*) ' ' 
    1199       IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments' 
    1200       IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    1201       ! 
    1202       IF(lwp) WRITE(numout,*) 'Number of east  segments     : ', nbdysege 
    1203       IF(lwp) WRITE(numout,*) 'Number of west  segments     : ', nbdysegw 
    1204       IF(lwp) WRITE(numout,*) 'Number of north segments     : ', nbdysegn 
    1205       IF(lwp) WRITE(numout,*) 'Number of south segments     : ', nbdysegs 
    1206       ! 1. Check bounds 
    1207       !---------------- 
    1208       DO ib = 1, nbdysegn 
    1209          IF (lwp) WRITE(numout,*) '**check north seg bounds pckg: ', npckgn(ib) 
    1210          IF ((jpjnob(ib).ge.jpjglo-1).or.&  
    1211             &(jpjnob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' ) 
    1212          IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 
    1213          IF (jpindt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' ) 
    1214          IF (jpinft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' ) 
    1215       END DO 
    1216       ! 
    1217       DO ib = 1, nbdysegs 
    1218          IF (lwp) WRITE(numout,*) '**check south seg bounds pckg: ', npckgs(ib) 
    1219          IF ((jpjsob(ib).ge.jpjglo-1).or.&  
    1220             &(jpjsob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' ) 
    1221          IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 
    1222          IF (jpisdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' ) 
    1223          IF (jpisft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' ) 
    1224       END DO 
    1225       ! 
    1226       DO ib = 1, nbdysege 
    1227          IF (lwp) WRITE(numout,*) '**check east  seg bounds pckg: ', npckge(ib) 
    1228          IF ((jpieob(ib).ge.jpiglo-1).or.&  
    1229             &(jpieob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' ) 
    1230          IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 
    1231          IF (jpjedt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' ) 
    1232          IF (jpjeft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' ) 
    1233       END DO 
    1234       ! 
    1235       DO ib = 1, nbdysegw 
    1236          IF (lwp) WRITE(numout,*) '**check west  seg bounds pckg: ', npckgw(ib) 
    1237          IF ((jpiwob(ib).ge.jpiglo-1).or.&  
    1238             &(jpiwob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' ) 
    1239          IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' ) 
    1240          IF (jpjwdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' ) 
    1241          IF (jpjwft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' ) 
    1242       ENDDO 
    1243       ! 
    1244       !       
    1245       ! 2. Look for segment crossings 
    1246       !------------------------------  
    1247       IF (lwp) WRITE(numout,*) '**Look for segments corners  :' 
    1248       ! 
    1249       itest = 0 ! corner number 
    1250       ! 
    1251       ! flag to detect if start or end of open boundary belongs to a corner 
    1252       ! if not (=0), it must be on land. 
    1253       ! if a corner is detected, save bdy package number for further tests 
    1254       icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0. 
    1255       ! South/West crossings 
    1256       IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN 
    1257          DO ib1 = 1, nbdysegw         
    1258             DO ib2 = 1, nbdysegs 
    1259                IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. & 
    1260                 &  ( jpisft(ib2)>=jpiwob(ib1)).AND. & 
    1261                 &  ( jpjwdt(ib1)<=jpjsob(ib2)).AND. & 
    1262                 &  ( jpjwft(ib1)>=jpjsob(ib2))) THEN 
    1263                   IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN  
    1264                      ! We have a possible South-West corner                       
    1265 !                     WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)  
    1266 !                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2) 
    1267                      icornw(ib1,1) = npckgs(ib2) 
    1268                      icorns(ib2,1) = npckgw(ib1) 
    1269                   ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN 
    1270                      IF(lwp) WRITE(numout,*) 
    1271                      IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
    1272                      &                                     jpisft(ib2), jpjwft(ib1) 
    1273                      IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet' 
    1274                      IF(lwp) WRITE(numout,*) '             Crossing problem with West segment: ',npckgw(ib1), &  
    1275                      &                                                    ' and South segment: ',npckgs(ib2) 
    1276                      IF(lwp) WRITE(numout,*) 
    1277                      nstop = nstop + 1 
    1278                   ELSE 
    1279                      IF(lwp) WRITE(numout,*) 
    1280                      IF(lwp) WRITE(numout,*) ' E R R O R : Check South and West Open boundary indices' 
    1281                      IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1) , & 
    1282                      &                                                    ' and South segment: ',npckgs(ib2) 
    1283                      IF(lwp) WRITE(numout,*) 
    1284                      nstop = nstop+1 
    1285                   END IF 
    1286                END IF 
    1287             END DO 
    1288          END DO 
    1289       END IF 
    1290       ! 
    1291       ! South/East crossings 
    1292       IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN 
    1293          DO ib1 = 1, nbdysege 
    1294             DO ib2 = 1, nbdysegs 
    1295                IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. & 
    1296                 &  ( jpisft(ib2)>=jpieob(ib1)+1).AND. & 
    1297                 &  ( jpjedt(ib1)<=jpjsob(ib2)  ).AND. & 
    1298                 &  ( jpjeft(ib1)>=jpjsob(ib2)  )) THEN 
    1299                   IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN 
    1300                      ! We have a possible South-East corner  
    1301 !                     WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)  
    1302 !                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2) 
    1303                      icorne(ib1,1) = npckgs(ib2) 
    1304                      icorns(ib2,2) = npckge(ib1) 
    1305                   ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN 
    1306                      IF(lwp) WRITE(numout,*) 
    1307                      IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
    1308                      &                                     jpisdt(ib2), jpjeft(ib1) 
    1309                      IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet' 
    1310                      IF(lwp) WRITE(numout,*) '             Crossing problem with East segment: ',npckge(ib1), & 
    1311                      &                                                    ' and South segment: ',npckgs(ib2) 
    1312                      IF(lwp) WRITE(numout,*) 
    1313                      nstop = nstop + 1 
    1314                   ELSE 
    1315                      IF(lwp) WRITE(numout,*) 
    1316                      IF(lwp) WRITE(numout,*) ' E R R O R : Check South and East Open boundary indices' 
    1317                      IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), & 
    1318                      &                                                    ' and South segment: ',npckgs(ib2) 
    1319                      IF(lwp) WRITE(numout,*) 
    1320                      nstop = nstop + 1 
    1321                   END IF 
    1322                END IF 
    1323             END DO 
    1324          END DO 
    1325       END IF 
    1326       ! 
    1327       ! North/West crossings 
    1328       IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN 
    1329          DO ib1 = 1, nbdysegw         
    1330             DO ib2 = 1, nbdysegn 
    1331                IF (( jpindt(ib2)<=jpiwob(ib1)  ).AND. & 
    1332                 &  ( jpinft(ib2)>=jpiwob(ib1)  ).AND. & 
    1333                 &  ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. & 
    1334                 &  ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN 
    1335                   IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN 
    1336                      ! We have a possible North-West corner  
    1337 !                     WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)  
    1338 !                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2) 
    1339                      icornw(ib1,2) = npckgn(ib2) 
    1340                      icornn(ib2,1) = npckgw(ib1) 
    1341                   ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN 
    1342                      IF(lwp) WRITE(numout,*) 
    1343                      IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
    1344                      &                                     jpinft(ib2), jpjwdt(ib1) 
    1345                      IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet' 
    1346                      IF(lwp) WRITE(numout,*) '             Crossing problem with West segment: ',npckgw(ib1), & 
    1347                      &                                                    ' and North segment: ',npckgn(ib2) 
    1348                      IF(lwp) WRITE(numout,*) 
    1349                      nstop = nstop + 1 
    1350                   ELSE 
    1351                      IF(lwp) WRITE(numout,*) 
    1352                      IF(lwp) WRITE(numout,*) ' E R R O R : Check North and West Open boundary indices' 
    1353                      IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1), & 
    1354                      &                                                    ' and North segment: ',npckgn(ib2) 
    1355                      IF(lwp) WRITE(numout,*) 
    1356                      nstop = nstop + 1 
    1357                   END IF 
    1358                END IF 
    1359             END DO 
    1360          END DO 
    1361       END IF 
    1362       ! 
    1363       ! North/East crossings 
    1364       IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN 
    1365          DO ib1 = 1, nbdysege         
    1366             DO ib2 = 1, nbdysegn 
    1367                IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. & 
    1368                 &  ( jpinft(ib2)>=jpieob(ib1)+1).AND. & 
    1369                 &  ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. & 
    1370                 &  ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN 
    1371                   IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN 
    1372                      ! We have a possible North-East corner  
    1373 !                     WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1) 
    1374 !                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2) 
    1375                      icorne(ib1,2) = npckgn(ib2) 
    1376                      icornn(ib2,2) = npckge(ib1) 
    1377                   ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN 
    1378                      IF(lwp) WRITE(numout,*) 
    1379                      IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', & 
    1380                      &                                     jpindt(ib2), jpjedt(ib1) 
    1381                      IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet' 
    1382                      IF(lwp) WRITE(numout,*) '             Crossing problem with East segment: ',npckge(ib1), & 
    1383                      &                                                    ' and North segment: ',npckgn(ib2) 
    1384                      IF(lwp) WRITE(numout,*) 
    1385                      nstop = nstop + 1 
    1386                   ELSE 
    1387                      IF(lwp) WRITE(numout,*) 
    1388                      IF(lwp) WRITE(numout,*) ' E R R O R : Check North and East Open boundary indices' 
    1389                      IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), & 
    1390                      &                                                    ' and North segment: ',npckgn(ib2) 
    1391                      IF(lwp) WRITE(numout,*) 
    1392                      nstop = nstop + 1 
    1393                   END IF 
    1394                END IF 
    1395             END DO 
    1396          END DO 
    1397       END IF 
    1398       ! 
    1399       ! 3. Check if segment extremities are on land 
    1400       !--------------------------------------------  
    1401       ! 
    1402       ! West segments 
    1403       DO ib = 1, nbdysegw 
    1404          ! get mask at boundary extremities: 
    1405          ztestmask(1:2)=0. 
    1406          DO ji = 1, jpi 
    1407             DO jj = 1, jpj              
    1408               IF (((ji + nimpp - 1) == jpiwob(ib)).AND. &  
    1409                &  ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1) 
    1410               IF (((ji + nimpp - 1) == jpiwob(ib)).AND. &  
    1411                &  ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1)   
    1412             END DO 
    1413          END DO 
    1414          IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain 
    1415  
    1416          IF (ztestmask(1)==1) THEN  
    1417             IF (icornw(ib,1)==0) THEN 
    1418                IF(lwp) WRITE(numout,*) 
    1419                IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib) 
    1420                IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                   
    1421                IF(lwp) WRITE(numout,*) 
    1422                nstop = nstop + 1 
    1423             ELSE 
    1424                ! This is a corner 
    1425                WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib) 
    1426                CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1)) 
    1427                itest=itest+1 
    1428             ENDIF 
    1429          ENDIF 
    1430          IF (ztestmask(2)==1) THEN 
    1431             IF (icornw(ib,2)==0) THEN 
    1432                IF(lwp) WRITE(numout,*) 
    1433                IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib) 
    1434                IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                   
    1435                IF(lwp) WRITE(numout,*) 
    1436                nstop = nstop + 1 
    1437             ELSE 
    1438                ! This is a corner 
    1439                WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib) 
    1440                CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2)) 
    1441                itest=itest+1 
    1442             ENDIF 
    1443          ENDIF 
    1444       END DO 
    1445       ! 
    1446       ! East segments 
    1447       DO ib = 1, nbdysege 
    1448          ! get mask at boundary extremities: 
    1449          ztestmask(1:2)=0. 
    1450          DO ji = 1, jpi 
    1451             DO jj = 1, jpj              
    1452               IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. &  
    1453                &  ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1) 
    1454               IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. &  
    1455                &  ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1)   
    1456             END DO 
    1457          END DO 
    1458          IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain 
    1459  
    1460          IF (ztestmask(1)==1) THEN 
    1461             IF (icorne(ib,1)==0) THEN 
    1462                IF(lwp) WRITE(numout,*) 
    1463                IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib) 
    1464                IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                   
    1465                IF(lwp) WRITE(numout,*) 
    1466                nstop = nstop + 1  
    1467             ELSE 
    1468                ! This is a corner 
    1469                WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib) 
    1470                CALL bdy_ctl_corn(npckge(ib), icorne(ib,1)) 
    1471                itest=itest+1 
    1472             ENDIF 
    1473          ENDIF 
    1474          IF (ztestmask(2)==1) THEN 
    1475             IF (icorne(ib,2)==0) THEN 
    1476                IF(lwp) WRITE(numout,*) 
    1477                IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib) 
    1478                IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                   
    1479                IF(lwp) WRITE(numout,*) 
    1480                nstop = nstop + 1 
    1481             ELSE 
    1482                ! This is a corner 
    1483                WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib) 
    1484                CALL bdy_ctl_corn(npckge(ib), icorne(ib,2)) 
    1485                itest=itest+1 
    1486             ENDIF 
    1487          ENDIF 
    1488       END DO 
    1489       ! 
    1490       ! South segments 
    1491       DO ib = 1, nbdysegs 
    1492          ! get mask at boundary extremities: 
    1493          ztestmask(1:2)=0. 
    1494          DO ji = 1, jpi 
    1495             DO jj = 1, jpj              
    1496               IF (((jj + njmpp - 1) == jpjsob(ib)).AND. &  
    1497                &  ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1) 
    1498               IF (((jj + njmpp - 1) == jpjsob(ib)).AND. &  
    1499                &  ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1)   
    1500             END DO 
    1501          END DO 
    1502          IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain 
    1503  
    1504          IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN 
    1505             IF(lwp) WRITE(numout,*) 
    1506             IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib) 
    1507             IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                   
    1508             IF(lwp) WRITE(numout,*) 
    1509             nstop = nstop + 1 
    1510          ENDIF 
    1511          IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN 
    1512             IF(lwp) WRITE(numout,*) 
    1513             IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib) 
    1514             IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                   
    1515             IF(lwp) WRITE(numout,*) 
    1516             nstop = nstop + 1 
    1517          ENDIF 
    1518       END DO 
    1519       ! 
    1520       ! North segments 
    1521       DO ib = 1, nbdysegn 
    1522          ! get mask at boundary extremities: 
    1523          ztestmask(1:2)=0. 
    1524          DO ji = 1, jpi 
    1525             DO jj = 1, jpj              
    1526               IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. &  
    1527                &  ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1) 
    1528               IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. &  
    1529                &  ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1)   
    1530             END DO 
    1531          END DO 
    1532          IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain 
    1533  
    1534          IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN 
    1535             IF(lwp) WRITE(numout,*) 
    1536             IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib) 
    1537             IF(lwp) WRITE(numout,*) ' ==========  does not start on land'                                                   
    1538             IF(lwp) WRITE(numout,*) 
    1539             nstop = nstop + 1 
    1540          ENDIF 
    1541          IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN 
    1542             IF(lwp) WRITE(numout,*) 
    1543             IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib) 
    1544             IF(lwp) WRITE(numout,*) ' ==========  does not end on land'                                                   
    1545             IF(lwp) WRITE(numout,*) 
    1546             nstop = nstop + 1 
    1547          ENDIF 
    1548       END DO 
    1549       ! 
    1550       IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found' 
    1551       ! 
    1552       ! Other tests TBD:  
    1553       ! segments completly on land 
    1554       ! optimized open boundary array length according to landmask 
    1555       ! Nudging layers that overlap with interior domain 
    1556       ! 
    1557    END SUBROUTINE bdy_ctl_seg 
    1558  
    1559    SUBROUTINE bdy_ctl_corn( ib1, ib2 ) 
    1560       !!---------------------------------------------------------------------- 
    1561       !!                 ***  ROUTINE bdy_ctl_corn  *** 
    1562       !! 
    1563       !! ** Purpose :   Check numerical schemes consistency between 
    1564       !!                segments having a common corner 
    1565       !! 
    1566       !! ** Method  :    
    1567       !!---------------------------------------------------------------------- 
    1568       INTEGER, INTENT(in)  ::   ib1, ib2 
    1569       INTEGER :: itest 
    1570       !!---------------------------------------------------------------------- 
    1571       itest = 0 
    1572  
    1573       IF (nn_dyn2d(ib1)/=nn_dyn2d(ib2)) itest = itest + 1 
    1574       IF (nn_dyn3d(ib1)/=nn_dyn3d(ib2)) itest = itest + 1 
    1575       IF (nn_tra(ib1)/=nn_tra(ib2)) itest = itest + 1 
    1576       ! 
    1577       IF (nn_dyn2d_dta(ib1)/=nn_dyn2d_dta(ib2)) itest = itest + 1 
    1578       IF (nn_dyn3d_dta(ib1)/=nn_dyn3d_dta(ib2)) itest = itest + 1 
    1579       IF (nn_tra_dta(ib1)/=nn_tra_dta(ib2)) itest = itest + 1 
    1580       ! 
    1581       IF (nn_rimwidth(ib1)/=nn_rimwidth(ib2)) itest = itest + 1    
    1582       ! 
    1583       IF ( itest>0 ) THEN 
    1584          IF(lwp) WRITE(numout,*) ' E R R O R : Segments ', ib1, 'and ', ib2 
    1585          IF(lwp) WRITE(numout,*) ' ==========  have different open bdy schemes'                                                   
    1586          IF(lwp) WRITE(numout,*) 
    1587          nstop = nstop + 1 
    1588       ENDIF 
    1589       ! 
    1590    END SUBROUTINE bdy_ctl_corn 
    1591819 
    1592820#else 

  • branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

    r3651 r6736  
    88   !!            3.0  !  2008-04  (NEMO team)  add in the reference version 
    99   !!            3.3  !  2010-09  (D.Storkey and E.O'Dea)  bug fixes 
    10    !!            3.4  !  2012-09  (G. Reffray and J. Chanut) New inputs + mods 
     10   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
     11   !!            3.4  !  2013     (J. Harle) rewite to used tide_mod for phase and nodal 
     12   !!                             corrections every day 
    1113   !!---------------------------------------------------------------------- 
    1214#if defined key_bdy 
    … …  
    1517   !!---------------------------------------------------------------------- 
    1618   !!   PUBLIC 
    17    !!      bdytide_init     : read of namelist and initialisation of tidal harmonics data 
     19   !!      tide_init     : read of namelist and initialisation of tidal harmonics data 
    1820   !!      tide_update   : calculation of tidal forcing at each timestep 
    1921   !!---------------------------------------------------------------------- 
    … …  
    2729   USE bdy_par         ! Unstructured boundary parameters 
    2830   USE bdy_oce         ! ocean open boundary conditions 
     31   USE fldread, ONLY: fld_map 
    2932   USE daymod          ! calendar 
    30    USE wrk_nemo        ! Memory allocation 
    31    USE tideini 
    32 !   USE tide_mod       ! Useless ?? 
    33    USE fldread, ONLY: fld_map 
     33   USE tide_mod 
     34   USE ioipsl, ONLY :   ymds2ju   ! for calendar 
    3435 
    3536   IMPLICIT NONE 
    3637   PRIVATE 
    3738 
    38    PUBLIC   bdytide_init     ! routine called in bdy_init 
    39    PUBLIC   bdytide_update   ! routine called in bdy_dta 
     39   PUBLIC   tide_init     ! routine called in nemo_init 
     40   PUBLIC   tide_update   ! routine called in bdydyn 
    4041 
    4142   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data 
    42       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh0       !: Tidal constituents : SSH0 (read in file) 
    43       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u0         !: Tidal constituents : U0   (read in file) 
    44       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v0         !: Tidal constituents : V0   (read in file) 
    45       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh        !: Tidal constituents : SSH  (after nodal cor.) 
    46       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u          !: Tidal constituents : U    (after nodal cor.) 
    47       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V    (after nodal cor.) 
     43      INTEGER                                ::   ncpt       !: Actual number of tidal components 
     44      REAL(wp), POINTER, DIMENSION(:)        ::   speed      !: Phase speed of tidal constituent (deg/hr) 
     45      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh        !: Tidal constituents : SSH 
     46      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u          !: Tidal constituents : U 
     47      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V 
     48      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   sshr       !: Tidal constituents : SSH (reference) 
     49      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ur         !: Tidal constituents : U (reference) 
     50      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   vr         !: Tidal constituents : V (reference) 
    4851   END TYPE TIDES_DATA 
    4952 
    50    TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data 
    51  
     53   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET ::   tides                 !: External tidal harmonics data 
     54 
     55   INTEGER, ALLOCATABLE, DIMENSION(:)  :: bdy_ntide 
     56   REAL(wp), ALLOCATABLE, DIMENSION(:) :: bdy_omega_tide 
     57   REAL(wp), ALLOCATABLE, DIMENSION(:) :: bdy_v0tide,      & 
     58                                          bdy_blank,       & 
     59                                          bdy_utide,       & 
     60                                          bdy_ftide,       & 
     61                                          rbdy_ftide 
     62   LOGICAL                             ::   ln_tide_date !: =T correct tide phases and amplitude for model start date 
     63   LOGICAL                             ::   ln_tide_v0 !: =T correct tide phases and amplitude for model start date 
     64   INTEGER                             ::   nn_tide_date !: yyyymmdd reference date of tidal data 
     65   INTEGER ::   bdy_nn_tide 
     66   INTEGER ::   bdy_kt_tide      ! Main tide timestep counter 
     67   INTEGER ::   bdy_tide_offset      ! Main tide timestep counter 
     68    
    5269   !!---------------------------------------------------------------------- 
    5370   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    … …  
    5774CONTAINS 
    5875 
    59    SUBROUTINE bdytide_init 
    60       !!---------------------------------------------------------------------- 
    61       !!                    ***  SUBROUTINE bdytide_init  *** 
     76   SUBROUTINE tide_init 
     77      !!---------------------------------------------------------------------- 
     78      !!                    ***  SUBROUTINE tide_init  *** 
    6279      !!                      
    6380      !! ** Purpose : - Read in namelist for tides and initialise external 
    … …  
    6784      !! namelist variables 
    6885      !!------------------- 
    69       CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files 
    70       LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data 
    71       LOGICAL                                   ::   ln_bdytide_conj     !: If true, assume complex conjugate tidal data 
     86      CHARACTER(len=80)                         ::   filtide      !: Filename root for tidal input files 
     87      CHARACTER(len= 4), DIMENSION(jpmax_harmo) ::   tide_cpt     !: Names of tidal components used. 
    7288      !! 
    73       INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices 
    74       INTEGER                                   ::   ii, ij              !: dummy loop indices 
     89      INTEGER                                   ::   ib_bdy, itide, ib, ji  !: dummy loop indices 
    7590      INTEGER                                   ::   inum, igrd 
    76       INTEGER, DIMENSION(3)                     ::   ilen0       !: length of boundary data (from OBC arrays) 
    77       INTEGER, POINTER, DIMENSION(:)            ::   nblen, nblenrim     ! short cuts 
    78       CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file  
    79       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            !: work space to read in tidal harmonics data 
    80       REAL(wp), POINTER, DIMENSION(:,:)         ::   ztr, zti            !:  "     "    "   "   "   "        "      "  
     91      INTEGER  :: lcl_ryear, lcl_rmonth, lcl_rday 
     92      INTEGER, DIMENSION(3)                     ::   ilen0                  !: length of boundary data (from OBC arrays) 
     93      CHARACTER(len=80)                         ::   clfile                 !: full file name for tidal input file  
     94      REAL(wp)                                  ::   z_arg, z_atde, z_btde, z1t, z2t, fdayn, fdayr 
     95      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read           !: work space to read in tidal harmonics data 
    8196      !! 
    82       TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut    
     97      TYPE(TIDES_DATA),  POINTER                ::   td                 !: local short cut    
    8398      !! 
    84       NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj 
    85       !!---------------------------------------------------------------------- 
    86  
    87       IF( nn_timing == 1 ) CALL timing_start('bdytide_init') 
    88  
    89       IF (nb_bdy>0) THEN 
    90          IF(lwp) WRITE(numout,*) 
    91          IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries' 
    92          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 
    93       ENDIF 
    94  
    95       ln_bdytide_2ddta = .FALSE. 
    96       ln_bdytide_conj  = .FALSE. 
     99      NAMELIST/nambdy_tide/filtide, tide_cpt, ln_tide_date, nn_tide_date, ln_tide_v0 
     100      !!---------------------------------------------------------------------- 
     101 
     102      IF( nn_timing == 1 ) CALL timing_start('tide_init') 
     103 
     104      IF(lwp) WRITE(numout,*) 
     105      IF(lwp) WRITE(numout,*) 'tide_init : initialization of tidal harmonic forcing at open boundaries' 
     106      IF(lwp) WRITE(numout,*) '~~~~~~~~~' 
    97107 
    98108      REWIND(numnam) 
    … …  
    101111 
    102112            td => tides(ib_bdy) 
    103             nblen => idx_bdy(ib_bdy)%nblen 
    104             nblenrim => idx_bdy(ib_bdy)%nblenrim 
    105113 
    106114            ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries 
     115            ln_tide_date = .false. 
     116            ln_tide_v0 = .false. 
    107117            filtide(:) = '' 
     118            tide_cpt(:) = '' 
     119 
     120            ! Initialise bdy_ky_tide: updated in tide_update if using time correction otherwise defaults to 1 
     121            bdy_kt_tide=1 
    108122 
    109123            ! Don't REWIND here - may need to read more than one of these namelists. 
    110124            READ  ( numnam, nambdy_tide ) 
     125            !                                               ! Count number of components specified 
     126            td%ncpt = 0 
     127            DO itide = 1, jpmax_harmo 
     128              IF( tide_cpt(itide) /= '' ) THEN 
     129                 td%ncpt = td%ncpt + 1 
     130              ENDIF 
     131            END DO 
     132 
     133            CALL tide_init_Wave 
     134 
     135            ! Find constituents in standard list 
     136            ALLOCATE(bdy_ntide     (td%ncpt)) 
     137     
     138            DO itide=1,td%ncpt 
     139               bdy_ntide(itide)=0 
     140               DO ji=1,jpmax_harmo 
     141                  IF ( TRIM( tide_cpt(itide) ) .eq. Wave(ji)%cname_tide) THEN 
     142                     bdy_ntide(itide) = ji 
     143                     EXIT 
     144                  END IF 
     145               END DO 
     146               IF (bdy_ntide(itide).eq.0) THEN 
     147                  CALL ctl_stop( 'BDYTIDE tidal components do not match up with tide.h90' ) 
     148               ENDIF 
     149            END DO 
     150     
     151            ! Fill in phase speeds from tide_pulse 
     152            ALLOCATE(bdy_omega_tide(td%ncpt)) 
     153            CALL tide_pulse( bdy_omega_tide, bdy_ntide ,td%ncpt) 
     154 
     155            ALLOCATE( td%speed(td%ncpt) ) 
     156            td%speed = bdy_omega_tide(1:td%ncpt) 
     157 
    111158            !                                               ! Parameter control and print 
    112             IF(lwp) WRITE(numout,*) '  ' 
    113             IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries' 
    114             IF(lwp) WRITE(numout,*) '             read tidal data in 2d files: ', ln_bdytide_2ddta 
    115             IF(lwp) WRITE(numout,*) '             assume complex conjugate   : ', ln_bdytide_conj 
    116             IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo 
    117             IF(lwp) THEN  
    118                     WRITE(numout,*) '             Tidal cpt name    -     Phase speed (deg/hr)'             
    119                DO itide = 1, nb_harmo 
    120                   WRITE(numout,*)  '             ', Wave(ntide(itide))%cname_tide, omega_tide(itide)    
    121                END DO 
    122             ENDIF  
    123             IF(lwp) WRITE(numout,*) ' ' 
    124  
    125             ! Allocate space for tidal harmonics data - get size from OBC data arrays 
    126             ! ----------------------------------------------------------------------- 
    127  
    128             ! JC: If FRS scheme is used, we assume that tidal is needed over the whole 
    129             ! relaxation area       
    130             IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
    131                ilen0(:)=nblen(:) 
     159            IF( td%ncpt < 1 ) THEN  
     160               CALL ctl_stop( '          Did not find any tidal components in namelist nambdy_tide' ) 
    132161            ELSE 
    133                ilen0(:)=nblenrim(:) 
     162               IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries' 
     163               IF(lwp) WRITE(numout,*) '             tidal components specified ', td%ncpt 
     164               IF(lwp) WRITE(numout,*) '                ', tide_cpt(1:td%ncpt) 
     165               IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : ' 
     166               IF(lwp) WRITE(numout,*) '                ', td%speed(1:td%ncpt) 
    134167            ENDIF 
    135168 
    136             ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) ) 
    137             ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) ) 
    138  
    139             ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) ) 
    140             ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) ) 
    141  
    142             ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) ) 
    143             ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) ) 
    144  
    145             td%ssh0(:,:,:) = 0.e0 
    146             td%ssh(:,:,:) = 0.e0 
    147             td%u0(:,:,:) = 0.e0 
    148             td%u(:,:,:) = 0.e0 
    149             td%v0(:,:,:) = 0.e0 
    150             td%v(:,:,:) = 0.e0 
    151  
    152             IF (ln_bdytide_2ddta) THEN 
    153                ! It is assumed that each data file contains all complex harmonic amplitudes 
    154                ! given on the data domain (ie global, jpidta x jpjdta) 
    155                ! 
    156                CALL wrk_alloc( jpi, jpj, zti, ztr ) 
    157                ! 
    158                ! SSH fields 
    159                clfile = TRIM(filtide)//'_grid_T.nc' 
    160                CALL iom_open (clfile , inum )  
    161                igrd = 1                       ! Everything is at T-points here 
    162                DO itide = 1, nb_harmo 
    163                   CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) ) 
    164                   CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) )  
    165                   DO ib = 1, ilen0(igrd) 
    166                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
    167                      ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
    168                      td%ssh0(ib,itide,1) = ztr(ii,ij) 
    169                      td%ssh0(ib,itide,2) = zti(ii,ij) 
    170                   END DO 
    171                END DO  
     169            ! Allocate space for tidal harmonics data -  
     170            ! get size from OBC data arrays 
     171            ! --------------------------------------- 
     172 
     173            ilen0(1) = SIZE( dta_bdy(ib_bdy)%ssh )  
     174            ALLOCATE( td%ssh( ilen0(1), td%ncpt, 2 ) ) 
     175            ALLOCATE( td%sshr( ilen0(1), td%ncpt, 2 ) ) 
     176 
     177            ilen0(2) = SIZE( dta_bdy(ib_bdy)%u2d )  
     178            ALLOCATE( td%u( ilen0(2), td%ncpt, 2 ) ) 
     179            ALLOCATE( td%ur( ilen0(2), td%ncpt, 2 ) ) 
     180 
     181            ilen0(3) = SIZE( dta_bdy(ib_bdy)%v2d )  
     182            ALLOCATE( td%v( ilen0(3), td%ncpt, 2 ) ) 
     183            ALLOCATE( td%vr( ilen0(3), td%ncpt, 2 ) ) 
     184 
     185            ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) ) 
     186  
     187            ! Set day length in timesteps for use if making phase and nodal corrections 
     188            bdy_nn_tide=NINT(rday/rdt) 
     189             
     190  
     191            ALLOCATE(bdy_v0tide   (td%ncpt)) 
     192            ALLOCATE(bdy_blank   (td%ncpt)) 
     193            ALLOCATE(bdy_utide    (td%ncpt)) 
     194            ALLOCATE(bdy_ftide    (td%ncpt)) 
     195            ALLOCATE(rbdy_ftide    (td%ncpt)) 
     196       
     197            ! Open files and read in tidal forcing data 
     198            ! ----------------------------------------- 
     199 
     200            DO itide = 1, td%ncpt 
     201               !                                                              ! SSH fields 
     202               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc' 
     203               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
     204               CALL iom_open( clfile, inum ) 
     205               CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
     206               td%ssh(:,itide,1) = dta_read(1:ilen0(1),1,1) 
     207               CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
     208               td%ssh(:,itide,2) = dta_read(1:ilen0(1),1,1) 
     209               CALL iom_close( inum ) 
     210               !                                                              ! U fields 
     211               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc' 
     212               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
     213               CALL iom_open( clfile, inum ) 
     214               CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
     215               td%u(:,itide,1) = dta_read(1:ilen0(2),1,1) 
     216               CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
     217               td%u(:,itide,2) = dta_read(1:ilen0(2),1,1) 
     218               CALL iom_close( inum ) 
     219               !                                                              ! V fields 
     220               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc' 
     221               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 
     222               CALL iom_open( clfile, inum ) 
     223               CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
     224               td%v(:,itide,1) = dta_read(1:ilen0(3),1,1) 
     225               CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
     226               td%v(:,itide,2) = dta_read(1:ilen0(3),1,1) 
    172227               CALL iom_close( inum ) 
    173228               ! 
    174             &