Changeset 5620


Ignore:
Timestamp:
2015-07-21T10:55:28+02:00 (5 years ago)
Author:
jamesharle
Message:

Merge with r5619 of trunk, update to unstructured BDY interpolation in
fldread.F90. Structured BDY interpolation incomplete.

Location:
branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC
Files:
155 edited
1 copied

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90

    • Property svn:keywords set to Id
    r5038 r5620  
    5757   !!---------------------------------------------------------------------- 
    5858   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    59    !! $Id: $ 
     59   !! $Id$ 
    6060   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    6161   !!---------------------------------------------------------------------- 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r5038 r5620  
    658658 
    659659      DO jk = 1, jpkm1 
    660          fzptnz(:,:,jk) = eos_fzp( tsn(:,:,jk,jp_sal), fsdept(:,:,jk) ) 
     660        CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), fsdept(:,:,jk) ) 
    661661      END DO 
    662662 
     
    746746 
    747747 
    748             IF( ln_zps .AND. .NOT. lk_c1d ) & 
    749                &  CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
    750                &                rhd, gru , grv, aru, arv, gzu, gzv, ge3ru, ge3rv,  &             ! 
    751                &                gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
     748            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)      & 
     749               &  CALL zps_hde    ( kt, jpts, tsb, gtsu, gtsv,        &  ! Partial steps: before horizontal gradient 
     750               &                              rhd, gru , grv          )  ! of t, s, rd at the last ocean level 
     751            IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)      & 
     752               &  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,    &    ! Partial steps for top cell (ISF) 
     753               &                                  rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     754               &                           gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
    752755 
    753756#if defined key_zdfkpp 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90

    r4792 r5620  
    4949      LOGICAL                         ::  ll_tem 
    5050      LOGICAL                         ::  ll_sal 
     51      LOGICAL                         ::  ll_fvl 
    5152      REAL(wp), POINTER, DIMENSION(:)     ::  ssh 
    5253      REAL(wp), POINTER, DIMENSION(:)     ::  u2d 
     
    130131   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays (unstr.  bdy) 
    131132   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global_z      !: workspace for reading in global depth arrays (unstr.  bdy) 
     133   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global_dz      !: workspace for reading in global depth arrays (unstr.  bdy) 
    132134   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2       !: workspace for reading in global data arrays (struct. bdy) 
    133135   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2_z     !: workspace for reading in global depth arrays (struct. bdy) 
     136   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2_dz     !: workspace for reading in global depth arrays (struct. bdy) 
    134137!$AGRIF_DO_NOT_TREAT 
    135138   TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r4792 r5620  
    3333   USE ice_2 
    3434#elif defined key_lim3 
    35    USE par_ice 
    3635   USE ice 
    37    USE limcat_1D          ! redistribute ice input into categories 
     36   USE limvar          ! redistribute ice input into categories 
    3837#endif 
    3938   USE sbcapr 
     
    275274 
    276275                        jend = jstart + dta%nread(2) - 1 
    277                         CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  & 
    278                                      & kit=jit, kt_offset=time_offset ) 
     276                        IF( ln_full_vel_array(ib_bdy) ) THEN 
     277                           CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  & 
     278                                     & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(ib_bdy)  ) 
     279                        ELSE 
     280                           CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  & 
     281                                     & kit=jit, kt_offset=time_offset  ) 
     282                        ENDIF 
    279283 
    280284                        ! If full velocities in boundary data then extract barotropic velocities from 3D fields 
     
    341345                     jend = jstart + dta%nread(1) - 1 
    342346                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), & 
    343                                   & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy ) 
     347                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(ib_bdy) ) 
    344348                  ENDIF 
    345349                  ! If full velocities in boundary data then split into barotropic and baroclinic data 
     
    380384#if defined key_lim3 
    381385               IF( .NOT. ll_bdylim3 .AND. cn_ice_lim(ib_bdy) /= 'none' .AND. nn_ice_lim_dta(ib_bdy) == 1 ) THEN ! bdy ice input (case input is lim2 type) 
    382                 CALL lim_cat_1D ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), & 
     386                CALL lim_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), & 
    383387                                  & dta_bdy(ib_bdy)%ht_i,     dta_bdy(ib_bdy)%ht_s,     dta_bdy(ib_bdy)%a_i     ) 
    384388               ENDIF 
     
    536540            cn_dir_array(ib_bdy) = cn_dir 
    537541            ln_full_vel_array(ib_bdy) = ln_full_vel 
     542            !dta%ll_fvl = ln_full_vel ! jdha need this in fldread routine to work out what type of correction to apply to interpolated bdy data (maybe we replace all instances of ln_full_vel_array with this rather than duplicate) 
    538543 
    539544            nblen => idx_bdy(ib_bdy)%nblen 
     
    734739         IF( blf_i(jfld)%ln_tint ) ALLOCATE( bf(jfld)%fdta(ilen1(jfld),1,ilen3(jfld),2) ) 
    735740         nbmap_ptr(jfld)%ptr => idx_bdy(ibdy(jfld))%nbmap(:,igrid(jfld)) 
     741         nbmap_ptr(jfld)%ll_unstruc = ln_coords_file(ibdy(jfld)) 
    736742      ENDDO 
    737743 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90

    • Property svn:keywords set to Id
    r5038 r5620  
    3636   !!---------------------------------------------------------------------- 
    3737   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    38    !! $Id: bdydyn.F90 2528 2010-12-27 17:33:53Z rblod $  
     38   !! $Id$  
    3939   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    4040   !!---------------------------------------------------------------------- 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90

    • Property svn:keywords set to Id
    r4354 r5620  
    3333   !!---------------------------------------------------------------------- 
    3434   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    35    !! $Id: bdydyn.F90 2528 2010-12-27 17:33:53Z rblod $  
     35   !! $Id$  
    3636   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    3737   !!---------------------------------------------------------------------- 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90

    • Property svn:keywords set to Id
    r5038 r5620  
    2626   USE dom_ice_2       ! sea-ice domain 
    2727#elif defined key_lim3 
    28    USE par_ice 
    2928   USE ice             ! LIM_3 ice variables 
    3029   USE dom_ice         ! sea-ice domain 
     30   USE limvar 
    3131#endif  
    3232   USE par_oce         ! ocean parameters 
     
    4242   PRIVATE 
    4343 
    44    PUBLIC   bdy_ice_lim    ! routine called in sbcmod 
     44   PUBLIC   bdy_ice_lim     ! routine called in sbcmod 
    4545   PUBLIC   bdy_ice_lim_dyn ! routine called in limrhg 
    4646 
    4747   !!---------------------------------------------------------------------- 
    4848   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    49    !! $Id: bdyice.F90 2715 2011-03-30 15:58:35Z rblod $ 
     49   !! $Id$ 
    5050   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    5151   !!---------------------------------------------------------------------- 
     
    6060      !!---------------------------------------------------------------------- 
    6161      INTEGER, INTENT( in ) :: kt     ! Main time step counter 
    62       !! 
    6362      INTEGER               :: ib_bdy ! Loop index 
     63 
     64#if defined key_lim3 
     65      CALL lim_var_glo2eqv 
     66#endif 
     67 
    6468      DO ib_bdy=1, nb_bdy 
    6569 
     
    7276            CALL ctl_stop( 'bdy_ice_lim : unrecognised option for open boundaries for ice fields' ) 
    7377         END SELECT 
    74       ENDDO 
     78 
     79      END DO 
     80 
     81#if defined key_lim3 
     82      CALL lim_var_zapsmall 
     83      CALL lim_var_agg(1) 
     84#endif 
    7585 
    7686   END SUBROUTINE bdy_ice_lim 
     
    8999      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data 
    90100      INTEGER,         INTENT(in) ::   kt   ! main time-step counter 
    91       INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index      !! 
     101      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index 
    92102 
    93103      INTEGER  ::   jpbound            ! 0 = incoming ice 
     
    169179            jpbound = 0; ii = ji; ij = jj; 
    170180 
    171             IF ( u_ice(ji+1,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. ) jpbound = 1; ii = ji+1; ij = jj 
    172             IF ( u_ice(ji-1,jj  ) > 0. .AND. umask(ji+1,jj  ,1) == 0. ) jpbound = 1; ii = ji-1; ij = jj 
    173             IF ( v_ice(ji  ,jj+1) < 0. .AND. vmask(ji  ,jj-1,1) == 0. ) jpbound = 1; ii = ji  ; ij = jj+1 
    174             IF ( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj+1,1) == 0. ) jpbound = 1; ii = ji  ; ij = jj-1 
    175  
    176             rswitch = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ii,ij) + 0.01 ) ) ! 0 if no ice 
     181            IF( u_ice(ji+1,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. ) jpbound = 1; ii = ji+1; ij = jj 
     182            IF( u_ice(ji-1,jj  ) > 0. .AND. umask(ji+1,jj  ,1) == 0. ) jpbound = 1; ii = ji-1; ij = jj 
     183            IF( v_ice(ji  ,jj+1) < 0. .AND. vmask(ji  ,jj-1,1) == 0. ) jpbound = 1; ii = ji  ; ij = jj+1 
     184            IF( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj+1,1) == 0. ) jpbound = 1; ii = ji  ; ij = jj-1 
     185 
     186            IF( nn_ice_lim_dta(ib_bdy) == 0 ) jpbound = 0; ii = ji; ij = jj   ! case ice boundaries = initial conditions 
     187                                                                              !      do not make state variables dependent on velocity 
     188                
     189 
     190            rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ii,ij) - 0.01 ) ) ! 0 if no ice 
    177191 
    178192            ! concentration and thickness 
     
    190204 
    191205               ! Ice salinity, age, temperature 
    192 <<<<<<< .working 
    193                sm_i(ji,jj,jl)   = zinda * rn_ice_sal(ib_bdy)  + ( 1.0 - zinda ) * s_i_min 
    194                o_i(ji,jj,jl)    = zinda * rn_ice_age(ib_bdy)  + ( 1.0 - zinda ) 
    195                t_su(ji,jj,jl)   = zinda * rn_ice_tem(ib_bdy)  + ( 1.0 - zinda ) * rn_ice_tem(ib_bdy) 
    196 ======= 
    197                sm_i(ji,jj,jl)   = rswitch * rn_ice_sal(ib_bdy)  + ( 1.0 - rswitch ) * s_i_min 
    198                o_i(ji,jj,jl)    = rswitch * rn_ice_age(ib_bdy)  + ( 1.0 - rswitch ) 
     206               sm_i(ji,jj,jl)   = rswitch * rn_ice_sal(ib_bdy)  + ( 1.0 - rswitch ) * rn_simin 
     207               oa_i(ji,jj,jl)   = rswitch * rn_ice_age(ib_bdy) * a_i(ji,jj,jl) 
    199208               t_su(ji,jj,jl)   = rswitch * rn_ice_tem(ib_bdy)  + ( 1.0 - rswitch ) * rn_ice_tem(ib_bdy) 
    200 >>>>>>> .merge-right.r5035 
    201209               DO jk = 1, nlay_s 
    202 <<<<<<< .working 
    203                   t_s(ji,jj,jk,jl) = zinda * rn_ice_tem(ib_bdy) + ( 1.0 - zinda ) * rtt 
    204 ======= 
    205                   t_s(ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rtt 
    206 >>>>>>> .merge-right.r5035 
     210                  t_s(ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rt0 
    207211               END DO 
    208212               DO jk = 1, nlay_i 
    209 <<<<<<< .working 
    210                   t_i(ji,jj,jk,jl) = zinda * rn_ice_tem(ib_bdy) + ( 1.0 - zinda ) * rtt  
    211                   s_i(ji,jj,jk,jl) = zinda * rn_ice_sal(ib_bdy) + ( 1.0 - zinda ) * s_i_min 
    212 ======= 
    213                   t_i(ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rtt  
    214                   s_i(ji,jj,jk,jl) = rswitch * rn_ice_sal(ib_bdy) + ( 1.0 - rswitch ) * s_i_min 
    215 >>>>>>> .merge-right.r5035 
     213                  t_i(ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rt0  
     214                  s_i(ji,jj,jk,jl) = rswitch * rn_ice_sal(ib_bdy) + ( 1.0 - rswitch ) * rn_simin 
    216215               END DO 
    217216                
     
    219218  
    220219               ! Ice salinity, age, temperature 
    221                sm_i(ji,jj,jl)   = rswitch * sm_i(ii,ij,jl)  + ( 1.0 - rswitch ) * s_i_min 
    222                o_i(ji,jj,jl)    = rswitch * o_i(ii,ij,jl)   + ( 1.0 - rswitch ) 
    223                t_su(ji,jj,jl)   = rswitch * t_su(ii,ij,jl)  + ( 1.0 - rswitch ) * rtt 
     220               sm_i(ji,jj,jl)   = rswitch * sm_i(ii,ij,jl)  + ( 1.0 - rswitch ) * rn_simin 
     221               oa_i(ji,jj,jl)   = rswitch * oa_i(ii,ij,jl) 
     222               t_su(ji,jj,jl)   = rswitch * t_su(ii,ij,jl)  + ( 1.0 - rswitch ) * rt0 
    224223               DO jk = 1, nlay_s 
    225                   t_s(ji,jj,jk,jl) = rswitch * t_s(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rtt 
     224                  t_s(ji,jj,jk,jl) = rswitch * t_s(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rt0 
    226225               END DO 
    227226               DO jk = 1, nlay_i 
    228                   t_i(ji,jj,jk,jl) = rswitch * t_i(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rtt 
    229                   s_i(ji,jj,jk,jl) = rswitch * s_i(ii,ij,jk,jl) + ( 1.0 - rswitch ) * s_i_min 
     227                  t_i(ji,jj,jk,jl) = rswitch * t_i(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rt0 
     228                  s_i(ji,jj,jk,jl) = rswitch * s_i(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rn_simin 
    230229               END DO 
    231230 
     
    233232 
    234233            ! if salinity is constant, then overwrite rn_ice_sal 
    235             IF( num_sal == 1 ) THEN 
    236                sm_i(ji,jj,jl)   = bulk_sal 
    237                s_i (ji,jj,:,jl) = bulk_sal 
     234            IF( nn_icesal == 1 ) THEN 
     235               sm_i(ji,jj,jl)   = rn_icesal 
     236               s_i (ji,jj,:,jl) = rn_icesal 
    238237            ENDIF 
    239238 
    240239            ! contents 
    241240            smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 
    242             oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl) 
    243241            DO jk = 1, nlay_s 
    244242               ! Snow energy of melting 
    245                e_s(ji,jj,jk,jl) = rswitch * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 
    246                ! Change dimensions 
    247                e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 
    248                ! Multiply by volume, so that heat content in 10^9 Joules 
    249                e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s 
     243               e_s(ji,jj,jk,jl) = rswitch * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus ) 
     244               ! Multiply by volume, so that heat content in J/m2 
     245               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s 
    250246            END DO 
    251247            DO jk = 1, nlay_i 
    252                ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K                   
     248               ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K                   
    253249               ! heat content per unit volume 
    254250               e_i(ji,jj,jk,jl) = rswitch * rhoic * & 
    255251                  (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
    256                   +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & 
    257                   - rcp      * ( ztmelts - rtt ) ) 
    258                ! Correct dimensions to avoid big values 
    259                e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac  
    260                ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
    261                e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / nlay_i 
     252                  +   lfus    * ( 1.0 - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) & 
     253                  - rcp      * ( ztmelts - rt0 ) ) 
     254               ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2 
     255               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) * r1_nlay_i 
    262256            END DO 
    263257 
    264  
    265          END DO !jb 
     258         END DO 
    266259  
    267          CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy )                                         ! lateral boundary conditions 
     260         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy ) 
    268261         CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy ) 
    269262         CALL lbc_bdy_lnk( ht_s(:,:,jl), 'T', 1., ib_bdy ) 
     
    274267         CALL lbc_bdy_lnk(  sm_i(:,:,jl), 'T', 1., ib_bdy ) 
    275268         CALL lbc_bdy_lnk(  oa_i(:,:,jl), 'T', 1., ib_bdy ) 
    276          CALL lbc_bdy_lnk(   o_i(:,:,jl), 'T', 1., ib_bdy ) 
    277269         CALL lbc_bdy_lnk(  t_su(:,:,jl), 'T', 1., ib_bdy ) 
    278270         DO jk = 1, nlay_s 
     
    306298      !! 
    307299      CHARACTER(len=1), INTENT(in)  ::   cd_type   ! nature of velocity grid-points 
    308       INTEGER  ::   jb, jgrd   ! dummy loop indices 
     300      INTEGER  ::   jb, jgrd           ! dummy loop indices 
    309301      INTEGER  ::   ji, jj             ! local scalar 
    310       INTEGER  ::   ib_bdy ! Loop index 
     302      INTEGER  ::   ib_bdy             ! Loop index 
    311303      REAL(wp) ::   zmsk1, zmsk2, zflag 
    312304     !!------------------------------------------------------------------------------ 
     
    324316         CASE('frs') 
    325317             
    326  
     318            IF( nn_ice_lim_dta(ib_bdy) == 0 ) CYCLE            ! case ice boundaries = initial conditions  
     319                                                               !      do not change ice velocity (it is only computed by rheology) 
     320  
    327321            SELECT CASE ( cd_type ) 
    328  
     322                
    329323            CASE ( 'U' ) 
    330324                
     
    341335                      
    342336                     ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce) 
    343                      u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + & 
    344                         &            u_ice(ji-1,jj) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + & 
     337                     u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + & 
     338                        &            u_ice(ji-1,jj) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + & 
    345339                        &            u_oce(ji  ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 
    346340                  ELSE                             ! everywhere else 
     
    349343                  ENDIF 
    350344                  ! mask ice velocities 
    351                   rswitch = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ji,jj) + 0.01 ) ) ! 0 if no ice 
     345                  rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ji,jj) - 0.01_wp ) ) ! 0 if no ice 
    352346                  u_ice(ji,jj) = rswitch * u_ice(ji,jj) 
    353347                   
    354348               ENDDO 
    355  
     349                
    356350               CALL lbc_bdy_lnk( u_ice(:,:), 'U', -1., ib_bdy ) 
    357351                
     
    370364                      
    371365                     ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce) 
    372                      v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5 * ABS( zflag + 1._wp ) * zmsk1 + & 
    373                         &            v_ice(ji,jj-1) * 0.5 * ABS( zflag - 1._wp ) * zmsk2 + & 
     366                     v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + & 
     367                        &            v_ice(ji,jj-1) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + & 
    374368                        &            v_oce(ji,jj  ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) ) 
    375369                  ELSE                             ! everywhere else 
     
    378372                  ENDIF 
    379373                  ! mask ice velocities 
    380                   rswitch = 1.0 - MAX( 0.0_wp , SIGN ( 1.0_wp , - at_i(ji,jj) + 0.01 ) ) ! 0 if no ice 
     374                  rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ji,jj) - 0.01 ) ) ! 0 if no ice 
    381375                  v_ice(ji,jj) = rswitch * v_ice(ji,jj) 
    382376                   
     
    384378                
    385379               CALL lbc_bdy_lnk( v_ice(:,:), 'V', -1., ib_bdy ) 
    386                 
     380                   
    387381            END SELECT 
    388382             
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r5038 r5620  
    155155        ENDIF 
    156156        IF(lwp) WRITE(numout,*) 
    157  
     157         
    158158        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  ' 
    159159        SELECT CASE( cn_dyn2d(ib_bdy) )                   
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90

    • Property svn:keywords set to Id
    r4292 r5620  
    2929   !!---------------------------------------------------------------------- 
    3030   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    31    !! $Id: bdydyn.F90 2528 2010-12-27 17:33:53Z rblod $  
     31   !! $Id$  
    3232   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    3333   !!---------------------------------------------------------------------- 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

    r4792 r5620  
    3232   USE tideini 
    3333!   USE tide_mod       ! Useless ?? 
    34    USE fldread, ONLY: fld_map 
     34   USE fldread 
    3535   USE dynspg_oce, ONLY: lk_dynspg_ts 
    3636 
     
    8888      !! 
    8989      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut    
     90      TYPE(MAP_POINTER), DIMENSION(jpbgrd)      ::   ibmap_ptr           !: array of pointers to nbmap 
    9091      !! 
    9192      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj 
     
    125126            IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo 
    126127            IF(lwp) THEN  
    127                     WRITE(numout,*) '             Tidal cpt name    -     Phase speed (deg/hr)'             
     128                    WRITE(numout,*) '             Tidal components: '  
    128129               DO itide = 1, nb_harmo 
    129                   WRITE(numout,*)  '             ', Wave(ntide(itide))%cname_tide, omega_tide(itide)   
     130                  WRITE(numout,*)  '                 ', Wave(ntide(itide))%cname_tide  
    130131               END DO 
    131132            ENDIF  
     
    220221               !  
    221222               ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) ) 
     223               ! 
     224               ! Set map structure 
     225               ibmap_ptr(1)%ptr => idx_bdy(ib_bdy)%nbmap(:,1) 
     226               ibmap_ptr(1)%ll_unstruc = ln_coords_file(ib_bdy) 
     227               ibmap_ptr(2)%ptr => idx_bdy(ib_bdy)%nbmap(:,2) 
     228               ibmap_ptr(2)%ll_unstruc = ln_coords_file(ib_bdy) 
     229               ibmap_ptr(3)%ptr => idx_bdy(ib_bdy)%nbmap(:,3) 
     230               ibmap_ptr(3)%ll_unstruc = ln_coords_file(ib_bdy) 
    222231 
    223232               ! Open files and read in tidal forcing data 
     
    228237                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc' 
    229238                  CALL iom_open( clfile, inum ) 
    230                   CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
     239                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1,  ibmap_ptr(1) ) 
    231240                  td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1) 
    232                   CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
     241                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1,  ibmap_ptr(1) ) 
    233242                  td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1) 
    234243                  CALL iom_close( inum ) 
     
    236245                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc' 
    237246                  CALL iom_open( clfile, inum ) 
    238                   CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
     247                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, ibmap_ptr(2) ) 
    239248                  td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1) 
    240                   CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
     249                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, ibmap_ptr(2) ) 
    241250                  td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1) 
    242251                  CALL iom_close( inum ) 
     
    244253                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc' 
    245254                  CALL iom_open( clfile, inum ) 
    246                   CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
     255                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, ibmap_ptr(3) ) 
    247256                  td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1) 
    248                   CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
     257                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, ibmap_ptr(3) ) 
    249258                  td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1) 
    250259                  CALL iom_close( inum ) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/C1D/c1d.F90

    • Property svn:keywords set to Id
    r4792 r5620  
    3131   !!---------------------------------------------------------------------- 
    3232   !! NEMO/C1D 3.3 , NEMO Consortium (2010) 
    33    !! $Id: c1d.F90 2382 2010-11-13 13:08:12Z gm $  
     33   !! $Id$  
    3434   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    3535   !!====================================================================== 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/C1D/domc1d.F90

    • Property svn:keywords set to Id
    r4792 r5620  
    2626   !!---------------------------------------------------------------------- 
    2727   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    28    !! $Id: domc1d.F90 3851 2013-04-30 10:30:51Z hadcv $  
     28   !! $Id$  
    2929   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    3030   !!---------------------------------------------------------------------- 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/C1D/dtauvd.F90

    • Property svn:keywords set to Id
    r4792 r5620  
    3535   !!---------------------------------------------------------------------- 
    3636   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    37    !! $Id: dtauvd.F90 2392 2010-11-15 21:20:05Z gm $  
     37   !! $Id$  
    3838   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    3939   !!---------------------------------------------------------------------- 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/C1D/dyncor_c1d.F90

    • Property svn:keywords set to Id
    r2409 r5620  
    3030   !!---------------------------------------------------------------------- 
    3131   !! NEMO/C1D 3.3 , NEMO Consortium (2010) 
    32    !! $Id: dyncor_c1d.F90 2382 2010-11-13 13:08:12Z gm $  
     32   !! $Id$  
    3333   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    3434   !!---------------------------------------------------------------------- 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/C1D/dyndmp.F90

    • Property svn:keywords set to Id
    r5038 r5620  
    33   !!                       ***  MODULE  dyndmp  *** 
    44   !! Ocean dynamics: internal restoring trend on momentum (U and V current) 
     5   !!                 This should only be used for C1D case in current form 
    56   !!====================================================================== 
    67   !! History :  3.5  ! 2013-08  (D. Calvert)  Original code 
     8   !!            3.6  ! 2014-08  (T. Graham) Modified to use netcdf file of 
     9   !!                             restoration coefficients supplied to tradmp 
    710   !!---------------------------------------------------------------------- 
    811 
     
    2528   USE wrk_nemo       ! Memory allocation 
    2629   USE timing         ! Timing 
     30   USE iom            ! I/O manager 
    2731 
    2832   IMPLICIT NONE 
     
    4347   !!---------------------------------------------------------------------- 
    4448   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    45    !! $Id: dyndmp.F90 3294 2012-01-28 16:44:18Z rblod $  
     49   !! $Id$  
    4650   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4751   !!---------------------------------------------------------------------- 
     
    7377      NAMELIST/namc1d_dyndmp/ ln_dyndmp 
    7478      INTEGER :: ios 
     79      INTEGER :: imask 
    7580      !!---------------------------------------------------------------------- 
    7681 
     
    9196         WRITE(numout,*) '      add a damping term or not       ln_dyndmp = ', ln_dyndmp 
    9297         WRITE(numout,*) '   Namelist namtra_dmp    : Set damping parameters' 
    93          WRITE(numout,*) '      horizontal damping option       nn_hdmp   = ', nn_hdmp 
    94          WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp, '(non-C1D zoom: forced to 0)' 
    95          WRITE(numout,*) '      surface time scale (days)       rn_surf   = ', rn_surf 
    96          WRITE(numout,*) '      bottom time scale (days)        rn_bot    = ', rn_bot 
    97          WRITE(numout,*) '      depth of transition (meters)    rn_dep    = ', rn_dep 
    98          WRITE(numout,*) '      create a damping.coeff file     nn_file   = ', nn_file 
     98         WRITE(numout,*) '      Apply relaxation   or not       ln_tradmp = ', ln_tradmp 
     99         WRITE(numout,*) '      mixed layer damping option      nn_zdmp   = ', nn_zdmp 
     100         WRITE(numout,*) '      Damping file name               cn_resto  = ', cn_resto 
    99101         WRITE(numout,*) 
    100102      ENDIF 
     
    104106         IF( dyn_dmp_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dyn_dmp_init: unable to allocate arrays' ) 
    105107         ! 
    106 #if ! defined key_c1d 
    107          SELECT CASE ( nn_hdmp )             !==   control print of horizontal option   ==! 
    108          CASE (  -1  )   ;   IF(lwp) WRITE(numout,*) '   momentum damping in the Med & Red seas only' 
    109          CASE ( 1:90 )   ;   IF(lwp) WRITE(numout,*) '   momentum damping poleward of', nn_hdmp, ' degrees' 
    110          CASE DEFAULT 
    111             WRITE(ctmp1,*) '          bad flag value for nn_hdmp = ', nn_hdmp 
    112             CALL ctl_stop(ctmp1) 
    113          END SELECT 
    114          ! 
    115 #endif 
    116108         SELECT CASE ( nn_zdmp )             !==   control print of vertical option   ==! 
    117109         CASE ( 0    )   ;   IF(lwp) WRITE(numout,*) '   momentum damping throughout the water column' 
     
    130122         utrdmp(:,:,:) = 0._wp               ! internal damping trends 
    131123         vtrdmp(:,:,:) = 0._wp 
    132          !                                   !==   Damping coefficients calculation:                           ==! 
    133          !                                   !==   use tradmp.F90 subroutines dtacof, dtacof_zoom and cofdis   ==! 
    134          !                                   !!!   NOTE: these need to be altered for use in this module if  
    135          !                                   !!!       they are to be used outside the C1D context  
    136          !                                   !!!       (use of U,V grid variables) 
    137          IF( lzoom .AND. .NOT. lk_c1d ) THEN   ;   CALL dtacof_zoom( resto_uv ) 
    138          ELSE               ;   CALL dtacof( nn_hdmp, rn_surf, rn_bot, rn_dep, nn_file, 'DYN', resto_uv ) 
    139          ENDIF 
    140          ! 
     124         ! 
     125         !Read in mask from file 
     126         CALL iom_open ( cn_resto, imask) 
     127         CALL iom_get  ( imask, jpdom_autoglo, 'resto', resto) 
     128         CALL iom_close( imask ) 
    141129      ENDIF 
    142130      ! 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/C1D/dynnxt_c1d.F90

    • Property svn:keywords set to Id
    r2409 r5620  
    2525   !!---------------------------------------------------------------------- 
    2626   !! NEMO/C1D 3.3 , NEMO Consortium (2010) 
    27    !! $Id: dynnxt_c1d.F90 2382 2010-11-13 13:08:12Z gm $  
     27   !! $Id$  
    2828   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    2929   !!---------------------------------------------------------------------- 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/C1D/step_c1d.F90

    • Property svn:keywords set to Id
    r5038 r5620  
    3232   !!---------------------------------------------------------------------- 
    3333   !! NEMO/C1D 3.3 , NEMO Consortium (2010) 
    34    !! $Id: step_c1d.F90 2382 2010-11-13 13:08:12Z gm $ 
     34   !! $Id$ 
    3535   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    3636   !!---------------------------------------------------------------------- 
     
    7272      ! Ocean physics update                (ua, va, ta, sa used as workspace) 
    7373      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     74                         CALL eos_rab( tsb, rab_b )   ! before local thermal/haline expension ratio at T-points 
     75                         CALL eos_rab( tsn, rab_n )   ! now    local thermal/haline expension ratio at T-points 
    7476                         CALL bn2( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency 
    7577                         CALL bn2( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency 
     
    132134                        CALL tra_nxt( kstp )       ! tracer fields at next time step 
    133135 
     136 
     137 
    134138      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    135139      ! Dynamics                                    (ta, sa used as workspace) 
     
    150154      IF( lrst_oce       )   CALL rst_write( kstp )        ! write output ocean restart file 
    151155      ! 
     156#if defined key_iomput 
     157      IF( kstp == nitend .OR. indic < 0 )   CALL xios_context_finalize()   ! needed for XIOS 
     158      ! 
     159#endif 
    152160   END SUBROUTINE stp_c1d 
    153161 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/CRS/crs.F90

    • Property svn:keywords set to Id
    r4064 r5620  
    164164 
    165165 
     166   !! $Id$ 
    166167CONTAINS 
    167168    
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/CRS/crsdom.F90

    • Property svn:keywords set to Id
    r4314 r5620  
    5757#  include "domzgr_substitute.h90" 
    5858    
     59   !! $Id$ 
    5960CONTAINS 
    6061 
     
    18821883      CALL crs_lbc_lnk( p_surf_crs_msk, cd_type, 1.0, pval=1.0 ) 
    18831884 
    1884       CALL wrk_dealloc( jpi, jpj, jpk, zsurf, zsurfmsk ) 
     1885      CALL wrk_dealloc( jpi, jpj, jpk, zsurfmsk, zsurf ) 
    18851886 
    18861887   END SUBROUTINE crs_dom_sfc 
     
    22742275      ENDDO 
    22752276      
    2276       CALL wrk_alloc( jpi_crs, jpj_crs, zmbk ) 
    2277  
    22782277      zmbk(:,:) = 0.0 
    22792278      zmbk(:,:) = REAL( mbathy_crs(:,:), wp ) ;   CALL crs_lbc_lnk(zmbk,'T',1.0)   ;   mbathy_crs(:,:) = INT( zmbk(:,:) ) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/CRS/crsdomwri.F90

    • Property svn:keywords set to Id
    r4294 r5620  
    3333   PUBLIC crs_dom_wri        ! routine called by crsini.F90 
    3434 
     35   !! $Id$ 
    3536CONTAINS 
    3637 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/CRS/crsfld.F90

    • Property svn:keywords set to Id
    r4149 r5620  
    3838   !!---------------------------------------------------------------------- 
    3939   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    40    !! $Id $ 
     40   !! $Id$ 
    4141   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4242   !!---------------------------------------------------------------------- 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/CRS/crsini.F90

    • Property svn:keywords set to Id
    r4792 r5620  
    2929#  include "domzgr_substitute.h90" 
    3030 
     31   !! $Id$ 
    3132CONTAINS 
    3233    
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/CRS/crslbclnk.F90

    • Property svn:keywords set to Id
    r4015 r5620  
    2222   PUBLIC crs_lbc_lnk 
    2323    
     24   !! $Id$ 
    2425CONTAINS 
    2526 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

    r5038 r5620  
    2121   USE timing         ! preformance summary 
    2222   USE wrk_nemo       ! working arrays 
     23   USE fldread        ! type FLD_N 
     24   USE phycst         ! physical constant 
     25   USE in_out_manager  ! I/O manager 
    2326 
    2427   IMPLICIT NONE 
     
    8386      CALL wrk_alloc( jpi , jpj , jpk , jpts , ztsn                 ) 
    8487 
    85       CALL iom_put( 'cellthc', fse3t(:,:,:) ) 
    86  
    8788      zarea_ssh(:,:) = area(:,:) * sshn(:,:) 
    8889 
     
    105106      END DO 
    106107      IF( .NOT.lk_vvl ) THEN 
    107          DO ji=1,jpi 
    108             DO jj=1,jpj 
    109                zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
    110             END DO 
    111          END DO 
     108         IF ( ln_isfcav ) THEN 
     109            DO ji=1,jpi 
     110               DO jj=1,jpj 
     111                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
     112               END DO 
     113            END DO 
     114         ELSE 
     115            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     116         END IF 
    112117      END IF 
    113118      !                                          
     
    127132      END DO 
    128133      IF( .NOT.lk_vvl ) THEN 
    129          DO ji=1,jpi 
    130             DO jj=1,jpj 
    131                zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
    132             END DO 
    133          END DO 
     134         IF ( ln_isfcav ) THEN 
     135            DO ji=1,jpi 
     136               DO jj=1,jpj 
     137                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
     138               END DO 
     139            END DO 
     140         ELSE 
     141            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     142         END IF 
    134143      END IF 
    135144      !     
     
    157166      END DO 
    158167      IF( .NOT.lk_vvl ) THEN 
    159          DO ji=1,jpi 
    160             DO jj=1,jpj 
    161                ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem)  
    162                zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal)  
    163             END DO 
    164          END DO 
     168         IF ( ln_isfcav ) THEN 
     169            DO ji=1,jpi 
     170               DO jj=1,jpj 
     171                  ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem)  
     172                  zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal)  
     173               END DO 
     174            END DO 
     175         ELSE 
     176            ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) ) 
     177            zsal  = zsal  + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) ) 
     178         END IF 
    165179      ENDIF 
    166180      IF( lk_mpp ) THEN   
     
    197211      REAL(wp) ::   zztmp   
    198212      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zsaldta   ! Jan/Dec levitus salinity 
     213      ! reading initial file 
     214      LOGICAL  ::   ln_tsd_init      !: T & S data flag 
     215      LOGICAL  ::   ln_tsd_tradmp    !: internal damping toward input data flag 
     216      CHARACTER(len=100)            ::   cn_dir 
     217      TYPE(FLD_N)                   ::  sn_tem,sn_sal 
     218      INTEGER  ::   ios=0 
     219 
     220      NAMELIST/namtsd/ ln_tsd_init,ln_tsd_tradmp,cn_dir,sn_tem,sn_sal 
     221      ! 
     222 
     223      REWIND( numnam_ref )              ! Namelist namtsd in reference namelist : 
     224      READ  ( numnam_ref, namtsd, IOSTAT = ios, ERR = 901) 
     225901   IF( ios /= 0 ) CALL ctl_nam ( ios , ' namtsd in reference namelist for dia_ar5', lwp ) 
     226      REWIND( numnam_cfg )              ! Namelist namtsd in configuration namelist : Parameters of the run 
     227      READ  ( numnam_cfg, namtsd, IOSTAT = ios, ERR = 902 ) 
     228902   IF( ios /= 0 ) CALL ctl_nam ( ios , ' namtsd in configuration namelist for dia_ar5', lwp ) 
     229      IF(lwm) WRITE ( numond, namtsd ) 
     230      ! 
    199231      !!---------------------------------------------------------------------- 
    200232      ! 
     
    216248      END DO 
    217249      IF( lk_mpp )   CALL mpp_sum( vol0 ) 
    218        
    219       CALL iom_open ( 'data_1m_salinity_nomask', inum ) 
    220       CALL iom_get  ( inum, jpdom_data, 'vosaline', zsaldta(:,:,:,1), 1  ) 
    221       CALL iom_get  ( inum, jpdom_data, 'vosaline', zsaldta(:,:,:,2), 12 ) 
     250 
     251      CALL iom_open ( TRIM( cn_dir )//TRIM(sn_sal%clname), inum ) 
     252      CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,1), 1  ) 
     253      CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,2), 12 ) 
    222254      CALL iom_close( inum ) 
    223255      sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )         
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90

    • Property svn:keywords set to Id
    r4792 r5620  
    4242#endif 
    4343#if defined key_lim3 
    44   USE par_ice 
    4544  USE ice 
    4645#endif 
     
    113112  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d   
    114113 
     114   !! $Id$ 
    115115CONTAINS 
    116116 
     
    176176 
    177177     !open output file 
    178      IF( lwp ) THEN 
     178     IF( lwm ) THEN 
    179179        CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
    180180        CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. ) 
     
    283283           DO jsec=1,nb_sec 
    284284 
    285               IF( lwp )CALL dia_dct_wri(kt,jsec,secs(jsec)) 
     285              IF( lwm )CALL dia_dct_wri(kt,jsec,secs(jsec)) 
    286286             
    287287              !nullify transports values after writing 
     
    12981298   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag 
    12991299   PUBLIC  
     1300   !! $Id$ 
    13001301CONTAINS 
    13011302 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DIA/diafwb.F90

    r5038 r5620  
    5151      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    5252      !! 
    53       INTEGER :: inum             ! temporary logical unit 
    54       INTEGER :: ji, jj, jk, jt   ! dummy loop indices 
    55       INTEGER :: ii0, ii1, ij0, ij1 
    56       REAL(wp) ::   zarea, zvol, zwei 
    57       REAL(wp) ::  ztemi(4), ztemo(4), zsali(4), zsalo(4), zflxi(4), zflxo(4) 
    58       REAL(wp) ::  zt, zs, zu   
    59       REAL(wp) ::  zsm0, zfwfnew 
     53      INTEGER  :: inum             ! temporary logical unit 
     54      INTEGER  :: ji, jj, jk, jt   ! dummy loop indices 
     55      INTEGER  :: ii0, ii1, ij0, ij1 
     56      INTEGER  :: isrow         ! index for ORCA1 starting row 
     57      REAL(wp) :: zarea, zvol, zwei 
     58      REAL(wp) :: ztemi(4), ztemo(4), zsali(4), zsalo(4), zflxi(4), zflxo(4) 
     59      REAL(wp) :: zt, zs, zu   
     60      REAL(wp) :: zsm0, zfwfnew 
    6061      IF( cp_cfg == "orca" .AND. jp_cfg == 1 .OR. jp_cfg == 2 .OR. jp_cfg == 4 ) THEN 
    6162      !!---------------------------------------------------------------------- 
     
    165166         CASE ( 1 )                                  !  ORCA_R1 configurations 
    166167            !                                        ! ======================= 
    167             ii0 = 283   ;   ii1 = 283 
    168             ij0 = 200   ;   ij1 = 200 
     168            ! This dirty section will be suppressed by simplification process: 
     169            ! all this will come back in input files 
     170            ! Currently these hard-wired indices relate to configuration with 
     171            ! extend grid (jpjglo=332) 
     172            isrow = 332 - jpjglo 
     173            ! 
     174            ii0 = 283           ;   ii1 = 283 
     175            ij0 = 241 - isrow   ;   ij1 = 241 - isrow 
    169176            !                                        ! ======================= 
    170177         CASE DEFAULT                                !    ORCA R05 or R025 
     
    212219         CASE ( 1 )                                  !  ORCA_R1 configurations 
    213220            !                                        ! ======================= 
    214             ii0 = 282   ;   ii1 = 282 
    215             ij0 = 200   ;   ij1 = 200 
     221            ! This dirty section will be suppressed by simplification process: 
     222            ! all this will come back in input files 
     223            ! Currently these hard-wired indices relate to configuration with 
     224            ! extend grid (jpjglo=332) 
     225            isrow = 332 - jpjglo 
     226            ii0 = 282           ;   ii1 = 282 
     227            ij0 = 240 - isrow   ;   ij1 = 240 - isrow 
    216228            !                                        ! ======================= 
    217229         CASE DEFAULT                                !    ORCA R05 or R025 
     
    259271         CASE ( 1 )                                  !  ORCA_R1 configurations 
    260272            !                                        ! ======================= 
    261             ii0 = 331   ;   ii1 = 331 
    262             ij0 = 176   ;   ij1 = 176 
     273            ! This dirty section will be suppressed by simplification process: 
     274            ! all this will come back in input files 
     275            ! Currently these hard-wired indices relate to configuration with 
     276            ! extend grid (jpjglo=332) 
     277            isrow = 332 - jpjglo 
     278            ii0 = 331           ;   ii1 = 331 
     279            ij0 = 215 - isrow   ;   ij1 = 215 - isrow 
    263280            !                                        ! ======================= 
    264281         CASE DEFAULT                                !    ORCA R05 or R025 
     
    306323         CASE ( 1 )                                  !  ORCA_R1 configurations 
    307324            !                                        ! ======================= 
    308             ii0 = 297   ;   ii1 = 297  
    309             ij0 = 230   ;   ij1 = 230 
     325            ! This dirty section will be suppressed by simplification process: 
     326            ! all this will come back in input files 
     327            ! Currently these hard-wired indices relate to configuration with 
     328            ! extend grid (jpjglo=332) 
     329            isrow = 332 - jpjglo 
     330            ii0 = 297           ;   ii1 = 297 
     331            ij0 = 269 - isrow   ;   ij1 = 269 - isrow 
    310332            !                                        ! ======================= 
    311333         CASE DEFAULT                                !    ORCA R05 or R025 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90

    • Property svn:keywords set to Id
    r5038 r5620  
    6060   !!---------------------------------------------------------------------- 
    6161   !! NEMO/OPA 3.5 , NEMO Consortium (2013) 
    62    !! $Id:$ 
     62   !! $Id$ 
    6363   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    6464   !!---------------------------------------------------------------------- 
     
    196196                  DO ji = 1,jpi 
    197197                     ! Elevation 
    198                      ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)           *tmask_i(ji,jj)         
    199 #if defined key_dynspg_ts 
    200                      ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*hur(ji,jj)*umask_i(ji,jj) 
    201                      ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*hvr(ji,jj)*vmask_i(ji,jj) 
    202 #endif 
     198                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*tmask_i(ji,jj)         
     199                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*umask_i(ji,jj) 
     200                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*vmask_i(ji,jj) 
    203201                  END DO 
    204202               END DO 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90

    r5038 r5620  
    9696      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) )                               ! heat fluxes 
    9797      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) )                               ! salt fluxes 
    98       ! Add runoff heat & salt input 
     98      ! Add runoff    heat & salt input 
    9999      IF( ln_rnf    )   z_frc_trd_t = z_frc_trd_t + glob_sum( rnf_tsc(:,:,jp_tem) * surf(:,:) ) 
    100100      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * surf(:,:) ) 
    101       ! Add geothermal ice shelf 
     101      ! Add ice shelf heat & salt input 
    102102      IF( nn_isf .GE. 1 )  THEN 
    103103          z_frc_trd_t = z_frc_trd_t & 
     
    112112      ! 
    113113      IF( .NOT. lk_vvl ) THEN 
    114          z2d0=0.0_wp ; z2d1=0.0_wp 
    115          DO ji=1,jpi 
    116             DO jj=1,jpj 
    117               z2d0(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_tem) 
    118               z2d1(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_sal) 
     114         IF ( ln_isfcav ) THEN 
     115            DO ji=1,jpi 
     116               DO jj=1,jpj 
     117                  z2d0(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_tem) 
     118                  z2d1(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_sal) 
     119               ENDDO 
    119120            ENDDO 
    120          ENDDO 
     121         ELSE 
     122            z2d0(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem) 
     123            z2d1(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal) 
     124         END IF 
    121125         z_wn_trd_t = - glob_sum( z2d0 )  
    122126         z_wn_trd_s = - glob_sum( z2d1 ) 
     
    144148      ! heat & salt content variation (associated with ssh) 
    145149      IF( .NOT. lk_vvl ) THEN 
    146          z2d0 = 0._wp   ;   z2d1 = 0._wp 
    147          DO ji = 1, jpi 
    148             DO jj = 1, jpj 
    149               z2d0(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj) - ssh_hc_loc_ini(ji,jj) )  
    150               z2d1(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj) - ssh_sc_loc_ini(ji,jj) )  
     150         IF ( ln_isfcav ) THEN 
     151            DO ji = 1, jpi 
     152               DO jj = 1, jpj 
     153                  z2d0(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj) - ssh_hc_loc_ini(ji,jj) )  
     154                  z2d1(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj) - ssh_sc_loc_ini(ji,jj) )  
     155               END DO 
    151156            END DO 
    152          END DO 
     157         ELSE 
     158            z2d0(:,:) = surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) )  
     159            z2d1(:,:) = surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) )  
     160         END IF 
    153161         z_ssh_hc = glob_sum( z2d0 )  
    154162         z_ssh_sc = glob_sum( z2d1 )  
     
    277285          frc_s = 0._wp                                           ! salt content   -    -   -    -         
    278286          IF( .NOT. lk_vvl ) THEN 
    279              DO ji=1,jpi 
    280                 DO jj=1,jpj 
    281                    ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj)   ! initial heat content in ssh 
    282                    ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj)   ! initial salt content in ssh 
     287             IF ( ln_isfcav ) THEN 
     288                DO ji=1,jpi 
     289                   DO jj=1,jpj 
     290                      ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj)   ! initial heat content in ssh 
     291                      ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj)   ! initial salt content in ssh 
     292                   ENDDO 
    283293                ENDDO 
    284              ENDDO 
     294             ELSE 
     295                ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh 
     296                ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh 
     297             END IF 
    285298             frc_wn_t = 0._wp                                       ! initial heat content misfit due to free surface 
    286299             frc_wn_s = 0._wp                                       ! initial salt content misfit due to free surface 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90

    r4292 r5620  
    245245      CALL iom_put( "mldr10_3", zrho10_3     )   ! MLD delta rho(10m) = 0.03 
    246246      CALL iom_put( "pycndep" , zpycn        )   ! MLD delta rho equi. delta T(10m) = 0.2 
    247       CALL iom_put( "BLT"     , ztm2 - zpycn )   ! Barrier Layer Thickness 
    248247      CALL iom_put( "tinv"    , ztinv        )   ! max. temp. inv. (t10 ref)  
    249248      CALL iom_put( "depti"   , zdepinv      )   ! depth of max. temp. inv. (t10 ref)  
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90

    r5038 r5620  
    88   !!            3.2  ! 2010-03  (O. Marti, S. Flavoni) Add fields 
    99   !!            3.3  ! 2010-10  (G. Madec)  dynamical allocation 
     10   !!            3.6  ! 2014-12  (C. Ethe) use of IOM 
    1011   !!---------------------------------------------------------------------- 
    1112 
     
    1314   !!   dia_ptr      : Poleward Transport Diagnostics module 
    1415   !!   dia_ptr_init : Initialization, namelist read 
    15    !!   dia_ptr_wri  : Output of poleward fluxes 
    16    !!   ptr_vjk      : "zonal" sum computation of a "meridional" flux array 
    17    !!   ptr_tjk      : "zonal" mean computation of a tracer field 
    18    !!   ptr_vj       : "zonal" and vertical sum computation of a "meridional" flux array 
    19    !!                   (Generic interface to ptr_vj_3d, ptr_vj_2d) 
     16   !!   ptr_sjk      : "zonal" mean computation of a field - tracer or flux array 
     17   !!   ptr_sj       : "zonal" and vertical sum computation of a "meridional" flux array 
     18   !!                   (Generic interface to ptr_sj_3d, ptr_sj_2d) 
    2019   !!---------------------------------------------------------------------- 
    2120   USE oce              ! ocean dynamics and active tracers 
    2221   USE dom_oce          ! ocean space and time domain 
    2322   USE phycst           ! physical constants 
    24    USE ldftra_oce       ! ocean active tracers: lateral physics 
    25    USE dianam           ! 
     23   ! 
    2624   USE iom              ! IOM library 
    27    USE ioipsl           ! IO-IPSL library 
    2825   USE in_out_manager   ! I/O manager 
    2926   USE lib_mpp          ! MPP library 
    30    USE lbclnk           ! lateral boundary condition - processor exchanges 
    3127   USE timing           ! preformance summary 
    32    USE wrk_nemo         ! working arrays 
    3328 
    3429   IMPLICIT NONE 
    3530   PRIVATE 
    3631 
    37    INTERFACE ptr_vj 
    38       MODULE PROCEDURE ptr_vj_3d, ptr_vj_2d 
     32   INTERFACE ptr_sj 
     33      MODULE PROCEDURE ptr_sj_3d, ptr_sj_2d 
    3934   END INTERFACE 
    4035 
    41    PUBLIC   dia_ptr_init   ! call in opa module 
     36   PUBLIC   ptr_sj         ! call by tra_ldf & tra_adv routines 
     37   PUBLIC   ptr_sjk        !  
     38   PUBLIC   dia_ptr_init   ! call in step module 
    4239   PUBLIC   dia_ptr        ! call in step module 
    43    PUBLIC   ptr_vj         ! call by tra_ldf & tra_adv routines 
    44    PUBLIC   ptr_vjk        ! call by tra_ldf & tra_adv routines 
    4540 
    4641   !                                  !!** namelist  namptr  ** 
    47    LOGICAL , PUBLIC ::   ln_diaptr     !: Poleward transport flag (T) or not (F) 
    48    LOGICAL , PUBLIC ::   ln_subbas     !: Atlantic/Pacific/Indian basins calculation 
    49    LOGICAL , PUBLIC ::   ln_diaznl     !: Add zonal means and meridional stream functions 
    50    LOGICAL , PUBLIC ::   ln_ptrcomp    !: Add decomposition : overturning (and gyre, soon ...) 
    51    INTEGER , PUBLIC ::   nn_fptr       !: frequency of ptr computation  [time step] 
    52    INTEGER , PUBLIC ::   nn_fwri       !: frequency of ptr outputs      [time step] 
    53  
    54    REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   htr_adv, htr_ldf, htr_ove   !: Heat TRansports (adv, diff, overturn.) 
    55    REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   str_adv, str_ldf, str_ove   !: Salt TRansports (adv, diff, overturn.) 
     42   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   htr_adv, htr_ldf   !: Heat TRansports (adv, diff, overturn.) 
     43   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:) ::   str_adv, str_ldf   !: Salt TRansports (adv, diff, overturn.) 
    5644    
    57    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   btmsk                  ! T-point basin interior masks 
    58    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   btm30                  ! mask out Southern Ocean (=0 south of 30°S) 
    59    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htr  , str             ! adv heat and salt transports (approx) 
    60    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn_jk, sn_jk , v_msf   ! i-mean T and S, j-Stream-Function 
    61    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sjk  , r1_sjk          ! i-mean i-k-surface and its inverse         
    62    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htr_eiv, str_eiv       ! bolus adv heat ans salt transports ('key_diaeiv') 
    63    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_msf_eiv              ! bolus j-streamfuction              ('key_diaeiv') 
    64  
    65  
    66    INTEGER ::   niter       ! 
    67    INTEGER ::   nidom_ptr   ! 
    68    INTEGER ::   numptr      ! logical unit for Poleward TRansports 
    69    INTEGER ::   nptr        ! = 1 (ln_subbas=F) or = 5 (glo, atl, pac, ind, ipc) (ln_subbas=T)  
     45 
     46   LOGICAL, PUBLIC ::   ln_diaptr   !  Poleward transport flag (T) or not (F) 
     47   LOGICAL, PUBLIC ::   ln_subbas   !  Atlantic/Pacific/Indian basins calculation 
     48   INTEGER         ::   nptr        ! = 1 (l_subbas=F) or = 5 (glo, atl, pac, ind, ipc) (l_subbas=T)  
    7049 
    7150   REAL(wp) ::   rc_sv    = 1.e-6_wp   ! conversion from m3/s to Sverdrup 
     
    7352   REAL(wp) ::   rc_ggram = 1.e-6_wp   ! conversion from g    to Pg 
    7453 
    75    REAL(wp), TARGET, DIMENSION(:),   ALLOCATABLE, SAVE :: p_fval1d 
    76    REAL(wp), TARGET, DIMENSION(:,:), ALLOCATABLE, SAVE :: p_fval2d 
    77  
    78    !! Integer, 1D workspace arrays. Not common enough to be implemented in  
    79    !! wrk_nemo module. 
    80    INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndex  , ndex_atl     , ndex_pac     , ndex_ind     , ndex_ipc 
    81    INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) ::         ndex_atl_30  , ndex_pac_30  , ndex_ind_30  , ndex_ipc_30 
    82    INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndex_h, ndex_h_atl_30, ndex_h_pac_30, ndex_h_ind_30, ndex_h_ipc_30 
     54   CHARACTER(len=3), ALLOCATABLE, SAVE, DIMENSION(:)     :: clsubb 
     55   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: btmsk   ! T-point basin interior masks 
     56   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   :: btm30   ! mask out Southern Ocean (=0 south of 30°S) 
     57 
     58   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:)     :: p_fval1d 
     59   REAL(wp), TARGET, ALLOCATABLE, SAVE, DIMENSION(:,:)   :: p_fval2d 
     60 
    8361 
    8462   !! * Substitutions 
     
    9270CONTAINS 
    9371 
    94    FUNCTION dia_ptr_alloc() 
    95       !!---------------------------------------------------------------------- 
    96       !!                    ***  ROUTINE dia_ptr_alloc  *** 
    97       !!---------------------------------------------------------------------- 
    98       INTEGER               ::   dia_ptr_alloc   ! return value 
    99       INTEGER, DIMENSION(6) ::   ierr 
    100       !!---------------------------------------------------------------------- 
    101       ierr(:) = 0 
    102       ! 
    103       ALLOCATE( btmsk(jpi,jpj,nptr) ,           & 
    104          &      htr_adv(jpj) , str_adv(jpj) ,   & 
    105          &      htr_ldf(jpj) , str_ldf(jpj) ,   & 
    106          &      htr_ove(jpj) , str_ove(jpj),    & 
    107          &      htr(jpj,nptr) , str(jpj,nptr) , & 
    108          &      tn_jk(jpj,jpk,nptr) , sn_jk (jpj,jpk,nptr) , v_msf(jpj,jpk,nptr) , & 
    109          &      sjk  (jpj,jpk,nptr) , r1_sjk(jpj,jpk,nptr) , STAT=ierr(1)  ) 
    110          ! 
    111 #if defined key_diaeiv 
    112       ALLOCATE( htr_eiv(jpj,nptr) , str_eiv(jpj,nptr) , & 
    113          &      v_msf_eiv(jpj,jpk,nptr) , STAT=ierr(2) ) 
    114 #endif 
    115       ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(3)) 
    116       ! 
    117       ALLOCATE(ndex(jpj*jpk),        ndex_atl(jpj*jpk), ndex_pac(jpj*jpk), & 
    118          &     ndex_ind(jpj*jpk),    ndex_ipc(jpj*jpk),                    & 
    119          &     ndex_atl_30(jpj*jpk), ndex_pac_30(jpj*jpk), Stat=ierr(4)) 
    120  
    121       ALLOCATE(ndex_ind_30(jpj*jpk), ndex_ipc_30(jpj*jpk),                   & 
    122          &     ndex_h(jpj),          ndex_h_atl_30(jpj), ndex_h_pac_30(jpj), & 
    123          &     ndex_h_ind_30(jpj),   ndex_h_ipc_30(jpj), Stat=ierr(5) ) 
    124          ! 
    125      ALLOCATE( btm30(jpi,jpj) , STAT=ierr(6)  ) 
    126          ! 
    127       dia_ptr_alloc = MAXVAL( ierr ) 
    128       IF(lk_mpp)   CALL mpp_sum( dia_ptr_alloc ) 
    129       ! 
    130    END FUNCTION dia_ptr_alloc 
    131  
    132  
    133    FUNCTION ptr_vj_3d( pva )   RESULT ( p_fval ) 
    134       !!---------------------------------------------------------------------- 
    135       !!                    ***  ROUTINE ptr_vj_3d  *** 
    136       !! 
    137       !! ** Purpose :   i-k sum computation of a j-flux array 
    138       !! 
    139       !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i). 
    140       !!              pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v) 
    141       !! 
    142       !! ** Action  : - p_fval: i-k-mean poleward flux of pva 
    143       !!---------------------------------------------------------------------- 
    144       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pva   ! mask flux array at V-point 
    145       !! 
    146       INTEGER                  ::   ji, jj, jk   ! dummy loop arguments 
    147       INTEGER                  ::   ijpj         ! ??? 
    148       REAL(wp), POINTER, DIMENSION(:) :: p_fval  ! function value 
    149       !!-------------------------------------------------------------------- 
    150       ! 
    151       p_fval => p_fval1d 
    152  
    153       ijpj = jpj 
    154       p_fval(:) = 0._wp 
    155       DO jk = 1, jpkm1 
    156          DO jj = 2, jpjm1 
    157             DO ji = fs_2, fs_jpim1   ! Vector opt. 
    158                p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj)  
    159             END DO 
    160          END DO 
    161       END DO 
    162 #if defined key_mpp_mpi 
    163       IF(lk_mpp)   CALL mpp_sum( p_fval, ijpj, ncomm_znl) 
    164 #endif 
    165       ! 
    166    END FUNCTION ptr_vj_3d 
    167  
    168  
    169    FUNCTION ptr_vj_2d( pva )   RESULT ( p_fval ) 
    170       !!---------------------------------------------------------------------- 
    171       !!                    ***  ROUTINE ptr_vj_2d  *** 
    172       !! 
    173       !! ** Purpose :   "zonal" and vertical sum computation of a i-flux array 
    174       !! 
    175       !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i). 
    176       !!      pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v) 
    177       !! 
    178       !! ** Action  : - p_fval: i-k-mean poleward flux of pva 
    179       !!---------------------------------------------------------------------- 
    180       IMPLICIT none 
    181       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj) ::   pva   ! mask flux array at V-point 
    182       !! 
    183       INTEGER                  ::   ji,jj       ! dummy loop arguments 
    184       INTEGER                  ::   ijpj        ! ??? 
    185       REAL(wp), POINTER, DIMENSION(:) :: p_fval ! function value 
    186       !!-------------------------------------------------------------------- 
    187       !  
    188       p_fval => p_fval1d 
    189  
    190       ijpj = jpj 
    191       p_fval(:) = 0._wp 
    192       DO jj = 2, jpjm1 
    193          DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ? 
    194             p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj) 
    195          END DO 
    196       END DO 
    197 #if defined key_mpp_mpi 
    198       CALL mpp_sum( p_fval, ijpj, ncomm_znl ) 
    199 #endif 
    200       !  
    201    END FUNCTION ptr_vj_2d 
    202  
    203  
    204    FUNCTION ptr_vjk( pva, pmsk )   RESULT ( p_fval ) 
    205       !!---------------------------------------------------------------------- 
    206       !!                    ***  ROUTINE ptr_vjk  *** 
    207       !! 
    208       !! ** Purpose :   i-sum computation of a j-velocity array 
    209       !! 
    210       !! ** Method  : - i-sum of pva using the interior 2D vmask (vmask_i). 
    211       !!              pva is supposed to be a masked flux (i.e. * vmask) 
    212       !! 
    213       !! ** Action  : - p_fval: i-mean poleward flux of pva 
    214       !!---------------------------------------------------------------------- 
    215       !! 
    216       IMPLICIT none 
    217       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk)           ::   pva    ! mask flux array at V-point 
    218       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)    , OPTIONAL ::   pmsk   ! Optional 2D basin mask 
    219       !! 
    220       INTEGER                           :: ji, jj, jk ! dummy loop arguments 
    221       REAL(wp), POINTER, DIMENSION(:,:) :: p_fval     ! return function value 
    222 #if defined key_mpp_mpi 
    223       INTEGER, DIMENSION(1) ::   ish 
    224       INTEGER, DIMENSION(2) ::   ish2 
    225       INTEGER               ::   ijpjjpk 
    226 #endif 
    227 #if defined key_mpp_mpi 
    228       REAL(wp), POINTER, DIMENSION(:) ::   zwork    ! mask flux array at V-point 
    229 #endif 
    230       !!-------------------------------------------------------------------- 
    231       ! 
    232 #if defined key_mpp_mpi 
    233       ijpjjpk = jpj*jpk 
    234       CALL wrk_alloc( jpj*jpk, zwork ) 
    235 #endif 
    236  
    237       p_fval => p_fval2d 
    238  
    239       p_fval(:,:) = 0._wp 
    240       ! 
    241       IF( PRESENT( pmsk ) ) THEN  
    242          DO jk = 1, jpkm1 
    243             DO jj = 2, jpjm1 
    244 !!gm here, use of tmask_i  ==> no need of loop over nldi, nlei.... 
    245                DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ? 
    246                   p_fval(jj,jk) = p_fval(jj,jk) + pva(ji,jj,jk) * e1v(ji,jj) * fse3v(ji,jj,jk) * pmsk(ji,jj) 
     72   SUBROUTINE dia_ptr( pvtr ) 
     73      !!---------------------------------------------------------------------- 
     74      !!                  ***  ROUTINE dia_ptr  *** 
     75      !!---------------------------------------------------------------------- 
     76      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport 
     77      ! 
     78      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     79      REAL(wp) ::   zv, zsfc               ! local scalar 
     80      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d   ! 2D workspace 
     81      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  z3d   ! 3D workspace 
     82      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zmask   ! 3D workspace 
     83      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::  zts   ! 3D workspace 
     84      CHARACTER( len = 10 )  :: cl1 
     85      !!---------------------------------------------------------------------- 
     86      ! 
     87      IF( nn_timing == 1 )   CALL timing_start('dia_ptr') 
     88 
     89      ! 
     90      IF( PRESENT( pvtr ) ) THEN 
     91         IF( iom_use("zomsfglo") ) THEN    ! effective MSF 
     92            z3d(1,:,:) = ptr_sjk( pvtr(:,:,:) )  ! zonal cumulative effective transport 
     93            DO jk = 2, jpkm1  
     94              z3d(1,:,jk) = z3d(1,:,jk-1) + z3d(1,:,jk)   ! effective j-Stream-Function (MSF) 
     95            END DO 
     96            DO ji = 1, jpi 
     97               z3d(ji,:,:) = z3d(1,:,:) 
     98            ENDDO 
     99            cl1 = TRIM('zomsf'//clsubb(1) ) 
     100            CALL iom_put( cl1, z3d * rc_sv ) 
     101            DO jn = 2, nptr                                    ! by sub-basins 
     102               z3d(1,:,:) =  ptr_sjk( pvtr(:,:,:), btmsk(:,:,jn)*btm30(:,:) )  
     103               DO jk = 2, jpkm1  
     104                  z3d(1,:,jk) = z3d(1,:,jk-1) + z3d(1,:,jk)    ! effective j-Stream-Function (MSF) 
    247105               END DO 
    248             END DO 
    249          END DO 
    250       ELSE  
    251          DO jk = 1, jpkm1 
    252             DO jj = 2, jpjm1 
    253                DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ? 
    254                   p_fval(jj,jk) = p_fval(jj,jk) + pva(ji,jj,jk) * e1v(ji,jj) * fse3v(ji,jj,jk) * tmask_i(ji,jj) 
    255                END DO 
    256             END DO 
    257          END DO 
    258       END IF 
    259       ! 
    260 #if defined key_mpp_mpi 
    261       ijpjjpk = jpj*jpk 
    262       ish(1) = ijpjjpk  ;   ish2(1) = jpj   ;   ish2(2) = jpk 
    263       zwork(1:ijpjjpk) = RESHAPE( p_fval, ish ) 
    264       CALL mpp_sum( zwork, ijpjjpk, ncomm_znl ) 
    265       p_fval(:,:) = RESHAPE( zwork, ish2 ) 
    266 #endif 
    267       ! 
    268 #if defined key_mpp_mpi 
    269       CALL wrk_dealloc( jpj*jpk, zwork ) 
    270 #endif 
    271       ! 
    272    END FUNCTION ptr_vjk 
    273  
    274  
    275    FUNCTION ptr_tjk( pta, pmsk )   RESULT ( p_fval ) 
    276       !!---------------------------------------------------------------------- 
    277       !!                    ***  ROUTINE ptr_tjk  *** 
    278       !! 
    279       !! ** Purpose :   i-sum computation of e1t*e3t * a tracer field 
    280       !! 
    281       !! ** Method  : - i-sum of mj(pta) using tmask 
    282       !! 
    283       !! ** Action  : - p_fval: i-sum of e1t*e3t*pta 
    284       !!---------------------------------------------------------------------- 
    285       !! 
    286       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pta    ! tracer flux array at T-point 
    287       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)     ::   pmsk   ! Optional 2D basin mask 
    288       !! 
    289       INTEGER                           :: ji, jj, jk   ! dummy loop arguments 
    290       REAL(wp), POINTER, DIMENSION(:,:) :: p_fval       ! return function value 
    291 #if defined key_mpp_mpi 
    292       INTEGER, DIMENSION(1) ::   ish 
    293       INTEGER, DIMENSION(2) ::   ish2 
    294       INTEGER               ::   ijpjjpk 
    295 #endif 
    296 #if defined key_mpp_mpi 
    297       REAL(wp), POINTER, DIMENSION(:) ::   zwork    ! mask flux array at V-point 
    298 #endif 
    299       !!--------------------------------------------------------------------  
    300       ! 
    301 #if defined key_mpp_mpi 
    302       ijpjjpk = jpj*jpk 
    303       CALL wrk_alloc( jpj*jpk, zwork ) 
    304 #endif 
    305  
    306       p_fval => p_fval2d 
    307  
    308       p_fval(:,:) = 0._wp 
    309       DO jk = 1, jpkm1 
    310          DO jj = 2, jpjm1 
    311             DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ? 
    312                p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * e1t(ji,jj) * fse3t(ji,jj,jk) * pmsk(ji,jj) 
    313             END DO 
    314          END DO 
    315       END DO 
    316 #if defined key_mpp_mpi 
    317       ijpjjpk = jpj*jpk 
    318       ish(1) = jpj*jpk   ;   ish2(1) = jpj   ;   ish2(2) = jpk 
    319       zwork(1:ijpjjpk)= RESHAPE( p_fval, ish ) 
    320       CALL mpp_sum( zwork, ijpjjpk, ncomm_znl ) 
    321       p_fval(:,:)= RESHAPE( zwork, ish2 ) 
    322 #endif 
    323       ! 
    324 #if defined key_mpp_mpi 
    325       CALL wrk_dealloc( jpj*jpk, zwork ) 
    326 #endif 
    327       !     
    328    END FUNCTION ptr_tjk 
    329  
    330  
    331    SUBROUTINE dia_ptr( kt ) 
    332       !!---------------------------------------------------------------------- 
    333       !!                  ***  ROUTINE dia_ptr  *** 
    334       !!---------------------------------------------------------------------- 
    335       USE oce,     vt  =>   ua   ! use ua as workspace 
    336       USE oce,     vs  =>   va   ! use va as workspace 
    337       IMPLICIT none 
    338       !! 
    339       INTEGER, INTENT(in) ::   kt   ! ocean time step index 
    340       ! 
    341       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    342       REAL(wp) ::   zv               ! local scalar 
    343       !!---------------------------------------------------------------------- 
    344       ! 
    345       IF( nn_timing == 1 )   CALL timing_start('dia_ptr') 
    346       ! 
    347       IF( kt == nit000 .OR. MOD( kt, nn_fptr ) == 0 )   THEN 
    348          ! 
    349          IF( MOD( kt, nn_fptr ) == 0 ) THEN  
    350             ! 
    351             IF( ln_diaznl ) THEN               ! i-mean temperature and salinity 
    352                DO jn = 1, nptr 
    353                   tn_jk(:,:,jn) = ptr_tjk( tsn(:,:,:,jp_tem), btmsk(:,:,jn) ) * r1_sjk(:,:,jn) 
    354                   sn_jk(:,:,jn) = ptr_tjk( tsn(:,:,:,jp_sal), btmsk(:,:,jn) ) * r1_sjk(:,:,jn) 
    355                END DO 
    356             ENDIF 
    357             ! 
    358             !                          ! horizontal integral and vertical dz  
    359             !                                ! eulerian velocity 
    360             v_msf(:,:,1) = ptr_vjk( vn(:,:,:) )  
    361             DO jn = 2, nptr 
    362                v_msf(:,:,jn) = ptr_vjk( vn(:,:,:), btmsk(:,:,jn)*btm30(:,:) )  
    363             END DO 
    364 #if defined key_diaeiv 
    365             DO jn = 1, nptr                  ! bolus velocity 
    366                v_msf_eiv(:,:,jn) = ptr_vjk( v_eiv(:,:,:), btmsk(:,:,jn) )   ! here no btm30 for MSFeiv 
    367             END DO 
    368             !                                ! add bolus stream-function to the eulerian one 
    369             v_msf(:,:,:) = v_msf(:,:,:) + v_msf_eiv(:,:,:) 
    370 #endif 
    371             ! 
    372             !                          ! Transports 
    373             !                                ! local heat & salt transports at T-points  ( tsn*mj[vn+v_eiv] ) 
    374             vt(:,:,jpk) = 0._wp   ;   vs(:,:,jpk) = 0._wp 
    375             DO jk= 1, jpkm1 
    376                DO jj = 2, jpj 
     106               DO ji = 1, jpi 
     107                  z3d(ji,:,:) = z3d(1,:,:) 
     108               ENDDO 
     109               cl1 = TRIM('zomsf'//clsubb(jn) ) 
     110               CALL iom_put( cl1, z3d * rc_sv ) 
     111            END DO 
     112         ENDIF 
     113         ! 
     114      ELSE 
     115         ! 
     116         IF( iom_use("zotemglo") ) THEN    ! i-mean i-k-surface  
     117            DO jk = 1, jpkm1 
     118               DO jj = 1, jpj 
    377119                  DO ji = 1, jpi 
    378 #if defined key_diaeiv  
    379                      zv = ( vn(ji,jj,jk) + vn(ji,jj-1,jk) + v_eiv(ji,jj,jk) + v_eiv(ji,jj-1,jk) ) * 0.5_wp 
    380 #else 
    381                      zv = ( vn(ji,jj,jk) + vn(ji,jj-1,jk) ) * 0.5_wp 
    382 #endif  
    383                      vt(ji,jj,jk) = zv * tsn(ji,jj,jk,jp_tem) 
    384                      vs(ji,jj,jk) = zv * tsn(ji,jj,jk,jp_sal) 
    385                   END DO 
    386                END DO 
    387             END DO 
    388 !!gm useless as overlap areas are not used in ptr_vjk 
    389             CALL lbc_lnk( vs, 'V', -1. )   ;   CALL lbc_lnk( vt, 'V', -1. ) 
    390 !!gm 
    391             !                                ! heat & salt advective transports (approximation) 
    392             htr(:,1) = SUM( ptr_vjk( vt(:,:,:) ) , 2 ) * rc_pwatt   ! SUM over jk + conversion 
    393             str(:,1) = SUM( ptr_vjk( vs(:,:,:) ) , 2 ) * rc_ggram 
    394             DO jn = 2, nptr  
    395                htr(:,jn) = SUM( ptr_vjk( vt(:,:,:), btmsk(:,:,jn)*btm30(:,:) ) , 2 ) * rc_pwatt   ! mask Southern Ocean 
    396                str(:,jn) = SUM( ptr_vjk( vs(:,:,:), btmsk(:,:,jn)*btm30(:,:) ) , 2 ) * rc_ggram   ! mask Southern Ocean 
    397             END DO 
    398  
    399             IF( ln_ptrcomp ) THEN            ! overturning transport 
    400                htr_ove(:) = SUM( v_msf(:,:,1) * tn_jk(:,:,1), 2 ) * rc_pwatt   ! SUM over jk + conversion 
    401                str_ove(:) = SUM( v_msf(:,:,1) * sn_jk(:,:,1), 2 ) * rc_ggram 
    402             END IF 
    403             !                                ! Advective and diffusive transport 
    404             htr_adv(:) = htr_adv(:) * rc_pwatt        ! these are computed in tra_adv... and tra_ldf... routines  
    405             htr_ldf(:) = htr_ldf(:) * rc_pwatt        ! here just the conversion in PW and Gg 
    406             str_adv(:) = str_adv(:) * rc_ggram 
    407             str_ldf(:) = str_ldf(:) * rc_ggram 
    408  
    409 #if defined key_diaeiv 
    410             DO jn = 1, nptr                  ! Bolus component 
    411                htr_eiv(:,jn) = SUM( v_msf_eiv(:,:,jn) * tn_jk(:,:,jn), 2 ) * rc_pwatt   ! SUM over jk 
    412                str_eiv(:,jn) = SUM( v_msf_eiv(:,:,jn) * sn_jk(:,:,jn), 2 ) * rc_ggram   ! SUM over jk 
    413             END DO 
    414 #endif 
    415             !                                ! "Meridional" Stream-Function 
     120                     zsfc = e1t(ji,jj) * fse3t(ji,jj,jk) 
     121                     zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc 
     122                     zts(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * zsfc 
     123                     zts(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * zsfc 
     124                  ENDDO 
     125               ENDDO 
     126            ENDDO 
    416127            DO jn = 1, nptr 
    417                DO jk = 2, jpk  
    418                   v_msf    (:,jk,jn) = v_msf    (:,jk-1,jn) + v_msf    (:,jk,jn)       ! Eulerian j-Stream-Function 
    419 #if defined key_diaeiv 
    420                   v_msf_eiv(:,jk,jn) = v_msf_eiv(:,jk-1,jn) + v_msf_eiv(:,jk,jn)       ! Bolus    j-Stream-Function 
    421  
    422 #endif 
    423                END DO 
    424             END DO 
    425             v_msf    (:,:,:) = v_msf    (:,:,:) * rc_sv       ! converte in Sverdrups 
    426 #if defined key_diaeiv 
    427             v_msf_eiv(:,:,:) = v_msf_eiv(:,:,:) * rc_sv 
    428 #endif 
    429          ENDIF 
    430          ! 
    431          CALL dia_ptr_wri( kt )                        ! outputs 
     128               zmask(1,:,:) = ptr_sjk( zmask(:,:,:), btmsk(:,:,jn) ) 
     129               cl1 = TRIM('zosrf'//clsubb(jn) ) 
     130               CALL iom_put( cl1, zmask ) 
     131               ! 
     132               z3d(1,:,:) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) & 
     133                  &            / MAX( zmask(1,:,:), 10.e-15 ) 
     134               DO ji = 1, jpi 
     135                  z3d(ji,:,:) = z3d(1,:,:) 
     136               ENDDO 
     137               cl1 = TRIM('zotem'//clsubb(jn) ) 
     138               CALL iom_put( cl1, z3d ) 
     139               ! 
     140               z3d(1,:,:) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) & 
     141                  &            / MAX( zmask(1,:,:), 10.e-15 ) 
     142               DO ji = 1, jpi 
     143                  z3d(ji,:,:) = z3d(1,:,:) 
     144               ENDDO 
     145               cl1 = TRIM('zosal'//clsubb(jn) ) 
     146               CALL iom_put( cl1, z3d ) 
     147            END DO 
     148         ENDIF 
     149         ! 
     150         !                                ! Advective and diffusive heat and salt transport 
     151         IF( iom_use("sophtadv") .OR. iom_use("sopstadv") ) THEN    
     152            z2d(1,:) = htr_adv(:) * rc_pwatt        !  (conversion in PW) 
     153            DO ji = 1, jpi 
     154               z2d(ji,:) = z2d(1,:) 
     155            ENDDO 
     156            cl1 = 'sophtadv'                  
     157            CALL iom_put( TRIM(cl1), z2d ) 
     158            z2d(1,:) = str_adv(:) * rc_ggram        ! (conversion in Gg) 
     159            DO ji = 1, jpi 
     160               z2d(ji,:) = z2d(1,:) 
     161            ENDDO 
     162            cl1 = 'sopstadv' 
     163            CALL iom_put( TRIM(cl1), z2d ) 
     164         ENDIF 
     165         ! 
     166         IF( iom_use("sophtldf") .OR. iom_use("sopstldf") ) THEN    
     167            z2d(1,:) = htr_ldf(:) * rc_pwatt        !  (conversion in PW)  
     168            DO ji = 1, jpi 
     169               z2d(ji,:) = z2d(1,:) 
     170            ENDDO 
     171            cl1 = 'sophtldf' 
     172            CALL iom_put( TRIM(cl1), z2d ) 
     173            z2d(1,:) = str_ldf(:) * rc_ggram        !  (conversion in Gg) 
     174            DO ji = 1, jpi 
     175               z2d(ji,:) = z2d(1,:) 
     176            ENDDO 
     177            cl1 = 'sopstldf' 
     178            CALL iom_put( TRIM(cl1), z2d ) 
     179         ENDIF 
    432180         ! 
    433181      ENDIF 
    434       ! 
    435 #if defined key_mpp_mpi 
    436       IF( kt == nitend .AND. l_znl_root )   CALL histclo( numptr )      ! Close the file 
    437 #else 
    438       IF( kt == nitend )                    CALL histclo( numptr )      ! Close the file 
    439 #endif 
    440182      ! 
    441183      IF( nn_timing == 1 )   CALL timing_stop('dia_ptr') 
     
    450192      !! ** Purpose :   Initialization, namelist read 
    451193      !!---------------------------------------------------------------------- 
    452       INTEGER ::   jn           ! dummy loop indices  
    453       INTEGER ::   inum, ierr   ! local integers 
    454       INTEGER ::   ios          ! Local integer output status for namelist read 
    455 #if defined key_mpp_mpi 
    456       INTEGER, DIMENSION(1) :: iglo, iloc, iabsf, iabsl, ihals, ihale, idid 
    457 #endif 
    458       !! 
    459       NAMELIST/namptr/ ln_diaptr, ln_diaznl, ln_subbas, ln_ptrcomp, nn_fptr, nn_fwri 
     194      INTEGER ::  jn           ! local integers 
     195      INTEGER ::  inum, ierr   ! local integers 
     196      INTEGER ::  ios          ! Local integer output status for namelist read 
     197      !! 
     198      NAMELIST/namptr/ ln_diaptr, ln_subbas 
    460199      !!---------------------------------------------------------------------- 
    461200 
     
    475214         WRITE(numout,*) '   Namelist namptr : set ptr parameters' 
    476215         WRITE(numout,*) '      Poleward heat & salt transport (T) or not (F)      ln_diaptr  = ', ln_diaptr 
    477          WRITE(numout,*) '      Overturning heat & salt transport                  ln_ptrcomp = ', ln_ptrcomp 
    478          WRITE(numout,*) '      T & S zonal mean and meridional stream function    ln_diaznl  = ', ln_diaznl  
    479216         WRITE(numout,*) '      Global (F) or glo/Atl/Pac/Ind/Indo-Pac basins      ln_subbas  = ', ln_subbas 
    480          WRITE(numout,*) '      Frequency of computation                           nn_fptr    = ', nn_fptr 
    481          WRITE(numout,*) '      Frequency of outputs                               nn_fwri    = ', nn_fwri 
    482217      ENDIF 
    483        
    484       IF( ln_diaptr) THEN   
    485       
    486          IF( nn_timing == 1 )   CALL timing_start('dia_ptr_init') 
    487        
    488          IF( ln_subbas ) THEN   ;   nptr = 5       ! Global, Atlantic, Pacific, Indian, Indo-Pacific 
    489          ELSE                   ;   nptr = 1       ! Global only 
     218 
     219      IF( ln_diaptr ) THEN   
     220         ! 
     221         IF( ln_subbas ) THEN  
     222            nptr = 5            ! Global, Atlantic, Pacific, Indian, Indo-Pacific 
     223            ALLOCATE( clsubb(nptr) ) 
     224            clsubb(1) = 'glo' ;  clsubb(2) = 'atl'  ;  clsubb(3) = 'pac'  ;  clsubb(4) = 'ind'  ;  clsubb(5) = 'ipc' 
     225         ELSE                
     226            nptr = 1       ! Global only 
     227            ALLOCATE( clsubb(nptr) ) 
     228            clsubb(1) = 'glo'  
    490229         ENDIF 
    491230 
     
    493232         IF( dia_ptr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ptr_init : unable to allocate arrays' ) 
    494233 
    495          rc_pwatt = rc_pwatt * rau0 * rcp          ! conversion from K.s-1 to PetaWatt 
     234         rc_pwatt = rc_pwatt * rau0_rcp          ! conversion from K.s-1 to PetaWatt 
    496235 
    497236         IF( lk_mpp )   CALL mpp_ini_znl( numout )     ! Define MPI communicator for zonal sum 
    498237 
    499238         IF( ln_subbas ) THEN                ! load sub-basin mask 
    500             CALL iom_open( 'subbasins', inum ) 
     239            CALL iom_open( 'subbasins', inum,  ldstop = .FALSE. ) 
    501240            CALL iom_get( inum, jpdom_data, 'atlmsk', btmsk(:,:,2) )   ! Atlantic basin 
    502241            CALL iom_get( inum, jpdom_data, 'pacmsk', btmsk(:,:,3) )   ! Pacific  basin 
     
    508247            END WHERE 
    509248         ENDIF 
     249    
    510250         btmsk(:,:,1) = tmask_i(:,:)                                   ! global ocean 
    511251       
     
    513253            btmsk(:,:,jn) = btmsk(:,:,jn) * tmask_i(:,:)               ! interior domain only 
    514254         END DO 
    515        
    516          IF( lk_vvl )   CALL ctl_stop( 'diaptr: error in vvl case as constant i-mean surface is used' ) 
    517  
    518          !                                   ! i-sum of e1v*e3v surface and its inverse 
    519          DO jn = 1, nptr 
    520             sjk(:,:,jn) = ptr_tjk( tmask(:,:,:), btmsk(:,:,jn) ) 
    521             r1_sjk(:,:,jn) = 0._wp 
    522             WHERE( sjk(:,:,jn) /= 0._wp )   r1_sjk(:,:,jn) = 1._wp / sjk(:,:,jn) 
    523          END DO 
    524  
    525       ! Initialise arrays to zero because diatpr is called before they are first calculated 
    526       ! Note that this means diagnostics will not be exactly correct when model run is restarted. 
    527       htr_adv(:) = 0._wp ; str_adv(:) =  0._wp ;  htr_ldf(:) = 0._wp ; str_ldf(:) =  0._wp 
    528  
    529 #if defined key_mpp_mpi  
    530          iglo (1) = jpjglo                   ! MPP case using MPI  ('key_mpp_mpi') 
    531          iloc (1) = nlcj 
    532          iabsf(1) = njmppt(narea) 
    533          iabsl(:) = iabsf(:) + iloc(:) - 1 
    534          ihals(1) = nldj - 1 
    535          ihale(1) = nlcj - nlej 
    536          idid (1) = 2 
    537          CALL flio_dom_set( jpnj, nproc/jpni, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom_ptr ) 
    538 #else 
    539          nidom_ptr = FLIO_DOM_NONE 
    540 #endif 
    541       IF( nn_timing == 1 )   CALL timing_stop('dia_ptr_init') 
    542       ! 
     255 
     256         ! Initialise arrays to zero because diatpr is called before they are first calculated 
     257         ! Note that this means diagnostics will not be exactly correct when model run is restarted. 
     258         htr_adv(:) = 0._wp  ;  str_adv(:) =  0._wp   
     259         htr_ldf(:) = 0._wp  ;  str_ldf(:) =  0._wp  
     260         ! 
    543261      ENDIF  
    544262      !  
     
    546264 
    547265 
    548    SUBROUTINE dia_ptr_wri( kt ) 
    549       !!--------------------------------------------------------------------- 
    550       !!                ***  ROUTINE dia_ptr_wri  *** 
    551       !! 
    552       !! ** Purpose :   output of poleward fluxes 
    553       !! 
    554       !! ** Method  :   NetCDF file 
    555       !!---------------------------------------------------------------------- 
    556       !! 
    557       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    558       !! 
    559       INTEGER, SAVE ::   nhoridz, ndepidzt, ndepidzw 
    560       INTEGER, SAVE ::   ndim  , ndim_atl     , ndim_pac     , ndim_ind     , ndim_ipc 
    561       INTEGER, SAVE ::           ndim_atl_30  , ndim_pac_30  , ndim_ind_30  , ndim_ipc_30 
    562       INTEGER, SAVE ::   ndim_h, ndim_h_atl_30, ndim_h_pac_30, ndim_h_ind_30, ndim_h_ipc_30 
    563       !! 
    564       CHARACTER (len=40) ::   clhstnam, clop, clop_once, cl_comment   ! temporary names 
    565       INTEGER            ::   iline, it, itmod, ji, jj, jk            ! 
    566 #if defined key_iomput 
    567       INTEGER            ::   inum                                    ! temporary logical unit 
     266   FUNCTION dia_ptr_alloc() 
     267      !!---------------------------------------------------------------------- 
     268      !!                    ***  ROUTINE dia_ptr_alloc  *** 
     269      !!---------------------------------------------------------------------- 
     270      INTEGER               ::   dia_ptr_alloc   ! return value 
     271      INTEGER, DIMENSION(3) ::   ierr 
     272      !!---------------------------------------------------------------------- 
     273      ierr(:) = 0 
     274      ! 
     275      ALLOCATE( btmsk(jpi,jpj,nptr) ,           & 
     276         &      htr_adv(jpj) , str_adv(jpj) ,   & 
     277         &      htr_ldf(jpj) , str_ldf(jpj) , STAT=ierr(1)  ) 
     278         ! 
     279      ALLOCATE( p_fval1d(jpj), p_fval2d(jpj,jpk), Stat=ierr(2)) 
     280      ! 
     281      ALLOCATE( btm30(jpi,jpj), STAT=ierr(3)  ) 
     282 
     283         ! 
     284      dia_ptr_alloc = MAXVAL( ierr ) 
     285      IF(lk_mpp)   CALL mpp_sum( dia_ptr_alloc ) 
     286      ! 
     287   END FUNCTION dia_ptr_alloc 
     288 
     289 
     290   FUNCTION ptr_sj_3d( pva, pmsk )   RESULT ( p_fval ) 
     291      !!---------------------------------------------------------------------- 
     292      !!                    ***  ROUTINE ptr_sj_3d  *** 
     293      !! 
     294      !! ** Purpose :   i-k sum computation of a j-flux array 
     295      !! 
     296      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i). 
     297      !!              pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v) 
     298      !! 
     299      !! ** Action  : - p_fval: i-k-mean poleward flux of pva 
     300      !!---------------------------------------------------------------------- 
     301      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)       ::   pva   ! mask flux array at V-point 
     302      REAL(wp), INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask 
     303      ! 
     304      INTEGER                  ::   ji, jj, jk   ! dummy loop arguments 
     305      INTEGER                  ::   ijpj         ! ??? 
     306      REAL(wp), POINTER, DIMENSION(:) :: p_fval  ! function value 
     307      !!-------------------------------------------------------------------- 
     308      ! 
     309      p_fval => p_fval1d 
     310 
     311      ijpj = jpj 
     312      p_fval(:) = 0._wp 
     313      IF( PRESENT( pmsk ) ) THEN  
     314         DO jk = 1, jpkm1 
     315            DO jj = 2, jpjm1 
     316               DO ji = fs_2, fs_jpim1   ! Vector opt. 
     317                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) * pmsk(ji,jj) 
     318               END DO 
     319            END DO 
     320         END DO 
     321      ELSE 
     322         DO jk = 1, jpkm1 
     323            DO jj = 2, jpjm1 
     324               DO ji = fs_2, fs_jpim1   ! Vector opt. 
     325                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj)  
     326               END DO 
     327            END DO 
     328         END DO 
     329      ENDIF 
     330#if defined key_mpp_mpi 
     331      IF(lk_mpp)   CALL mpp_sum( p_fval, ijpj, ncomm_znl) 
    568332#endif 
    569       REAL(wp)           ::   zsto, zout, zdt, zjulian                ! temporary scalars 
    570       !! 
    571       REAL(wp), POINTER, DIMENSION(:)   ::   zphi, zfoo    ! 1D workspace 
    572       REAL(wp), POINTER, DIMENSION(:,:) ::   z_1           ! 2D workspace 
    573       !!--------------------------------------------------------------------  
    574       ! 
    575       CALL wrk_alloc( jpj       , zphi , zfoo ) 
    576       CALL wrk_alloc( jpj , jpk , z_1  ) 
    577  
    578       ! define time axis 
    579       it    = kt / nn_fptr 
    580       itmod = kt - nit000 + 1 
    581        
    582       ! Initialization 
    583       ! -------------- 
    584       IF( kt == nit000 ) THEN 
    585          niter = ( nit000 - 1 ) / nn_fptr 
    586          zdt = rdt 
    587          IF( nacc == 1 )   zdt = rdtmin 
    588          ! 
    589          IF(lwp) THEN 
    590             WRITE(numout,*) 
    591             WRITE(numout,*) 'dia_ptr_wri : poleward transport and msf writing: initialization , niter = ', niter 
    592             WRITE(numout,*) '~~~~~~~~~~~~' 
    593          ENDIF 
    594  
    595          ! Reference latitude (used in plots) 
    596          ! ------------------ 
    597          !                                           ! ======================= 
    598          IF( cp_cfg == "orca" ) THEN                 !   ORCA configurations 
    599             !                                        ! ======================= 
    600             IF( jp_cfg == 05  )   iline = 192   ! i-line that passes near the North Pole 
    601             IF( jp_cfg == 025 )   iline = 384   ! i-line that passes near the North Pole 
    602             IF( jp_cfg == 1   )   iline =  96   ! i-line that passes near the North Pole 
    603             IF( jp_cfg == 2   )   iline =  48   ! i-line that passes near the North Pole 
    604             IF( jp_cfg == 4   )   iline =  24   ! i-line that passes near the North Pole 
    605             zphi(1:jpj) = 0._wp 
    606             DO ji = mi0(iline), mi1(iline)  
    607                zphi(1:jpj) = gphiv(ji,:)         ! if iline is in the local domain 
    608                ! Correct highest latitude for some configurations - will work if domain is parallelized in J ? 
    609                IF( jp_cfg == 05 ) THEN 
    610                   DO jj = mj0(jpjdta), mj1(jpjdta)  
    611                      zphi( jj ) = zphi(mj0(jpjdta-1)) + ( zphi(mj0(jpjdta-1))-zphi(mj0(jpjdta-2)) ) * 0.5_wp 
    612                      zphi( jj ) = MIN( zphi(jj), 90._wp ) 
    613                   END DO 
    614                END IF 
    615                IF( jp_cfg == 1 .OR. jp_cfg == 2 .OR. jp_cfg == 4 ) THEN 
    616                   DO jj = mj0(jpjdta-1), mj1(jpjdta-1)  
    617                      zphi( jj ) = 88.5_wp 
    618                   END DO 
    619                   DO jj = mj0(jpjdta  ), mj1(jpjdta  )  
    620                      zphi( jj ) = 89.5_wp 
    621                   END DO 
    622                END IF 
    623             END DO 
    624             ! provide the correct zphi to all local domains 
     333      ! 
     334   END FUNCTION ptr_sj_3d 
     335 
     336 
     337   FUNCTION ptr_sj_2d( pva, pmsk )   RESULT ( p_fval ) 
     338      !!---------------------------------------------------------------------- 
     339      !!                    ***  ROUTINE ptr_sj_2d  *** 
     340      !! 
     341      !! ** Purpose :   "zonal" and vertical sum computation of a i-flux array 
     342      !! 
     343      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i). 
     344      !!      pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v) 
     345      !! 
     346      !! ** Action  : - p_fval: i-k-mean poleward flux of pva 
     347      !!---------------------------------------------------------------------- 
     348      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)           ::   pva   ! mask flux array at V-point 
     349      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask 
     350      ! 
     351      INTEGER                  ::   ji,jj       ! dummy loop arguments 
     352      INTEGER                  ::   ijpj        ! ??? 
     353      REAL(wp), POINTER, DIMENSION(:) :: p_fval ! function value 
     354      !!-------------------------------------------------------------------- 
     355      !  
     356      p_fval => p_fval1d 
     357 
     358      ijpj = jpj 
     359      p_fval(:) = 0._wp 
     360      IF( PRESENT( pmsk ) ) THEN  
     361         DO jj = 2, jpjm1 
     362            DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ? 
     363               p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj) * pmsk(ji,jj) 
     364            END DO 
     365         END DO 
     366      ELSE 
     367         DO jj = 2, jpjm1 
     368            DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ? 
     369               p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj) 
     370            END DO 
     371         END DO 
     372      ENDIF 
    625373#if defined key_mpp_mpi 
    626             CALL mpp_sum( zphi, jpj, ncomm_znl )         
     374      CALL mpp_sum( p_fval, ijpj, ncomm_znl ) 
    627375#endif 
    628             !                                        ! ======================= 
    629          ELSE                                        !   OTHER configurations  
    630             !                                        ! ======================= 
    631             zphi(1:jpj) = gphiv(1,:)             ! assume lat/lon coordinate, select the first i-line 
    632             ! 
    633          ENDIF 
    634          ! 
    635          ! Work only on westmost processor (will not work if mppini2 is used) 
     376      !  
     377   END FUNCTION ptr_sj_2d 
     378 
     379 
     380   FUNCTION ptr_sjk( pta, pmsk )   RESULT ( p_fval ) 
     381      !!---------------------------------------------------------------------- 
     382      !!                    ***  ROUTINE ptr_sjk  *** 
     383      !! 
     384      !! ** Purpose :   i-sum computation of an array 
     385      !! 
     386      !! ** Method  : - i-sum of pva using the interior 2D vmask (vmask_i). 
     387      !! 
     388      !! ** Action  : - p_fval: i-mean poleward flux of pva 
     389      !!---------------------------------------------------------------------- 
     390      !! 
     391      IMPLICIT none 
     392      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk)           ::   pta    ! mask flux array at V-point 
     393      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)    , OPTIONAL ::   pmsk   ! Optional 2D basin mask 
     394      !! 
     395      INTEGER                           :: ji, jj, jk ! dummy loop arguments 
     396      REAL(wp), POINTER, DIMENSION(:,:) :: p_fval     ! return function value 
    636397#if defined key_mpp_mpi 
    637          IF( l_znl_root ) THEN  
     398      INTEGER, DIMENSION(1) ::   ish 
     399      INTEGER, DIMENSION(2) ::   ish2 
     400      INTEGER               ::   ijpjjpk 
     401      REAL(wp), DIMENSION(jpj*jpk) ::   zwork    ! mask flux array at V-point 
    638402#endif 
    639             ! 
    640             ! OPEN netcdf file  
    641             ! ---------------- 
    642             ! Define frequency of output and means 
    643             zsto = nn_fptr * zdt 
    644             IF( ln_mskland )   THEN    ! put 1.e+20 on land (very expensive!!) 
    645                clop      = "ave(only(x))" 
    646                clop_once = "once(only(x))" 
    647             ELSE                       ! no use of the mask value (require less cpu time) 
    648                clop      = "ave(x)"        
    649                clop_once = "once" 
    650             ENDIF 
    651  
    652             zout = nn_fwri * zdt 
    653             zfoo(1:jpj) = 0._wp 
    654  
    655             CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )  ! Compute julian date from starting date of the run 
    656             zjulian = zjulian - adatrj                         ! set calendar origin to the beginning of the experiment 
    657  
    658 #if defined key_iomput 
    659             ! Requested by IPSL people, use by their postpro... 
    660             IF(lwp) THEN 
    661                CALL dia_nam( clhstnam, nn_fwri,' ' ) 
    662                CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea ) 
    663                WRITE(inum,*) clhstnam 
    664                CLOSE(inum) 
    665             ENDIF 
     403      !!-------------------------------------------------------------------- 
     404      ! 
     405      p_fval => p_fval2d 
     406 
     407      p_fval(:,:) = 0._wp 
     408      ! 
     409      IF( PRESENT( pmsk ) ) THEN  
     410         DO jk = 1, jpkm1 
     411            DO jj = 2, jpjm1 
     412!!gm here, use of tmask_i  ==> no need of loop over nldi, nlei.... 
     413               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ? 
     414                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj) 
     415               END DO 
     416            END DO 
     417         END DO 
     418      ELSE  
     419         DO jk = 1, jpkm1 
     420            DO jj = 2, jpjm1 
     421               DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ? 
     422                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * tmask_i(ji,jj) 
     423               END DO 
     424            END DO 
     425         END DO 
     426      END IF 
     427      ! 
     428#if defined key_mpp_mpi 
     429      ijpjjpk = jpj*jpk 
     430      ish(1) = ijpjjpk  ;   ish2(1) = jpj   ;   ish2(2) = jpk 
     431      zwork(1:ijpjjpk) = RESHAPE( p_fval, ish ) 
     432      CALL mpp_sum( zwork, ijpjjpk, ncomm_znl ) 
     433      p_fval(:,:) = RESHAPE( zwork, ish2 ) 
    666434#endif 
    667  
    668             CALL dia_nam( clhstnam, nn_fwri, 'diaptr' ) 
    669             IF(lwp)WRITE( numout,*)" Name of diaptr NETCDF file : ", clhstnam 
    670  
    671             ! Horizontal grid : zphi() 
    672             CALL histbeg(clhstnam, 1, zfoo, jpj, zphi,   & 
    673                1, 1, 1, jpj, niter, zjulian, zdt*nn_fptr, nhoridz, numptr, domain_id=nidom_ptr) 
    674             ! Vertical grids : gdept_1d, gdepw_1d 
    675             CALL histvert( numptr, "deptht", "Vertical T levels",   & 
    676                &                   "m", jpk, gdept_1d, ndepidzt, "down" ) 
    677             CALL histvert( numptr, "depthw", "Vertical W levels",   & 
    678                &                   "m", jpk, gdepw_1d, ndepidzw, "down" ) 
    679             ! 
    680             CALL wheneq ( jpj*jpk, MIN(sjk(:,:,1), 1._wp), 1, 1., ndex  , ndim  )      ! Lat-Depth 
    681             CALL wheneq ( jpj    , MIN(sjk(:,1,1), 1._wp), 1, 1., ndex_h, ndim_h )     ! Lat 
    682  
    683             IF( ln_subbas ) THEN 
    684                z_1(:,1) = 1._wp 
    685                WHERE ( gphit(jpi/2,:) < -30._wp )   z_1(:,1) = 0._wp 
    686                DO jk = 2, jpk 
    687                   z_1(:,jk) = z_1(:,1) 
    688                END DO 
    689                !                       ! Atlantic (jn=2) 
    690                CALL wheneq ( jpj*jpk, MIN(sjk(:,:,2)         , 1._wp), 1, 1., ndex_atl     , ndim_atl      ) ! Lat-Depth 
    691                CALL wheneq ( jpj*jpk, MIN(sjk(:,:,2)*z_1(:,:), 1._wp), 1, 1., ndex_atl_30  , ndim_atl_30   ) ! Lat-Depth 
    692                CALL wheneq ( jpj    , MIN(sjk(:,1,2)*z_1(:,1), 1._wp), 1, 1., ndex_h_atl_30, ndim_h_atl_30 ) ! Lat 
    693                !                       ! Pacific (jn=3) 
    694                CALL wheneq ( jpj*jpk, MIN(sjk(:,:,3)         , 1._wp), 1, 1., ndex_pac     , ndim_pac      ) ! Lat-Depth 
    695                CALL wheneq ( jpj*jpk, MIN(sjk(:,:,3)*z_1(:,:), 1._wp), 1, 1., ndex_pac_30  , ndim_pac_30   ) ! Lat-Depth 
    696                CALL wheneq ( jpj    , MIN(sjk(:,1,3)*z_1(:,1), 1._wp), 1, 1., ndex_h_pac_30, ndim_h_pac_30 ) ! Lat 
    697                !                       ! Indian (jn=4) 
    698                CALL wheneq ( jpj*jpk, MIN(sjk(:,:,4)         , 1._wp), 1, 1., ndex_ind     , ndim_ind      ) ! Lat-Depth 
    699                CALL wheneq ( jpj*jpk, MIN(sjk(:,:,4)*z_1(:,:), 1._wp), 1, 1., ndex_ind_30  , ndim_ind_30   ) ! Lat-Depth 
    700                CALL wheneq ( jpj    , MIN(sjk(:,1,4)*z_1(:,1), 1._wp), 1, 1., ndex_h_ind_30, ndim_h_ind_30 ) ! Lat 
    701                !                       ! Indo-Pacific (jn=5) 
    702                CALL wheneq ( jpj*jpk, MIN(sjk(:,:,5)         , 1._wp), 1, 1., ndex_ipc     , ndim_ipc      ) ! Lat-Depth 
    703                CALL wheneq ( jpj*jpk, MIN(sjk(:,:,5)*z_1(:,:), 1._wp), 1, 1., ndex_ipc_30  , ndim_ipc_30   ) ! Lat-Depth 
    704                CALL wheneq ( jpj    , MIN(sjk(:,1,5)*z_1(:,1), 1._wp), 1, 1., ndex_h_ipc_30, ndim_h_ipc_30 ) ! Lat 
    705             ENDIF 
    706             !  
    707 #if defined key_diaeiv 
    708             cl_comment = ' (Bolus part included)' 
    709 #else 
    710             cl_comment = '                      ' 
    711 #endif 
    712             IF( ln_diaznl ) THEN             !  Zonal mean T and S 
    713                CALL histdef( numptr, "zotemglo", "Zonal Mean Temperature","C" ,   & 
    714                   1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    715                CALL histdef( numptr, "zosalglo", "Zonal Mean Salinity","PSU"  ,   & 
    716                   1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    717  
    718                CALL histdef( numptr, "zosrfglo", "Zonal Mean Surface","m^2"   ,   & 
    719                   1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout ) 
    720                ! 
    721                IF (ln_subbas) THEN  
    722                   CALL histdef( numptr, "zotematl", "Zonal Mean Temperature: Atlantic","C" ,   & 
    723                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    724                   CALL histdef( numptr, "zosalatl", "Zonal Mean Salinity: Atlantic","PSU"  ,   & 
    725                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    726                   CALL histdef( numptr, "zosrfatl", "Zonal Mean Surface: Atlantic","m^2"   ,   & 
    727                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout ) 
    728  
    729                   CALL histdef( numptr, "zotempac", "Zonal Mean Temperature: Pacific","C"  ,   & 
    730                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    731                   CALL histdef( numptr, "zosalpac", "Zonal Mean Salinity: Pacific","PSU"   ,   & 
    732                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    733                   CALL histdef( numptr, "zosrfpac", "Zonal Mean Surface: Pacific","m^2"    ,   & 
    734                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout ) 
    735  
    736                   CALL histdef( numptr, "zotemind", "Zonal Mean Temperature: Indian","C"   ,   & 
    737                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    738                   CALL histdef( numptr, "zosalind", "Zonal Mean Salinity: Indian","PSU"    ,   & 
    739                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    740                   CALL histdef( numptr, "zosrfind", "Zonal Mean Surface: Indian","m^2"     ,   & 
    741                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout ) 
    742  
    743                   CALL histdef( numptr, "zotemipc", "Zonal Mean Temperature: Pacific+Indian","C" ,   & 
    744                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    745                   CALL histdef( numptr, "zosalipc", "Zonal Mean Salinity: Pacific+Indian","PSU"  ,   & 
    746                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop, zsto, zout ) 
    747                   CALL histdef( numptr, "zosrfipc", "Zonal Mean Surface: Pacific+Indian","m^2"   ,   & 
    748                      1, jpj, nhoridz, jpk, 1, jpk, ndepidzt, 32, clop_once, zsto, zout ) 
    749                ENDIF 
    750             ENDIF 
    751             ! 
    752             !  Meridional Stream-Function (Eulerian and Bolus) 
    753             CALL histdef( numptr, "zomsfglo", "Meridional Stream-Function: Global"//TRIM(cl_comment),"Sv" ,   & 
    754                1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout ) 
    755             IF( ln_subbas .AND. ln_diaznl ) THEN 
    756                CALL histdef( numptr, "zomsfatl", "Meridional Stream-Function: Atlantic"//TRIM(cl_comment),"Sv" ,   & 
    757                   1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout ) 
    758                CALL histdef( numptr, "zomsfpac", "Meridional Stream-Function: Pacific"//TRIM(cl_comment),"Sv"  ,   & 
    759                   1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout ) 
    760                CALL histdef( numptr, "zomsfind", "Meridional Stream-Function: Indian"//TRIM(cl_comment),"Sv"   ,   & 
    761                   1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout ) 
    762                CALL histdef( numptr, "zomsfipc", "Meridional Stream-Function: Indo-Pacific"//TRIM(cl_comment),"Sv" ,& 
    763                   1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout ) 
    764             ENDIF 
    765             ! 
    766             !  Heat transport  
    767             CALL histdef( numptr, "sophtadv", "Advective Heat Transport"      ,   & 
    768                "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    769             CALL histdef( numptr, "sophtldf", "Diffusive Heat Transport"      ,   & 
    770                "PW",1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    771             IF ( ln_ptrcomp ) THEN  
    772                CALL histdef( numptr, "sophtove", "Overturning Heat Transport"    ,   & 
    773                   "PW",1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    774             END IF 
    775             IF( ln_subbas ) THEN 
    776                CALL histdef( numptr, "sohtatl", "Heat Transport Atlantic"//TRIM(cl_comment),  & 
    777                   "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    778                CALL histdef( numptr, "sohtpac", "Heat Transport Pacific"//TRIM(cl_comment) ,  & 
    779                   "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    780                CALL histdef( numptr, "sohtind", "Heat Transport Indian"//TRIM(cl_comment)  ,  & 
    781                   "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    782                CALL histdef( numptr, "sohtipc", "Heat Transport Pacific+Indian"//TRIM(cl_comment), & 
    783                   "PW", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    784             ENDIF 
    785             ! 
    786             !  Salt transport  
    787             CALL histdef( numptr, "sopstadv", "Advective Salt Transport"      ,   & 
    788                "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    789             CALL histdef( numptr, "sopstldf", "Diffusive Salt Transport"      ,   & 
    790                "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    791             IF ( ln_ptrcomp ) THEN  
    792                CALL histdef( numptr, "sopstove", "Overturning Salt Transport"    ,   & 
    793                   "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    794             END IF 
    795 #if defined key_diaeiv 
    796             ! Eddy induced velocity 
    797             CALL histdef( numptr, "zomsfeiv", "Bolus Meridional Stream-Function: global",   & 
    798                "Sv"      , 1, jpj, nhoridz, jpk, 1, jpk, ndepidzw, 32, clop, zsto, zout ) 
    799             CALL histdef( numptr, "sophteiv", "Bolus Advective Heat Transport",   & 
    800                "PW"      , 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    801             CALL histdef( numptr, "sopsteiv", "Bolus Advective Salt Transport",   & 
    802                "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    803 #endif 
    804             IF( ln_subbas ) THEN 
    805                CALL histdef( numptr, "sostatl", "Salt Transport Atlantic"//TRIM(cl_comment)      ,  & 
    806                   "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    807                CALL histdef( numptr, "sostpac", "Salt Transport Pacific"//TRIM(cl_comment)      ,   & 
    808                   "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    809                CALL histdef( numptr, "sostind", "Salt Transport Indian"//TRIM(cl_comment)      ,    & 
    810                   "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    811                CALL histdef( numptr, "sostipc", "Salt Transport Pacific+Indian"//TRIM(cl_comment),  & 
    812                   "Giga g/s", 1, jpj, nhoridz, 1, 1, 1, -99, 32, clop, zsto, zout ) 
    813             ENDIF 
    814             ! 
    815             CALL histend( numptr ) 
    816             ! 
    817          END IF 
    818 #if defined key_mpp_mpi 
    819       END IF 
    820 #endif 
    821  
    822 #if defined key_mpp_mpi 
    823       IF( MOD( itmod, nn_fptr ) == 0 .AND. l_znl_root ) THEN 
    824 #else 
    825       IF( MOD( itmod, nn_fptr ) == 0  ) THEN 
    826 #endif 
    827          niter = niter + 1 
    828  
    829          IF( ln_diaznl ) THEN  
    830             CALL histwrite( numptr, "zosrfglo", niter, sjk  (:,:,1) , ndim, ndex ) 
    831             CALL histwrite( numptr, "zotemglo", niter, tn_jk(:,:,1)  , ndim, ndex ) 
    832             CALL histwrite( numptr, "zosalglo", niter, sn_jk(:,:,1)  , ndim, ndex ) 
    833  
    834             IF (ln_subbas) THEN  
    835                CALL histwrite( numptr, "zosrfatl", niter, sjk(:,:,2), ndim_atl, ndex_atl ) 
    836                CALL histwrite( numptr, "zosrfpac", niter, sjk(:,:,3), ndim_pac, ndex_pac ) 
    837                CALL histwrite( numptr, "zosrfind", niter, sjk(:,:,4), ndim_ind, ndex_ind ) 
    838                CALL histwrite( numptr, "zosrfipc", niter, sjk(:,:,5), ndim_ipc, ndex_ipc ) 
    839  
    840                CALL histwrite( numptr, "zotematl", niter, tn_jk(:,:,2)  , ndim_atl, ndex_atl ) 
    841                CALL histwrite( numptr, "zosalatl", niter, sn_jk(:,:,2)  , ndim_atl, ndex_atl ) 
    842                CALL histwrite( numptr, "zotempac", niter, tn_jk(:,:,3)  , ndim_pac, ndex_pac ) 
    843                CALL histwrite( numptr, "zosalpac", niter, sn_jk(:,:,3)  , ndim_pac, ndex_pac ) 
    844                CALL histwrite( numptr, "zotemind", niter, tn_jk(:,:,4)  , ndim_ind, ndex_ind ) 
    845                CALL histwrite( numptr, "zosalind", niter, sn_jk(:,:,4)  , ndim_ind, ndex_ind ) 
    846                CALL histwrite( numptr, "zotemipc", niter, tn_jk(:,:,5)  , ndim_ipc, ndex_ipc ) 
    847                CALL histwrite( numptr, "zosalipc", niter, sn_jk(:,:,5)  , ndim_ipc, ndex_ipc ) 
    848             END IF 
    849          ENDIF 
    850  
    851          ! overturning outputs: 
    852          CALL histwrite( numptr, "zomsfglo", niter, v_msf(:,:,1), ndim, ndex ) 
    853          IF( ln_subbas .AND. ln_diaznl ) THEN 
    854             CALL histwrite( numptr, "zomsfatl", niter, v_msf(:,:,2) , ndim_atl_30, ndex_atl_30 ) 
    855             CALL histwrite( numptr, "zomsfpac", niter, v_msf(:,:,3) , ndim_pac_30, ndex_pac_30 ) 
    856             CALL histwrite( numptr, "zomsfind", niter, v_msf(:,:,4) , ndim_ind_30, ndex_ind_30 ) 
    857             CALL histwrite( numptr, "zomsfipc", niter, v_msf(:,:,5) , ndim_ipc_30, ndex_ipc_30 ) 
    858          ENDIF 
    859 #if defined key_diaeiv 
    860          CALL histwrite( numptr, "zomsfeiv", niter, v_msf_eiv(:,:,1), ndim  , ndex   ) 
    861 #endif 
    862  
    863          ! heat transport outputs: 
    864          IF( ln_subbas ) THEN 
    865             CALL histwrite( numptr, "sohtatl", niter, htr(:,2)  , ndim_h_atl_30, ndex_h_atl_30 ) 
    866             CALL histwrite( numptr, "sohtpac", niter, htr(:,3)  , ndim_h_pac_30, ndex_h_pac_30 ) 
    867             CALL histwrite( numptr, "sohtind", niter, htr(:,4)  , ndim_h_ind_30, ndex_h_ind_30 ) 
    868             CALL histwrite( numptr, "sohtipc", niter, htr(:,5)  , ndim_h_ipc_30, ndex_h_ipc_30 ) 
    869             CALL histwrite( numptr, "sostatl", niter, str(:,2)  , ndim_h_atl_30, ndex_h_atl_30 ) 
    870             CALL histwrite( numptr, "sostpac", niter, str(:,3)  , ndim_h_pac_30, ndex_h_pac_30 ) 
    871             CALL histwrite( numptr, "sostind", niter, str(:,4)  , ndim_h_ind_30, ndex_h_ind_30 ) 
    872             CALL histwrite( numptr, "sostipc", niter, str(:,5)  , ndim_h_ipc_30, ndex_h_ipc_30 ) 
    873          ENDIF 
    874  
    875          CALL histwrite( numptr, "sophtadv", niter, htr_adv     , ndim_h, ndex_h ) 
    876          CALL histwrite( numptr, "sophtldf", niter, htr_ldf     , ndim_h, ndex_h ) 
    877          CALL histwrite( numptr, "sopstadv", niter, str_adv     , ndim_h, ndex_h ) 
    878          CALL histwrite( numptr, "sopstldf", niter, str_ldf     , ndim_h, ndex_h ) 
    879          IF( ln_ptrcomp ) THEN  
    880             CALL histwrite( numptr, "sopstove", niter, str_ove(:) , ndim_h, ndex_h ) 
    881             CALL histwrite( numptr, "sophtove", niter, htr_ove(:) , ndim_h, ndex_h ) 
    882          ENDIF 
    883 #if defined key_diaeiv 
    884          CALL histwrite( numptr, "sophteiv", niter, htr_eiv(:,1)  , ndim_h, ndex_h ) 
    885          CALL histwrite( numptr, "sopsteiv", niter, str_eiv(:,1)  , ndim_h, ndex_h ) 
    886 #endif 
    887          ! 
    888       ENDIF 
    889       ! 
    890       CALL wrk_dealloc( jpj      , zphi , zfoo ) 
    891       CALL wrk_dealloc( jpj , jpk, z_1 ) 
    892       ! 
    893   END SUBROUTINE dia_ptr_wri 
     435      ! 
     436   END FUNCTION ptr_sjk 
     437 
    894438 
    895439   !!====================================================================== 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r5038 r5620  
    4444   USE in_out_manager  ! I/O manager 
    4545   USE diadimg         ! dimg direct access file format output 
    46 <<<<<<< .working 
    47    USE diaar5, ONLY :   lk_diaar5 
    48 ======= 
    49 >>>>>>> .merge-right.r5035 
    5046   USE iom 
    5147   USE ioipsl 
     48   USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities      
     49 
    5250#if defined key_lim2 
    5351   USE limwri_2  
     
    8280   !!---------------------------------------------------------------------- 
    8381   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    84    !! $Id $ 
     82   !! $Id$ 
    8583   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    8684   !!---------------------------------------------------------------------- 
     
    129127      !! 
    130128      INTEGER                      ::   ji, jj, jk              ! dummy loop indices 
     129      INTEGER                      ::   jkbot                   ! 
    131130      REAL(wp)                     ::   zztmp, zztmpx, zztmpy   !  
    132131      !! 
    133 <<<<<<< .working 
    134132      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace 
    135       REAL(wp), POINTER, DIMENSION(:,:)   :: z2ds     ! 2D workspace 
    136 ======= 
    137       REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace 
    138 >>>>>>> .merge-right.r5035 
    139133      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace 
    140134      !!---------------------------------------------------------------------- 
     
    142136      IF( nn_timing == 1 )   CALL timing_start('dia_wri') 
    143137      !  
    144       CALL wrk_alloc( jpi , jpj      , z2d , z2ds ) 
     138      CALL wrk_alloc( jpi , jpj      , z2d ) 
    145139      CALL wrk_alloc( jpi , jpj, jpk , z3d ) 
    146140      ! 
     
    151145      ENDIF 
    152146 
    153       IF( lk_vvl ) THEN 
    154          z3d(:,:,:) = tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) 
    155          CALL iom_put( "toce" , z3d                        )   ! heat content 
     147      IF( .NOT.lk_vvl ) THEN 
     148         CALL iom_put( "e3t" , fse3t_n(:,:,:) ) 
     149         CALL iom_put( "e3u" , fse3u_n(:,:,:) ) 
     150         CALL iom_put( "e3v" , fse3v_n(:,:,:) ) 
     151         CALL iom_put( "e3w" , fse3w_n(:,:,:) ) 
     152      ENDIF 
     153 
     154      CALL iom_put( "ssh" , sshn )                 ! sea surface height 
     155      if( iom_use('ssh2') )   CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height 
     156       
     157      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature 
     158      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature 
     159      IF ( iom_use("sbt") ) THEN 
    156160         DO jj = 1, jpj 
    157161            DO ji = 1, jpi 
    158                z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * fse3t_n(ji,jj,mikt(ji,jj)) 
    159             END DO 
    160          END DO   
    161          CALL iom_put( "sst"  , z2d(:,:)                 )   ! sea surface heat content       
     162               jkbot = mbkt(ji,jj) 
     163               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem) 
     164            END DO 
     165         END DO 
     166         CALL iom_put( "sbt", z2d )                ! bottom temperature 
     167      ENDIF 
     168       
     169      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity 
     170      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity 
     171      IF ( iom_use("sbs") ) THEN 
    162172         DO jj = 1, jpj 
    163173            DO ji = 1, jpi 
    164                z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem)**2 * fse3t_n(ji,jj,mikt(ji,jj)) 
    165             END DO 
    166          END DO   
    167          CALL iom_put( "sst2" , z2d(:,:)      )   ! sea surface content of squared temperature 
    168          z3d(:,:,:) = tsn(:,:,:,jp_sal) * fse3t_n(:,:,:)             
    169          CALL iom_put( "soce" , z3d                        )   ! salinity content 
     174               jkbot = mbkt(ji,jj) 
     175               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_sal) 
     176            END DO 
     177         END DO 
     178         CALL iom_put( "sbs", z2d )                ! bottom salinity 
     179      ENDIF 
     180 
     181      IF ( iom_use("taubot") ) THEN                ! bottom stress 
     182         z2d(:,:) = 0._wp 
     183         DO jj = 2, jpjm1 
     184            DO ji = fs_2, fs_jpim1   ! vector opt. 
     185               zztmpx = (  bfrua(ji  ,jj) * un(ji  ,jj,mbku(ji  ,jj))  & 
     186                      &  + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj))  )       
     187               zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  & 
     188                      &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  )  
     189               z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1)  
     190               ! 
     191            ENDDO 
     192         ENDDO 
     193         CALL lbc_lnk( z2d, 'T', 1. ) 
     194         CALL iom_put( "taubot", z2d )            
     195      ENDIF 
     196          
     197      CALL iom_put( "uoce", un(:,:,:)         )    ! 3D i-current 
     198      CALL iom_put(  "ssu", un(:,:,1)         )    ! surface i-current 
     199      IF ( iom_use("sbu") ) THEN 
    170200         DO jj = 1, jpj 
    171201            DO ji = 1, jpi 
    172                z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * fse3t_n(ji,jj,mikt(ji,jj)) 
    173             END DO 
    174          END DO   
    175          CALL iom_put( "sss"  , z2d(:,:)                 )   ! sea surface salinity content 
     202               jkbot = mbku(ji,jj) 
     203               z2d(ji,jj) = un(ji,jj,jkbot) 
     204            END DO 
     205         END DO 
     206         CALL iom_put( "sbu", z2d )                ! bottom i-current 
     207      ENDIF 
     208#if defined key_dynspg_ts 
     209      CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current 
     210#else 
     211      CALL iom_put(  "ubar", un_b(:,:)        )    ! barotropic i-current 
     212#endif 
     213       
     214      CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current 
     215      CALL iom_put(  "ssv", vn(:,:,1)         )    ! surface j-current 
     216      IF ( iom_use("sbv") ) THEN 
    176217         DO jj = 1, jpj 
    177218            DO ji = 1, jpi 
    178                z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal)**2 * fse3t_n(ji,jj,mikt(ji,jj)) 
    179             END DO 
    180          END DO   
    181          CALL iom_put( "sss2" , z2d(:,:)                 )   ! sea surface content of squared salinity 
    182       ELSE 
    183          CALL iom_put( "toce" , tsn(:,:,:,jp_tem)        )   ! temperature 
    184          IF ( iom_use("sst") ) THEN 
    185             DO jj = 1, jpj 
    186                DO ji = 1, jpi 
    187                   z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) 
    188                END DO 
    189             END DO 
    190             CALL iom_put( "sst"  , z2d(:,:)            ) ! sea surface temperature 
    191          ENDIF 
    192          IF ( iom_use("sst2") )   CALL iom_put( "sst2" , z2d(:,:) * z2d(:,:) ) ! square of sea surface temperature 
    193          CALL iom_put( "soce" , tsn(:,:,:,jp_sal)          )   ! salinity 
    194          IF ( iom_use("sss") ) THEN 
    195             DO jj = 1, jpj 
    196                DO ji = 1, jpi 
    197                   z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) 
    198                END DO 
    199             END DO 
    200             CALL iom_put( "sss"  , z2d(:,:)            ) ! sea surface salinity 
    201          ENDIF 
    202          CALL iom_put( "sss2" , z2d(:,:) * z2d(:,:) ) ! square of sea surface salinity 
    203       END IF 
    204       IF( lk_vvl .AND. (.NOT. ln_dynadv_vec) ) THEN 
    205          CALL iom_put( "uoce" , umask(:,:,:) * un(:,:,:) * fse3u_n(:,:,:) )    ! i-transport 
    206          CALL iom_put( "voce" , vmask(:,:,:) * vn(:,:,:) * fse3v_n(:,:,:) )    ! j-transport 
    207       ELSE 
    208          CALL iom_put( "uoce" , umask(:,:,:) * un(:,:,:)                  )    ! i-current 
    209          CALL iom_put( "voce" , vmask(:,:,:) * vn(:,:,:)                  )    ! j-current 
    210          IF ( iom_use("ssu") ) THEN 
    211             DO jj = 1, jpj 
    212                DO ji = 1, jpi 
    213                   z2d(ji,jj) = un(ji,jj,miku(ji,jj)) 
    214                END DO 
    215             END DO 
    216             CALL iom_put( "ssu"   , z2d                                    )    ! i-current 
    217          ENDIF 
    218          IF ( iom_use("ssv") ) THEN 
    219             DO jj = 1, jpj 
    220                DO ji = 1, jpi 
    221                   z2d(ji,jj) = vn(ji,jj,mikv(ji,jj)) 
    222                END DO 
    223             END DO 
    224             CALL iom_put( "ssv"   , z2d                                    )    ! j-current 
    225          ENDIF 
    226       ENDIF 
    227       CALL iom_put(    "avt"  , avt                        )    ! T vert. eddy diff. coef. 
    228       CALL iom_put(    "avm"  , avmu                       )    ! T vert. eddy visc. coef. 
    229       IF( lk_zdfddm ) THEN 
    230          CALL iom_put( "avs" , fsavs(:,:,:)                          )    ! S vert. eddy diff. coef. 
    231       ENDIF 
    232  
    233       IF ( iom_use("sstgrad2") .OR. iom_use("sstgrad2") ) THEN 
     219               jkbot = mbkv(ji,jj) 
     220               z2d(ji,jj) = vn(ji,jj,jkbot) 
     221            END DO 
     222         END DO 
     223         CALL iom_put( "sbv", z2d )                ! bottom j-current 
     224      ENDIF 
     225#if defined key_dynspg_ts 
     226      CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic j-current 
     227#else 
     228      CALL iom_put(  "vbar", vn_b(:,:)        )    ! barotropic j-current 
     229#endif 
     230 
     231      CALL iom_put( "woce", wn )                   ! vertical velocity 
     232      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value 
     233         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 
     234         z2d(:,:) = rau0 * e12t(:,:) 
     235         DO jk = 1, jpk 
     236            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) 
     237         END DO 
     238         CALL iom_put( "w_masstr" , z3d )   
     239         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 
     240      ENDIF 
     241 
     242      CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef. 
     243      CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef. 
     244      CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm) 
     245 
     246      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN 
    234247         DO jj = 2, jpjm1                                    ! sst gradient 
    235248            DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    243256         CALL lbc_lnk( z2d, 'T', 1. ) 
    244257         CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient 
    245          !CDIR NOVERRCHK< 
    246258         z2d(:,:) = SQRT( z2d(:,:) ) 
    247259         CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient 
     
    252264         z2d(:,:)  = 0._wp  
    253265         DO jk = 1, jpkm1 
    254             DO jj = 2, jpjm1 
    255                DO ji = fs_2, fs_jpim1   ! vector opt. 
     266            DO jj = 1, jpj 
     267               DO ji = 1, jpi 
    256268                  z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 
    257269               END DO 
    258270            END DO 
    259271         END DO 
    260          CALL lbc_lnk( z2d, 'T', 1. ) 
    261272         CALL iom_put( "heatc", (rau0 * rcp) * z2d )    ! vertically integrated heat content (J/m2) 
    262273      ENDIF 
    263274 
    264 <<<<<<< .working 
    265       ! clem: heat and salt content 
    266       z2d(:,:)  = 0._wp  
    267       z2ds(:,:) = 0._wp  
    268       DO jk = 1, jpkm1 
    269          DO jj = 2, jpjm1 
    270             DO ji = fs_2, fs_jpim1   ! vector opt. 
    271                z2d(ji,jj) = z2d(ji,jj) + rau0 * rcp * fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 
    272                z2ds(ji,jj) = z2ds(ji,jj) + rau0 * fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk) 
    273             END DO 
    274          END DO 
    275       END DO 
    276       CALL lbc_lnk( z2d, 'T', 1. ) 
    277       CALL lbc_lnk( z2ds, 'T', 1. ) 
    278       CALL iom_put( "heatc", z2d )    ! vertically integrated heat content (J/m2) 
    279       CALL iom_put( "saltc", z2ds )   ! vertically integrated salt content (PSU*kg/m2) 
    280        
    281  
    282       IF( lk_diaar5 ) THEN 
    283 ======= 
    284275      IF( iom_use("saltc") ) THEN 
    285276         z2d(:,:)  = 0._wp  
    286277         DO jk = 1, jpkm1 
    287             DO jj = 2, jpjm1 
    288                DO ji = fs_2, fs_jpim1   ! vector opt. 
     278            DO jj = 1, jpj 
     279               DO ji = 1, jpi 
    289280                  z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk) 
    290281               END DO 
    291282            END DO 
    292283         END DO 
    293          CALL lbc_lnk( z2d, 'T', 1. ) 
    294284         CALL iom_put( "saltc", rau0 * z2d )   ! vertically integrated salt content (PSU*kg/m2) 
    295285      ENDIF 
     
    319309          
    320310      IF( iom_use("u_masstr") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN 
    321 >>>>>>> .merge-right.r5035 
    322311         z3d(:,:,jpk) = 0.e0 
    323312         DO jk = 1, jpkm1 
     
    325314         END DO 
    326315         CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction 
    327 <<<<<<< .working 
    328  
    329          zztmp = 0.5 * rcp 
    330 ======= 
    331316      ENDIF 
    332317       
    333318      IF( iom_use("u_heattr") ) THEN 
    334 >>>>>>> .merge-right.r5035 
    335319         z2d(:,:) = 0.e0  
    336          z2ds(:,:) = 0.e0  
    337320         DO jk = 1, jpkm1 
    338321            DO jj = 2, jpjm1 
    339322               DO ji = fs_2, fs_jpim1   ! vector opt. 
    340323                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) ) 
    341                   z2ds(ji,jj) = z2ds(ji,jj) + z3d(ji,jj,jk) * 0.5_wp * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) ) 
    342324               END DO 
    343325            END DO 
    344326         END DO 
    345327         CALL lbc_lnk( z2d, 'U', -1. ) 
    346 <<<<<<< .working 
    347          CALL lbc_lnk( z2ds, 'U', -1. ) 
    348          CALL iom_put( "u_heattr", z2d )                  ! heat transport in i-direction 
    349 ======= 
    350328         CALL iom_put( "u_heattr", (0.5 * rcp) * z2d )    ! heat transport in i-direction 
    351329      ENDIF 
     
    353331      IF( iom_use("u_salttr") ) THEN 
    354332         z2d(:,:) = 0.e0  
    355 >>>>>>> .merge-right.r5035 
    356          CALL iom_put( "u_salttr", z2ds )                 ! salt transport in i-direction 
    357  
    358          z3d(:,:,jpk) = 0.e0 
    359333         DO jk = 1, jpkm1 
    360 <<<<<<< .working 
    361             z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk) 
    362 ======= 
    363334            DO jj = 2, jpjm1 
    364335               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    366337               END DO 
    367338            END DO 
    368 >>>>>>> .merge-right.r5035 
    369339         END DO 
    370340         CALL lbc_lnk( z2d, 'U', -1. ) 
     
    379349         END DO 
    380350         CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction 
    381 <<<<<<< .working 
    382  
    383 ======= 
    384351      ENDIF 
    385352       
    386353      IF( iom_use("v_heattr") ) THEN 
    387 >>>>>>> .merge-right.r5035 
    388354         z2d(:,:) = 0.e0  
    389          z2ds(:,:) = 0.e0  
    390355         DO jk = 1, jpkm1 
    391356            DO jj = 2, jpjm1 
    392357               DO ji = fs_2, fs_jpim1   ! vector opt. 
    393358                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) ) 
    394                   z2ds(ji,jj) = z2ds(ji,jj) + z3d(ji,jj,jk) * 0.5_wp * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) ) 
    395359               END DO 
    396360            END DO 
    397361         END DO 
    398362         CALL lbc_lnk( z2d, 'V', -1. ) 
    399 <<<<<<< .working 
    400          CALL lbc_lnk( z2ds, 'V', -1. ) 
    401          CALL iom_put( "v_heattr", z2d )                  !  heat transport in j-direction 
    402          CALL iom_put( "v_salttr", z2ds )                 !  salt transport in j-direction 
    403 ======= 
    404363         CALL iom_put( "v_heattr", (0.5 * rcp) * z2d )    !  heat transport in j-direction 
    405 >>>>>>> .merge-right.r5035 
    406364      ENDIF 
    407365 
     
    419377      ENDIF 
    420378      ! 
    421       CALL wrk_dealloc( jpi , jpj      , z2d , z2ds ) 
     379      CALL wrk_dealloc( jpi , jpj      , z2d ) 
    422380      CALL wrk_dealloc( jpi , jpj, jpk , z3d ) 
    423381      ! 
     
    480438      zdt = rdt 
    481439      IF( nacc == 1 ) zdt = rdtmin 
    482       IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!) 
    483       ELSE                      ;   clop = "x"         ! no use of the mask value (require less cpu time) 
    484       ENDIF 
     440      clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes) 
    485441#if defined key_diainstant 
    486442      zsto = nwrite * zdt 
     
    682638         ENDIF 
    683639 
    684          IF( .NOT. lk_cpl ) THEN 
     640         IF( .NOT. ln_cpl ) THEN 
    685641            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp 
    686642               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     
    691647         ENDIF 
    692648 
    693          IF( lk_cpl .AND. nn_ice <= 1 ) THEN 
     649         IF( ln_cpl .AND. nn_ice <= 1 ) THEN 
    694650            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp 
    695651               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     
    714670#endif 
    715671 
    716          IF( lk_cpl .AND. nn_ice == 2 ) THEN 
     672         IF( ln_cpl .AND. nn_ice == 2 ) THEN 
    717673            CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice 
    718674               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout ) 
     
    869825      ENDIF 
    870826 
    871       IF( .NOT. lk_cpl ) THEN 
     827      IF( .NOT. ln_cpl ) THEN 
    872828         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping 
    873829         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping 
     
    875831         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
    876832      ENDIF 
    877       IF( lk_cpl .AND. nn_ice <= 1 ) THEN 
     833      IF( ln_cpl .AND. nn_ice <= 1 ) THEN 
    878834         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping 
    879835         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping 
     
    891847#endif 
    892848 
    893       IF( lk_cpl .AND. nn_ice == 2 ) THEN 
     849      IF( ln_cpl .AND. nn_ice == 2 ) THEN 
    894850         CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature 
    895851         CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90

    r4162 r5620  
    7272      !!---------------------------------------------------------------------- 
    7373      INTEGER ::   jc            ! dummy loop indices 
     74      INTEGER :: isrow           ! local index 
    7475      !!---------------------------------------------------------------------- 
    7576       
     
    9192         CASE ( 1 )                                  ! ORCA_R1 configuration 
    9293            !                                        ! ======================= 
     94            ! This dirty section will be suppressed by simplification process: 
     95            ! all this will come back in input files 
     96            ! Currently these hard-wired indices relate to configuration with 
     97            ! extend grid (jpjglo=332) 
     98            isrow = 332 - jpjglo 
     99            ! 
    93100            ncsnr(1)   = 1    ; ncstt(1)   = 0           ! Caspian Sea 
    94             ncsi1(1)   = 332  ; ncsj1(1)   = 203 
    95             ncsi2(1)   = 344  ; ncsj2(1)   = 235 
     101            ncsi1(1)   = 332  ; ncsj1(1)   = 243 - isrow 
     102            ncsi2(1)   = 344  ; ncsj2(1)   = 275 - isrow 
    96103            ncsir(1,1) = 1    ; ncsjr(1,1) = 1 
    97104            !                                         
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90

    r5038 r5620  
    7373      !!---------------------------------------------------------------------- 
    7474      ! 
     75      ! max number of seconds between each restart 
     76      IF( REAL( nitend - nit000 + 1 ) * rdt > REAL( HUGE( nsec1jan000 ) ) ) THEN 
     77         CALL ctl_stop( 'The number of seconds between each restart exceeds the integer 4 max value: 2^31-1. ',   & 
     78            &           'You must do a restart at higher frequency (or remove this stop and recompile the code in I8)' ) 
     79      ENDIF 
    7580      ! all calendar staff is based on the fact that MOD( rday, rdttra(1) ) == 0 
    7681      IF( MOD( rday     , rdttra(1) ) /= 0. )   CALL ctl_stop( 'the time step must devide the number of second of in a day' ) 
     
    238243               nday_year = 1 
    239244               nsec_year = ndt05 
    240                IF( nsec1jan000 >= 2 * (2**30 - nsecd * nyear_len(1) / 2 ) ) THEN   ! test integer 4 max value 
    241                   CALL ctl_stop( 'The number of seconds between Jan. 1st 00h of nit000 year and Jan. 1st 00h ',   & 
    242                      &           'of the current year is exceeding the INTEGER 4 max VALUE: 2^31-1 -> 68.09 years in seconds', & 
    243                      & 'You must do a restart at higher frequency (or remove this STOP and recompile everything in I8)' ) 
    244                ENDIF 
    245245               nsec1jan000 = nsec1jan000 + nsecd * nyear_len(1) 
    246246               IF( nleapy == 1 )   CALL day_mth 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r5038 r5620  
    162162   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  gphit, gphiu   !: latitude  of t-, u-, v- and f-points (degre) 
    163163   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  gphiv, gphif   !: 
    164    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1t, e2t       !: horizontal scale factors at t-point (m) 
    165    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1u, e2u       !: horizontal scale factors at u-point (m) 
    166    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1v, e2v       !: horizontal scale factors at v-point (m) 
    167    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1f, e2f       !: horizontal scale factors at f-point (m) 
     164   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1t, e2t, r1_e1t, r1_e2t   !: horizontal scale factors and inverse at t-point (m) 
     165   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1u, e2u, r1_e1u, r1_e2u   !: horizontal scale factors and inverse at u-point (m) 
     166   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1v, e2v, r1_e1v, r1_e2v   !: horizontal scale factors and inverse at v-point (m) 
     167   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, TARGET, DIMENSION(:,:) ::  e1f, e2f, r1_e1f, r1_e2f   !: horizontal scale factors and inverse at f-point (m) 
    168168   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  e1e2t          !: surface at t-point (m2) 
    169169   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  ff             !: coriolis factor (2.*omega*sin(yphi) ) (s-1) 
     
    262262 
    263263   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts 
     264   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: wmask, wumask, wvmask        !: land/ocean mask at WT-, WU- and WV-pts 
    264265 
    265266   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tpol, fpol          !: north fold mask (jperio= 3 or 4) 
     
    332333   INTEGER FUNCTION dom_oce_alloc() 
    333334      !!---------------------------------------------------------------------- 
    334       INTEGER, DIMENSION(11) :: ierr 
     335      INTEGER, DIMENSION(12) :: ierr 
    335336      !!---------------------------------------------------------------------- 
    336337      ierr(:) = 0 
     
    345346         &      tpol(jpiglo)  , fpol(jpiglo)                               , STAT=ierr(2) ) 
    346347         ! 
    347       ALLOCATE( glamt(jpi,jpj) , gphit(jpi,jpj) , e1t(jpi,jpj) , e2t(jpi,jpj) ,                      &  
    348          &      glamu(jpi,jpj) , gphiu(jpi,jpj) , e1u(jpi,jpj) , e2u(jpi,jpj) ,                      &   
    349          &      glamv(jpi,jpj) , gphiv(jpi,jpj) , e1v(jpi,jpj) , e2v(jpi,jpj) , e1e2t(jpi,jpj) ,     &   
    350          &      glamf(jpi,jpj) , gphif(jpi,jpj) , e1f(jpi,jpj) , e2f(jpi,jpj) , ff   (jpi,jpj) , STAT=ierr(3) )      
     348      ALLOCATE( glamt(jpi,jpj) , gphit(jpi,jpj) , e1t(jpi,jpj) , e2t(jpi,jpj) , r1_e1t(jpi,jpj) , r1_e2t(jpi,jpj) ,   &  
     349         &      glamu(jpi,jpj) , gphiu(jpi,jpj) , e1u(jpi,jpj) , e2u(jpi,jpj) , r1_e1u(jpi,jpj) , r1_e2u(jpi,jpj) ,   &   
     350         &      glamv(jpi,jpj) , gphiv(jpi,jpj) , e1v(jpi,jpj) , e2v(jpi,jpj) , r1_e1v(jpi,jpj) , r1_e2v(jpi,jpj) ,   &   
     351         &      glamf(jpi,jpj) , gphif(jpi,jpj) , e1f(jpi,jpj) , e2f(jpi,jpj) , r1_e1f(jpi,jpj) , r1_e2f(jpi,jpj) ,   & 
     352         &      e1e2t(jpi,jpj) , ff   (jpi,jpj) , STAT=ierr(3) )      
    351353         ! 
    352354      ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 (jpi,jpj,jpk) ,                         & 
     
    400402         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk), STAT=ierr(11) ) 
    401403 
     404      ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) ) 
     405 
    402406#if defined key_noslip_accurate 
    403       ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(11) ) 
     407      ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(12) ) 
    404408#endif 
    405409      ! 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r5038 r5620  
    135135      !!---------------------------------------------------------------------- 
    136136      USE ioipsl 
    137       NAMELIST/namrun/ nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   & 
     137      NAMELIST/namrun/ cn_ocerst_indir, cn_ocerst_outdir, nn_stocklist, ln_rst_list,               & 
     138         &             nn_no   , cn_exp    , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl,   & 
    138139         &             nn_it000, nn_itend  , nn_date0    , nn_leapy     , nn_istate , nn_stock ,   & 
    139          &             nn_write, ln_dimgnnn, ln_mskland  , ln_clobber   , nn_chunksz, nn_euler 
     140         &             nn_write, ln_dimgnnn, ln_mskland  , ln_cfmeta    , ln_clobber, nn_chunksz, nn_euler 
    140141      NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin,   & 
    141142         &             nn_acc   , rn_atfp     , rn_rdt      , rn_rdtmin ,                  & 
     
    169170         WRITE(numout,*) '      experiment name for output      cn_exp     = ', cn_exp 
    170171         WRITE(numout,*) '      file prefix restart input       cn_ocerst_in= ', cn_ocerst_in 
     172         WRITE(numout,*) '      restart input directory         cn_ocerst_indir= ', cn_ocerst_indir 
    171173         WRITE(numout,*) '      file prefix restart output      cn_ocerst_out= ', cn_ocerst_out 
     174         WRITE(numout,*) '      restart output directory        cn_ocerst_outdir= ', cn_ocerst_outdir 
    172175         WRITE(numout,*) '      restart logical                 ln_rstart  = ', ln_rstart 
    173176         WRITE(numout,*) '      start with forward time step    nn_euler   = ', nn_euler 
     
    178181         WRITE(numout,*) '      leap year calendar (0/1)        nn_leapy   = ', nn_leapy 
    179182         WRITE(numout,*) '      initial state output            nn_istate  = ', nn_istate 
    180          WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock 
     183         IF( ln_rst_list ) THEN 
     184            WRITE(numout,*) '      list of restart dump times      nn_stocklist   =', nn_stocklist 
     185         ELSE 
     186            WRITE(numout,*) '      frequency of restart file       nn_stock   = ', nn_stock 
     187         ENDIF 
    181188         WRITE(numout,*) '      frequency of output file        nn_write   = ', nn_write 
    182189         WRITE(numout,*) '      multi file dimgout              ln_dimgnnn = ', ln_dimgnnn 
    183190         WRITE(numout,*) '      mask land points                ln_mskland = ', ln_mskland 
     191         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta  = ', ln_cfmeta 
    184192         WRITE(numout,*) '      overwrite an existing file      ln_clobber = ', ln_clobber 
    185193         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz = ', nn_chunksz 
     
    195203      ninist = nn_istate 
    196204      nstock = nn_stock 
     205      nstocklist = nn_stocklist 
    197206      nwrite = nn_write 
    198207      neuler = nn_euler 
    199       IF ( neuler == 1 .AND. .NOT.ln_rstart ) THEN 
     208      IF ( neuler == 1 .AND. .NOT. ln_rstart ) THEN 
    200209         WRITE(ctmp1,*) 'ln_rstart =.FALSE., nn_euler is forced to 0 ' 
    201210         CALL ctl_warn( ctmp1 ) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90

    r5038 r5620  
    105105      REAL(wp) ::   zlam1, zcos_alpha, zim1 , zjm1 , ze1, ze1deg 
    106106      REAL(wp) ::   zphi1, zsin_alpha, zim05, zjm05 
     107      INTEGER  ::   isrow                ! index for ORCA1 starting row 
     108 
    107109      !!---------------------------------------------------------------------- 
    108110      ! 
     
    159161         IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration 
    160162            !                                             ! ===================== 
    161  
    162             ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u = 20 km) 
    163             ij0 = 200   ;   ij1 = 200   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3 
     163            ! This dirty section will be suppressed by simplification process: all this will come back in input files 
     164            ! Currently these hard-wired indices relate to configuration with 
     165            ! extend grid (jpjglo=332) 
     166            ! which had a grid-size of 362x292. 
     167            !  
     168            isrow = 332 - jpjglo 
     169            ! 
     170            ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait (e2u = 20 km) 
     171            ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  20.e3 
    164172            IF(lwp) WRITE(numout,*) 
    165173            IF(lwp) WRITE(numout,*) '             orca_r1: Gibraltar : e2u reduced to 20 km' 
    166174 
    167             ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u = 10 km) 
    168             ij0 = 208   ;   ij1 = 208   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3 
     175            ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u = 10 km) 
     176            ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   e2u( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  10.e3 
    169177            IF(lwp) WRITE(numout,*) 
    170178            IF(lwp) WRITE(numout,*) '             orca_r1: Bhosporus : e2u reduced to 10 km' 
    171179 
    172             ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v = 13 km) 
    173             ij0 = 124   ;   ij1 = 125   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  13.e3 
     180            ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v = 13 km) 
     181            ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  13.e3 
    174182            IF(lwp) WRITE(numout,*) 
    175183            IF(lwp) WRITE(numout,*) '             orca_r1: Lombok : e1v reduced to 10 km' 
    176184 
    177             ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on] 
    178             ij0 = 124   ;   ij1 = 125   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  8.e3 
     185            ii0 =  48           ;   ii1 =  48        ! Sumba Strait (e1v = 8 km) [closed from bathy_11 on] 
     186            ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  8.e3 
    179187            IF(lwp) WRITE(numout,*) 
    180188            IF(lwp) WRITE(numout,*) '             orca_r1: Sumba : e1v reduced to 8 km' 
    181189 
    182             ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v = 13 km) 
    183             ij0 = 124   ;   ij1 = 125   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3 
     190            ii0 =  53           ;   ii1 =  53        ! Ombai Strait (e1v = 13 km) 
     191            ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 13.e3 
    184192            IF(lwp) WRITE(numout,*) 
    185193            IF(lwp) WRITE(numout,*) '             orca_r1: Ombai : e1v reduced to 13 km' 
    186194 
    187             ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v = 20 km) 
    188             ij0 = 124   ;   ij1 = 125   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3 
     195            ii0 =  56           ;   ii1 =  56        ! Timor Passage (e1v = 20 km) 
     196            ij0 = 164 - isrow   ;   ij1 = 145 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 20.e3 
    189197            IF(lwp) WRITE(numout,*) 
    190198            IF(lwp) WRITE(numout,*) '             orca_r1: Timor Passage : e1v reduced to 20 km' 
    191199 
    192             ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v = 30 km) 
    193             ij0 = 141   ;   ij1 = 142   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3 
     200            ii0 =  55           ;   ii1 =  55        ! West Halmahera Strait (e1v = 30 km) 
     201            ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 30.e3 
    194202            IF(lwp) WRITE(numout,*) 
    195203            IF(lwp) WRITE(numout,*) '             orca_r1: W Halmahera : e1v reduced to 30 km' 
    196204 
    197             ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v = 50 km) 
    198             ij0 = 141   ;   ij1 = 142   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3 
     205            ii0 =  58           ;   ii1 =  58        ! East Halmahera Strait (e1v = 50 km) 
     206            ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   e1v( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 50.e3 
    199207            IF(lwp) WRITE(numout,*) 
    200208            IF(lwp) WRITE(numout,*) '             orca_r1: E Halmahera : e1v reduced to 50 km' 
    201  
    202             ! 
    203  
    204             ! 
    205             ! 
    206209            ! 
    207210            ! 
     
    471474      re2u_e1u(:,:) = e2u(:,:) / e1u(:,:) 
    472475      re1v_e2v(:,:) = e1v(:,:) / e2v(:,:) 
     476      r1_e1t  (:,:) = 1._wp    / e1t(:,:) 
     477      r1_e1u  (:,:) = 1._wp    / e1u(:,:) 
     478      r1_e1v  (:,:) = 1._wp    / e1v(:,:) 
     479      r1_e1f  (:,:) = 1._wp    / e1f(:,:) 
     480      r1_e2t  (:,:) = 1._wp    / e2t(:,:) 
     481      r1_e2u  (:,:) = 1._wp    / e2u(:,:) 
     482      r1_e2v  (:,:) = 1._wp    / e2v(:,:) 
     483      r1_e2f  (:,:) = 1._wp    / e2f(:,:) 
    473484 
    474485      ! Control printing : Grid informations (if not restart) 
     
    616627      CALL iom_open( 'coordinates', inum ) 
    617628       
    618       CALL iom_get( inum, jpdom_data, 'glamt', glamt ) 
    619       CALL iom_get( inum, jpdom_data, 'glamu', glamu ) 
    620       CALL iom_get( inum, jpdom_data, 'glamv', glamv ) 
    621       CALL iom_get( inum, jpdom_data, 'glamf', glamf ) 
    622        
    623       CALL iom_get( inum, jpdom_data, 'gphit', gphit ) 
    624       CALL iom_get( inum, jpdom_data, 'gphiu', gphiu ) 
    625       CALL iom_get( inum, jpdom_data, 'gphiv', gphiv ) 
    626       CALL iom_get( inum, jpdom_data, 'gphif', gphif ) 
    627        
    628       CALL iom_get( inum, jpdom_data, 'e1t', e1t ) 
    629       CALL iom_get( inum, jpdom_data, 'e1u', e1u ) 
    630       CALL iom_get( inum, jpdom_data, 'e1v', e1v ) 
    631       CALL iom_get( inum, jpdom_data, 'e1f', e1f ) 
    632        
    633       CALL iom_get( inum, jpdom_data, 'e2t', e2t ) 
    634       CALL iom_get( inum, jpdom_data, 'e2u', e2u ) 
    635       CALL iom_get( inum, jpdom_data, 'e2v', e2v ) 
    636       CALL iom_get( inum, jpdom_data, 'e2f', e2f ) 
     629      CALL iom_get( inum, jpdom_data, 'glamt', glamt, lrowattr=ln_use_jattr ) 
     630      CALL iom_get( inum, jpdom_data, 'glamu', glamu, lrowattr=ln_use_jattr ) 
     631      CALL iom_get( inum, jpdom_data, 'glamv', glamv, lrowattr=ln_use_jattr ) 
     632      CALL iom_get( inum, jpdom_data, 'glamf', glamf, lrowattr=ln_use_jattr ) 
     633       
     634      CALL iom_get( inum, jpdom_data, 'gphit', gphit, lrowattr=ln_use_jattr ) 
     635      CALL iom_get( inum, jpdom_data, 'gphiu', gphiu, lrowattr=ln_use_jattr ) 
     636      CALL iom_get( inum, jpdom_data, 'gphiv', gphiv, lrowattr=ln_use_jattr ) 
     637      CALL iom_get( inum, jpdom_data, 'gphif', gphif, lrowattr=ln_use_jattr ) 
     638       
     639      CALL iom_get( inum, jpdom_data, 'e1t', e1t, lrowattr=ln_use_jattr ) 
     640      CALL iom_get( inum, jpdom_data, 'e1u', e1u, lrowattr=ln_use_jattr ) 
     641      CALL iom_get( inum, jpdom_data, 'e1v', e1v, lrowattr=ln_use_jattr ) 
     642      CALL iom_get( inum, jpdom_data, 'e1f', e1f, lrowattr=ln_use_jattr ) 
     643       
     644      CALL iom_get( inum, jpdom_data, 'e2t', e2t, lrowattr=ln_use_jattr ) 
     645      CALL iom_get( inum, jpdom_data, 'e2u', e2u, lrowattr=ln_use_jattr ) 
     646      CALL iom_get( inum, jpdom_data, 'e2v', e2v, lrowattr=ln_use_jattr ) 
     647      CALL iom_get( inum, jpdom_data, 'e2f', e2f, lrowattr=ln_use_jattr ) 
    637648       
    638649      CALL iom_close( inum ) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r5038 r5620  
    134134      INTEGER  ::   ijf, ijl, ij0, ij1       !   -       - 
    135135      INTEGER  ::   ios 
     136      INTEGER  ::   isrow                    ! index for ORCA1 starting row 
    136137      INTEGER , POINTER, DIMENSION(:,:) ::  imsk 
    137138      REAL(wp), POINTER, DIMENSION(:,:) ::  zwf 
     
    281282      CALL lbc_lnk( fmask_i, 'F', 1._wp ) 
    282283 
     284      ! 3. Ocean/land mask at wu-, wv- and w points  
     285      !---------------------------------------------- 
     286      wmask (:,:,1) = tmask(:,:,1) ! ???????? 
     287      wumask(:,:,1) = umask(:,:,1) ! ???????? 
     288      wvmask(:,:,1) = vmask(:,:,1) ! ???????? 
     289      DO jk=2,jpk 
     290         wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1) 
     291         wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)    
     292         wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1) 
     293      END DO 
    283294 
    284295      ! 4. ocean/land mask for the elliptic equation 
     
    391402      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN   ! ORCA R1 configuration 
    392403         !                                                 ! Increased lateral friction near of some straits 
     404         ! This dirty section will be suppressed by simplification process: 
     405         ! all this will come back in input files 
     406         ! Currently these hard-wired indices relate to configuration with 
     407         ! extend grid (jpjglo=332) 
     408         ! 
     409         isrow = 332 - jpjglo 
     410         ! 
    393411         IF(lwp) WRITE(numout,*) 
    394412         IF(lwp) WRITE(numout,*) '   orca_r1: increase friction near the following straits : ' 
    395413         IF(lwp) WRITE(numout,*) '      Gibraltar ' 
    396          ii0 = 283   ;   ii1 = 284        ! Gibraltar Strait  
    397          ij0 = 200   ;   ij1 = 200   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp   
     414         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait  
     415         ij0 = 241 - isrow   ;   ij1 = 241 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp   
    398416 
    399417         IF(lwp) WRITE(numout,*) '      Bhosporus ' 
    400          ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait  
    401          ij0 = 208   ;   ij1 = 208   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp   
     418         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait  
     419         ij0 = 248 - isrow   ;   ij1 = 248 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp   
    402420 
    403421         IF(lwp) WRITE(numout,*) '      Makassar (Top) ' 
    404          ii0 =  48   ;   ii1 =  48        ! Makassar Strait (Top)  
    405          ij0 = 149   ;   ij1 = 150   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp   
     422         ii0 =  48           ;   ii1 =  48        ! Makassar Strait (Top)  
     423         ij0 = 189 - isrow   ;   ij1 = 190 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp   
    406424 
    407425         IF(lwp) WRITE(numout,*) '      Lombok ' 
    408          ii0 =  44   ;   ii1 =  44        ! Lombok Strait  
    409          ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp   
     426         ii0 =  44           ;   ii1 =  44        ! Lombok Strait  
     427         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp   
    410428 
    411429         IF(lwp) WRITE(numout,*) '      Ombai ' 
    412          ii0 =  53   ;   ii1 =  53        ! Ombai Strait  
    413          ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp   
     430         ii0 =  53           ;   ii1 =  53        ! Ombai Strait  
     431         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp   
    414432 
    415433         IF(lwp) WRITE(numout,*) '      Timor Passage ' 
    416          ii0 =  56   ;   ii1 =  56        ! Timor Passage  
    417          ij0 = 124   ;   ij1 = 125   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 2._wp   
     434         ii0 =  56           ;   ii1 =  56        ! Timor Passage  
     435         ij0 = 164 - isrow   ;   ij1 = 165 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 2._wp   
    418436 
    419437         IF(lwp) WRITE(numout,*) '      West Halmahera ' 
    420          ii0 =  58   ;   ii1 =  58        ! West Halmahera Strait  
    421          ij0 = 141   ;   ij1 = 142   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp   
     438         ii0 =  58           ;   ii1 =  58        ! West Halmahera Strait  
     439         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp   
    422440 
    423441         IF(lwp) WRITE(numout,*) '      East Halmahera ' 
    424          ii0 =  55   ;   ii1 =  55        ! East Halmahera Strait  
    425          ij0 = 141   ;   ij1 = 142   ;   fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) = 3._wp   
     442         ii0 =  55           ;   ii1 =  55        ! East Halmahera Strait  
     443         ij0 = 181 - isrow   ;   ij1 = 182 - isrow   ;   fmask( mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1),1:jpk ) = 3._wp   
    426444         ! 
    427445      ENDIF 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r5038 r5620  
    88   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl: 
    99   !!                                          vvl option includes z_star and z_tilde coordinates 
     10   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
    1011   !!---------------------------------------------------------------------- 
    1112   !!   'key_vvl'                              variable volume 
     
    125126      INTEGER ::   ji,jj,jk 
    126127      INTEGER ::   ii0, ii1, ij0, ij1 
     128      REAL(wp)::   zcoef 
    127129      !!---------------------------------------------------------------------- 
    128130      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_init') 
     
    164166      ! t- and w- points depth 
    165167      ! ---------------------- 
     168      ! set the isf depth as it is in the initial step 
    166169      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 
    167170      fsdepw_n(:,:,1) = 0.0_wp 
     
    169172      fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1) 
    170173      fsdepw_b(:,:,1) = 0.0_wp 
    171       DO jj = 1,jpj 
    172          DO ji = 1,jpi 
    173             DO jk = 2,mikt(ji,jj)-1 
    174                fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 
    175                fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    176                fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 
    177                fsdept_b(ji,jj,jk) = gdept_0(ji,jj,jk) 
    178                fsdepw_b(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    179             END DO 
    180             IF (mikt(ji,jj) .GT. 1) THEN 
    181                jk = mikt(ji,jj) 
    182                fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk) 
    183                fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    184                fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
    185                fsdept_b(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_b(ji,jj,jk) 
    186                fsdepw_b(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    187             END IF 
    188             DO jk = mikt(ji,jj)+1, jpk 
    189                fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk) 
     174 
     175      DO jk = 2, jpk 
     176         DO jj = 1,jpj 
     177            DO ji = 1,jpi 
     178              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     179                                                     ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 
     180                                                     ! 0.5 where jk = mikt   
     181               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    190182               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 
    191                fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
    192                fsdept_b(ji,jj,jk) = fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk) 
     183               fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  & 
     184                   &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) +       fse3w_n(ji,jj,jk))  
     185               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 
    193186               fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1) 
     187               fsdept_b(ji,jj,jk) =      zcoef  * ( fsdepw_b(ji,jj,jk  ) + 0.5 * fse3w_b(ji,jj,jk))  & 
     188                   &                + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) +       fse3w_b(ji,jj,jk))  
    194189            END DO 
    195190         END DO 
     
    588583      INTEGER, INTENT( in )               :: kt       ! time step 
    589584      !! * Local declarations 
    590       REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def 
    591585      INTEGER                             :: ji,jj,jk       ! dummy loop indices 
     586      REAL(wp)                            :: zcoef 
    592587      !!---------------------------------------------------------------------- 
    593588 
    594589      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp') 
    595       ! 
    596       CALL wrk_alloc( jpi, jpj, jpk, z_e3t_def                ) 
    597590      ! 
    598591      IF( kt == nit000 )   THEN 
     
    638631      ! t- and w- points depth 
    639632      ! ---------------------- 
     633      ! set the isf depth as it is in the initial step 
    640634      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 
    641635      fsdepw_n(:,:,1) = 0.0_wp 
    642636      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 
    643       DO jj = 1,jpj 
    644          DO ji = 1,jpi 
    645             DO jk = 2,mikt(ji,jj)-1 
    646                fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 
    647                fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    648                fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 
    649             END DO 
    650             IF (mikt(ji,jj) .GT. 1) THEN 
    651                jk = mikt(ji,jj) 
    652                fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk) 
    653                fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    654                fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
    655             END IF 
    656             DO jk = mikt(ji,jj)+1, jpk 
    657                fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk) 
     637 
     638      DO jk = 2, jpk 
     639         DO jj = 1,jpj 
     640            DO ji = 1,jpi 
     641              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 
     642                                                                 ! 1 for jk = mikt 
     643               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
    658644               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 
    659                fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
     645               fsdept_n(ji,jj,jk) =      zcoef  * ( fsdepw_n(ji,jj,jk  ) + 0.5 * fse3w_n(ji,jj,jk))  & 
     646                   &                + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) +       fse3w_n(ji,jj,jk))  
     647               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 
    660648            END DO 
    661649         END DO 
    662650      END DO 
     651 
    663652      ! Local depth and Inverse of the local depth of the water column at u- and v- points 
    664653      ! ---------------------------------------------------------------------------------- 
     
    679668      ! Write outputs 
    680669      ! ============= 
    681       z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 
    682       CALL iom_put( "cellthc" , fse3t_n  (:,:,:) ) 
     670      CALL iom_put(     "e3t" , fse3t_n  (:,:,:) ) 
     671      CALL iom_put(     "e3u" , fse3u_n  (:,:,:) ) 
     672      CALL iom_put(     "e3v" , fse3v_n  (:,:,:) ) 
     673      CALL iom_put(     "e3w" , fse3w_n  (:,:,:) ) 
    683674      CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) ) 
    684       CALL iom_put( "e3tdef"  , z_e3t_def(:,:,:) ) 
     675      IF( iom_use("e3tdef") )   & 
     676         CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 
    685677 
    686678      ! write restart file 
    687679      ! ================== 
    688680      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) 
    689       ! 
    690       CALL wrk_dealloc( jpi, jpj, jpk, z_e3t_def ) 
    691681      ! 
    692682      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_swp') 
     
    10491039      INTEGER ::   ji, jj, jk                                          ! dummy loop indices 
    10501040      INTEGER ::   ij0, ij1, ii0, ii1                                  ! dummy loop indices 
     1041      INTEGER ::   isrow                                               ! index for ORCA1 starting row 
    10511042      !! acc 
    10521043      !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for 
     
    11321123      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration 
    11331124         !                                             ! ===================== 
    1134          ! 
    1135          ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u was modified) 
    1136          ij0 = 200   ;   ij1 = 200 
     1125         ! This dirty section will be suppressed by simplification process: 
     1126         ! all this will come back in input files 
     1127         ! Currently these hard-wired indices relate to configuration with 
     1128         ! extend grid (jpjglo=332) 
     1129         ! which had a grid-size of 362x292. 
     1130         isrow = 332 - jpjglo 
     1131         ! 
     1132         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait (e2u was modified) 
     1133         ij0 = 241 - isrow   ;   ij1 = 241 - isrow 
    11371134         DO jk = 1, jpkm1 
    11381135            DO jj = mj0(ij0), mj1(ij1) 
     
    11541151         END DO 
    11551152         ! 
    1156          ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u was modified) 
    1157          ij0 = 208   ;   ij1 = 208 
     1153         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u was modified) 
     1154         ij0 = 248 - isrow   ;   ij1 = 248 - isrow 
    11581155         DO jk = 1, jpkm1 
    11591156            DO jj = mj0(ij0), mj1(ij1) 
     
    11751172         END DO 
    11761173         ! 
    1177          ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v was modified) 
    1178          ij0 = 124   ;   ij1 = 125 
     1174         ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v was modified) 
     1175         ij0 = 164 - isrow   ;   ij1 = 165 - isrow 
    11791176         DO jk = 1, jpkm1 
    11801177            DO jj = mj0(ij0), mj1(ij1) 
     
    11911188         END DO 
    11921189         ! 
    1193          ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on] 
    1194          ij0 = 124   ;   ij1 = 125 
     1190         ii0 =  48           ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on] 
     1191         ij0 = 164 - isrow   ;   ij1 = 165 - isrow 
    11951192         DO jk = 1, jpkm1 
    11961193            DO jj = mj0(ij0), mj1(ij1) 
     
    12071204         END DO 
    12081205         ! 
    1209          ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v was modified) 
    1210          ij0 = 124   ;   ij1 = 125 
     1206         ii0 =  53          ;   ii1 =  53        ! Ombai Strait (e1v was modified) 
     1207         ij0 = 164 - isrow  ;   ij1 = 165  - isrow   
    12111208         DO jk = 1, jpkm1 
    12121209            DO jj = mj0(ij0), mj1(ij1) 
     
    12231220         END DO 
    12241221         ! 
    1225          ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v was modified) 
    1226          ij0 = 124   ;   ij1 = 125 
     1222         ii0 =  56            ;   ii1 =  56        ! Timor Passage (e1v was modified) 
     1223         ij0 = 164 - isrow    ;   ij1 = 165  - isrow   
    12271224         DO jk = 1, jpkm1 
    12281225            DO jj = mj0(ij0), mj1(ij1) 
     
    12391236         END DO 
    12401237         ! 
    1241          ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v was modified) 
    1242          ij0 = 141   ;   ij1 = 142 
     1238         ii0 =  55            ;   ii1 =  55        ! West Halmahera Strait (e1v was modified) 
     1239         ij0 = 181 - isrow    ;   ij1 = 182 - isrow   
    12431240         DO jk = 1, jpkm1 
    12441241            DO jj = mj0(ij0), mj1(ij1) 
     
    12551252         END DO 
    12561253         ! 
    1257          ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v was modified) 
    1258          ij0 = 141   ;   ij1 = 142 
     1254         ii0 =  58            ;   ii1 =  58        ! East Halmahera Strait (e1v was modified) 
     1255         ij0 = 181 - isrow    ;   ij1 = 182 - isrow   
    12591256         DO jk = 1, jpkm1 
    12601257            DO jj = mj0(ij0), mj1(ij1) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r5038 r5620  
    215215         CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d )  !    ! stretched system 
    216216         CALL iom_rstput( 0, 0, inum4, 'gdepw_1d' , gdepw_1d ) 
     217         CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r4 )      
     218         CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r4 ) 
    217219      ENDIF 
    218220       
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5038 r5620  
    1717   !!            3.4  ! 2012-08  (J. Siddorn) added Siddorn and Furner stretching function 
    1818   !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  modify C1D case   
     19   !!            3.6  ! 2014-11  (P. Mathiot and C. Harris) add ice shelf capabilitye   
    1920   !!---------------------------------------------------------------------- 
    2021 
     
    3536   USE oce               ! ocean variables 
    3637   USE dom_oce           ! ocean domain 
    37    USE sbc_oce           ! surface variable (isf) 
    3838   USE closea            ! closed seas 
    3939   USE c1d               ! 1D vertical configuration 
     
    298298      ENDIF 
    299299 
     300      IF ( ln_isfcav ) THEN 
    300301! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 
    301302! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 
    302       DO jk = 1, jpkm1 
    303          e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)  
    304       END DO 
    305       e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO 
    306  
    307       DO jk = 2, jpk 
    308          e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)  
    309       END DO 
    310       e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))  
     303         DO jk = 1, jpkm1 
     304            e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)  
     305         END DO 
     306         e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO 
     307 
     308         DO jk = 2, jpk 
     309            e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)  
     310         END DO 
     311         e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))  
     312      END IF 
    311313 
    312314!!gm BUG in s-coordinate this does not work! 
     
    365367      INTEGER  ::   ji, jj, jl, jk            ! dummy loop indices 
    366368      INTEGER  ::   inum                      ! temporary logical unit 
     369      INTEGER  ::   ierror                    ! error flag 
    367370      INTEGER  ::   ii_bump, ij_bump, ih      ! bump center position 
    368371      INTEGER  ::   ii0, ii1, ij0, ij1, ik    ! local indices 
    369372      REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics  
    370373      REAL(wp) ::   zi, zj, zh, zhmin         ! local scalars 
    371       INTEGER , POINTER, DIMENSION(:,:) ::   idta   ! global domain integer data 
    372       REAL(wp), POINTER, DIMENSION(:,:) ::   zdta   ! global domain scalar data 
     374      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   idta   ! global domain integer data 
     375      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdta   ! global domain scalar data 
    373376      !!---------------------------------------------------------------------- 
    374377      ! 
    375378      IF( nn_timing == 1 )  CALL timing_start('zgr_bat') 
    376       ! 
    377       CALL wrk_alloc( jpidta, jpjdta, idta ) 
    378       CALL wrk_alloc( jpidta, jpjdta, zdta ) 
    379379      ! 
    380380      IF(lwp) WRITE(numout,*) 
     
    385385         !                                            ! ================== ! 
    386386         !                                            ! global domain level and meter bathymetry (idta,zdta) 
     387         ! 
     388         ALLOCATE( idta(jpidta,jpjdta), STAT=ierror ) 
     389         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' ) 
     390         ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror ) 
     391         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' ) 
    387392         ! 
    388393         IF( ntopo == 0 ) THEN                        ! flat basin 
     
    468473         misfdep(:,:)=1 
    469474         ! 
    470          ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 
    471          IF( cp_cfg == "isomip" ) THEN  
    472            !  
    473            risfdep(:,:)=200.e0  
    474            misfdep(:,:)=1  
    475            ij0 = 1 ; ij1 = 40  
    476            DO jj = mj0(ij0), mj1(ij1)  
    477               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp  
    478                 END DO  
    479             WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp  
    480            !  
    481          ELSEIF ( cp_cfg == "isomip2" ) THEN 
    482          !  
    483             risfdep(:,:)=0.e0 
    484             misfdep(:,:)=1 
    485             ij0 = 1 ; ij1 = 40 
    486             DO jj = mj0(ij0), mj1(ij1) 
    487                risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
    488             END DO 
    489             WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp 
    490          END IF 
     475         DEALLOCATE( idta, zdta ) 
    491476         ! 
    492477         !                                            ! ================ ! 
     
    529514         IF( ln_zps .OR. ln_sco )   THEN              ! zps or sco : read meter bathymetry 
    530515            CALL iom_open ( 'bathy_meter.nc', inum )  
    531             CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy ) 
     516            IF ( ln_isfcav ) THEN 
     517               CALL iom_get  ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. ) 
     518            ELSE 
     519               CALL iom_get  ( inum, jpdom_data, 'Bathymetry'    , bathy, lrowattr=ln_use_jattr  ) 
     520            END IF 
    532521            CALL iom_close( inum ) 
    533             !   
     522            !                                                 
    534523            risfdep(:,:)=0._wp          
    535524            misfdep(:,:)=1              
     
    579568      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==! 
    580569         ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin 
    581          WHERE (bathy == risfdep) 
    582             bathy   = 0.0_wp ; risfdep = 0.0_wp 
    583          END WHERE 
     570         IF ( ln_isfcav ) THEN 
     571            WHERE (bathy == risfdep) 
     572               bathy   = 0.0_wp ; risfdep = 0.0_wp 
     573            END WHERE 
     574         END IF 
    584575         ! end patch 
    585576         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level 
     
    592583         IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik 
    593584      ENDIF 
    594       ! 
    595       CALL wrk_dealloc( jpidta, jpjdta, idta ) 
    596       CALL wrk_dealloc( jpidta, jpjdta, zdta ) 
    597585      ! 
    598586      IF( nn_timing == 1 )  CALL timing_stop('zgr_bat') 
     
    959947      !!---------------------------------------------------------------------- 
    960948      !! 
     949      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
     950      INTEGER  ::   ik, it, ikb, ikt ! temporary integers 
     951      LOGICAL  ::   ll_print         ! Allow  control print for debugging 
     952      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
     953      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
     954      REAL(wp) ::   zmax             ! Maximum depth 
     955      REAL(wp) ::   zdiff            ! temporary scalar 
     956      REAL(wp) ::   zrefdep          ! temporary scalar 
     957      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt 
     958      !!--------------------------------------------------------------------- 
     959      ! 
     960      IF( nn_timing == 1 )  CALL timing_start('zgr_zps') 
     961      ! 
     962      CALL wrk_alloc( jpi, jpj, jpk, zprt ) 
     963      ! 
     964      IF(lwp) WRITE(numout,*) 
     965      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps' 
     966      IF(lwp) WRITE(numout,*) '    ~~~~~~~ ' 
     967      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used' 
     968 
     969      ll_print = .FALSE.                   ! Local variable for debugging 
     970       
     971      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth 
     972         WRITE(numout,*) 
     973         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)' 
     974         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 
     975      ENDIF 
     976 
     977 
     978      ! bathymetry in level (from bathy_meter) 
     979      ! =================== 
     980      zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 
     981      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat) 
     982      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0 
     983      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level 
     984      END WHERE 
     985 
     986      ! Compute mbathy for ocean points (i.e. the number of ocean levels) 
     987      ! find the number of ocean levels such that the last level thickness 
     988      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 
     989      ! e3t_1d is the reference level thickness 
     990      DO jk = jpkm1, 1, -1 
     991         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     992         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
     993      END DO 
     994 
     995      IF ( ln_isfcav ) CALL zgr_isf 
     996 
     997      ! Scale factors and depth at T- and W-points 
     998      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate 
     999         gdept_0(:,:,jk) = gdept_1d(jk) 
     1000         gdepw_0(:,:,jk) = gdepw_1d(jk) 
     1001         e3t_0  (:,:,jk) = e3t_1d  (jk) 
     1002         e3w_0  (:,:,jk) = e3w_1d  (jk) 
     1003      END DO 
     1004      !  
     1005      DO jj = 1, jpj 
     1006         DO ji = 1, jpi 
     1007            ik = mbathy(ji,jj) 
     1008            IF( ik > 0 ) THEN               ! ocean point only 
     1009               ! max ocean level case 
     1010               IF( ik == jpkm1 ) THEN 
     1011                  zdepwp = bathy(ji,jj) 
     1012                  ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
     1013                  ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
     1014                  e3t_0(ji,jj,ik  ) = ze3tp 
     1015                  e3t_0(ji,jj,ik+1) = ze3tp 
     1016                  e3w_0(ji,jj,ik  ) = ze3wp 
     1017                  e3w_0(ji,jj,ik+1) = ze3tp 
     1018                  gdepw_0(ji,jj,ik+1) = zdepwp 
     1019                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
     1020                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
     1021                  ! 
     1022               ELSE                         ! standard case 
     1023                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
     1024                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
     1025                  ENDIF 
     1026!gm Bug?  check the gdepw_1d 
     1027                  !       ... on ik 
     1028                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
     1029                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
     1030                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
     1031                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
     1032                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
     1033                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
     1034                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
     1035                  !       ... on ik+1 
     1036                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1037                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1038                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
     1039               ENDIF 
     1040            ENDIF 
     1041         END DO 
     1042      END DO 
     1043      ! 
     1044      it = 0 
     1045      DO jj = 1, jpj 
     1046         DO ji = 1, jpi 
     1047            ik = mbathy(ji,jj) 
     1048            IF( ik > 0 ) THEN               ! ocean point only 
     1049               e3tp (ji,jj) = e3t_0(ji,jj,ik) 
     1050               e3wp (ji,jj) = e3w_0(ji,jj,ik) 
     1051               ! test 
     1052               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  ) 
     1053               IF( zdiff <= 0._wp .AND. lwp ) THEN  
     1054                  it = it + 1 
     1055                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
     1056                  WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
     1057                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
     1058                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  ) 
     1059               ENDIF 
     1060            ENDIF 
     1061         END DO 
     1062      END DO 
     1063      ! 
     1064      IF ( ln_isfcav ) THEN 
     1065      ! (ISF) Definition of e3t, u, v, w for ISF case 
     1066         DO jj = 1, jpj  
     1067            DO ji = 1, jpi  
     1068               ik = misfdep(ji,jj)  
     1069               IF( ik > 1 ) THEN               ! ice shelf point only  
     1070                  IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
     1071                  gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
     1072!gm Bug?  check the gdepw_0  
     1073               !       ... on ik  
     1074                  gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
     1075                     &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
     1076                     &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
     1077                  e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
     1078                  e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
     1079 
     1080                  IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
     1081                     e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
     1082                  ENDIF  
     1083               !       ... on ik / ik-1  
     1084                  e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
     1085                  e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
     1086! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
     1087                  gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
     1088               ENDIF  
     1089            END DO  
     1090         END DO  
     1091      !  
     1092         it = 0  
     1093         DO jj = 1, jpj  
     1094            DO ji = 1, jpi  
     1095               ik = misfdep(ji,jj)  
     1096               IF( ik > 1 ) THEN               ! ice shelf point only  
     1097                  e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
     1098                  e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
     1099               ! test  
     1100                  zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
     1101                  IF( zdiff <= 0. .AND. lwp ) THEN   
     1102                     it = it + 1  
     1103                     WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
     1104                     WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
     1105                     WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
     1106                     WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
     1107                  ENDIF  
     1108               ENDIF  
     1109            END DO  
     1110         END DO  
     1111      END IF 
     1112      ! END (ISF) 
     1113 
     1114      ! Scale factors and depth at U-, V-, UW and VW-points 
     1115      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     1116         e3u_0 (:,:,jk) = e3t_1d(jk) 
     1117         e3v_0 (:,:,jk) = e3t_1d(jk) 
     1118         e3uw_0(:,:,jk) = e3w_1d(jk) 
     1119         e3vw_0(:,:,jk) = e3w_1d(jk) 
     1120      END DO 
     1121      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
     1122         DO jj = 1, jpjm1 
     1123            DO ji = 1, fs_jpim1   ! vector opt. 
     1124               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
     1125               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
     1126               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
     1127               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
     1128            END DO 
     1129         END DO 
     1130      END DO 
     1131      IF ( ln_isfcav ) THEN 
     1132      ! (ISF) define e3uw (adapted for 2 cells in the water column) 
     1133         DO jj = 2, jpjm1  
     1134            DO ji = 2, fs_jpim1   ! vector opt.  
     1135               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 
     1136               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 
     1137               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) & 
     1138                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) ) 
     1139               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 
     1140               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 
     1141               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) & 
     1142                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) ) 
     1143            END DO 
     1144         END DO 
     1145      END IF 
     1146 
     1147      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
     1148      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
     1149      ! 
     1150      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1151         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
     1152         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk) 
     1153         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk) 
     1154         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk) 
     1155      END DO 
     1156       
     1157      ! Scale factor at F-point 
     1158      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     1159         e3f_0(:,:,jk) = e3t_1d(jk) 
     1160      END DO 
     1161      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
     1162         DO jj = 1, jpjm1 
     1163            DO ji = 1, fs_jpim1   ! vector opt. 
     1164               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
     1165            END DO 
     1166         END DO 
     1167      END DO 
     1168      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions 
     1169      ! 
     1170      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1171         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk) 
     1172      END DO 
     1173!!gm  bug ? :  must be a do loop with mj0,mj1 
     1174      !  
     1175      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
     1176      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)  
     1177      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)  
     1178      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)  
     1179      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)  
     1180 
     1181      ! Control of the sign 
     1182      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
     1183      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
     1184      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
     1185      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
     1186      
     1187      ! Compute gdep3w_0 (vertical sum of e3w) 
     1188      IF ( ln_isfcav ) THEN ! if cavity 
     1189         WHERE (misfdep == 0) misfdep = 1 
     1190         DO jj = 1,jpj 
     1191            DO ji = 1,jpi 
     1192               gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
     1193               DO jk = 2, misfdep(ji,jj) 
     1194                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1195               END DO 
     1196               IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 
     1197               DO jk = misfdep(ji,jj) + 1, jpk 
     1198                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1199               END DO 
     1200            END DO 
     1201         END DO 
     1202      ELSE ! no cavity 
     1203         gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 
     1204         DO jk = 2, jpk 
     1205            gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
     1206         END DO 
     1207      END IF 
     1208      !                                               ! ================= ! 
     1209      IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
     1210         !                                            ! ================= ! 
     1211         DO jj = 1,jpj 
     1212            DO ji = 1, jpi 
     1213               ik = MAX( mbathy(ji,jj), 1 ) 
     1214               zprt(ji,jj,1) = e3t_0   (ji,jj,ik) 
     1215               zprt(ji,jj,2) = e3w_0   (ji,jj,ik) 
     1216               zprt(ji,jj,3) = e3u_0   (ji,jj,ik) 
     1217               zprt(ji,jj,4) = e3v_0   (ji,jj,ik) 
     1218               zprt(ji,jj,5) = e3f_0   (ji,jj,ik) 
     1219               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 
     1220            END DO 
     1221         END DO 
     1222         WRITE(numout,*) 
     1223         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1224         WRITE(numout,*) 
     1225         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1226         WRITE(numout,*) 
     1227         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1228         WRITE(numout,*) 
     1229         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1230         WRITE(numout,*) 
     1231         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1232         WRITE(numout,*) 
     1233         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1234      ENDIF   
     1235      ! 
     1236      CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 
     1237      ! 
     1238      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps') 
     1239      ! 
     1240   END SUBROUTINE zgr_zps 
     1241 
     1242   SUBROUTINE zgr_isf 
     1243      !!---------------------------------------------------------------------- 
     1244      !!                    ***  ROUTINE zgr_isf  *** 
     1245      !!    
     1246      !! ** Purpose :   check the bathymetry in levels 
     1247      !!    
     1248      !! ** Method  :   THe water column have to contained at least 2 cells 
     1249      !!                Bathymetry and isfdraft are modified (dig/close) to respect 
     1250      !!                this criterion. 
     1251      !!                  
     1252      !!    
     1253      !! ** Action  : - test compatibility between isfdraft and bathy  
     1254      !!              - bathy and isfdraft are modified 
     1255      !!---------------------------------------------------------------------- 
     1256      !!    
    9611257      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    9621258      INTEGER  ::   ik, it           ! temporary integers 
     
    9691265      REAL(wp) ::   zdiff            ! temporary scalar 
    9701266      REAL(wp) ::   zrefdep          ! temporary scalar 
    971       REAL(wp) ::   zbathydiff, zrisfdepdiff  
    972       REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 3D workspace (ISH) 
    973       INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep   ! 3D workspace (ISH) 
    974       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt 
     1267      REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar 
     1268      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH) 
     1269      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH) 
    9751270      !!--------------------------------------------------------------------- 
    9761271      ! 
    977       IF( nn_timing == 1 )  CALL timing_start('zgr_zps') 
    978       ! 
    979       CALL wrk_alloc( jpi, jpj, jpk, zprt ) 
     1272      IF( nn_timing == 1 )  CALL timing_start('zgr_isf') 
     1273      ! 
    9801274      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 
    981       CALL wrk_alloc( jpi, jpj, zmbathy, zmisfdep) 
    982       ! 
    983       IF(lwp) WRITE(numout,*) 
    984       IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps' 
    985       IF(lwp) WRITE(numout,*) '    ~~~~~~~ ' 
    986       IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used' 
    987  
    988       ll_print = .FALSE.                   ! Local variable for debugging 
    989        
    990       IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth 
    991          WRITE(numout,*) 
    992          WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)' 
    993          CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 
    994       ENDIF 
    995  
    996       ! bathymetry in level (from bathy_meter) 
    997       ! =================== 
    998       zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 
    999       bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat) 
    1000       WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0 
    1001       ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level 
    1002       END WHERE 
    1003  
    1004       ! Compute mbathy for ocean points (i.e. the number of ocean levels) 
    1005       ! find the number of ocean levels such that the last level thickness 
    1006       ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 
    1007       ! e3t_1d is the reference level thickness 
    1008       DO jk = jpkm1, 1, -1 
    1009          zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
    1010          WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
    1011       END DO 
     1275      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy ) 
     1276 
     1277 
    10121278      ! (ISF) compute misfdep 
    10131279      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1   
     
    10531319            misfdep(jpi,:) = misfdep(  2  ,:)  
    10541320         ENDIF 
    1055   
     1321 
    10561322         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
    10571323            mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west 
    10581324            mbathy(jpi,:) = mbathy(  2  ,:) 
    10591325         ENDIF 
    1060   
     1326 
    10611327         ! split last cell if possible (only where water column is 2 cell or less) 
    10621328         DO jk = jpkm1, 1, -1 
     
    10761342            END WHERE 
    10771343         END DO 
    1078   
     1344 
    10791345  
    10801346 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 
     
    12521518  
    12531519 ! remove single point "bay" on isf coast line in the ice shelf draft' 
    1254          DO jk = 1, jpk 
     1520         DO jk = 2, jpk 
    12551521            WHERE (misfdep==0) misfdep=jpk 
    12561522            zmask=0 
     
    13571623               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0 
    13581624               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    1359                IF( ibtest == 0 ) THEN 
     1625               IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN 
    13601626                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 
    13611627               END IF 
     
    14731739      ENDIF  
    14741740 
    1475       ! Scale factors and depth at T- and W-points 
    1476       DO jk = 1, jpk                        ! intitialization to the reference z-coordinate 
    1477          gdept_0(:,:,jk) = gdept_1d(jk) 
    1478          gdepw_0(:,:,jk) = gdepw_1d(jk) 
    1479          e3t_0  (:,:,jk) = e3t_1d  (jk) 
    1480          e3w_0  (:,:,jk) = e3w_1d  (jk) 
    1481       END DO 
    1482       !  
    1483       DO jj = 1, jpj 
    1484          DO ji = 1, jpi 
    1485             ik = mbathy(ji,jj) 
    1486             IF( ik > 0 ) THEN               ! ocean point only 
    1487                ! max ocean level case 
    1488                IF( ik == jpkm1 ) THEN 
    1489                   zdepwp = bathy(ji,jj) 
    1490                   ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
    1491                   ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
    1492                   e3t_0(ji,jj,ik  ) = ze3tp 
    1493                   e3t_0(ji,jj,ik+1) = ze3tp 
    1494                   e3w_0(ji,jj,ik  ) = ze3wp 
    1495                   e3w_0(ji,jj,ik+1) = ze3tp 
    1496                   gdepw_0(ji,jj,ik+1) = zdepwp 
    1497                   gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
    1498                   gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
    1499                   ! 
    1500                ELSE                         ! standard case 
    1501                   IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
    1502                   ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
    1503                   ENDIF 
    1504 !gm Bug?  check the gdepw_1d 
    1505                   !       ... on ik 
    1506                   gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
    1507                      &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
    1508                      &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
    1509                   e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
    1510                      &                          / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
    1511                   e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
    1512                      &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
    1513                   !       ... on ik+1 
    1514                   e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1515                   e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1516                   gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
    1517                ENDIF 
    1518             ENDIF 
    1519          END DO 
    1520       END DO 
    1521       ! 
    1522       it = 0 
    1523       DO jj = 1, jpj 
    1524          DO ji = 1, jpi 
    1525             ik = mbathy(ji,jj) 
    1526             IF( ik > 0 ) THEN               ! ocean point only 
    1527                e3tp (ji,jj) = e3t_0(ji,jj,ik) 
    1528                e3wp (ji,jj) = e3w_0(ji,jj,ik) 
    1529                ! test 
    1530                zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  ) 
    1531                IF( zdiff <= 0._wp .AND. lwp ) THEN  
    1532                   it = it + 1 
    1533                   WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
    1534                   WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
    1535                   WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
    1536                   WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  ) 
    1537                ENDIF 
    1538             ENDIF 
    1539          END DO 
    1540       END DO 
    1541       ! 
    1542       ! (ISF) Definition of e3t, u, v, w for ISF case 
    1543       DO jj = 1, jpj  
    1544          DO ji = 1, jpi  
    1545             ik = misfdep(ji,jj)  
    1546             IF( ik > 1 ) THEN               ! ice shelf point only  
    1547                IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
    1548                gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
    1549 !gm Bug?  check the gdepw_0  
    1550                !       ... on ik  
    1551                gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
    1552                   &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
    1553                   &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
    1554                e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
    1555                e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
    1556  
    1557                IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
    1558                   e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
    1559                 ENDIF  
    1560                !       ... on ik / ik-1  
    1561                e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
    1562                e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
    1563 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
    1564                gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
    1565             ENDIF  
    1566          END DO  
    1567       END DO  
    1568       !  
    1569       it = 0  
    1570       DO jj = 1, jpj  
    1571          DO ji = 1, jpi  
    1572             ik = misfdep(ji,jj)  
    1573             IF( ik > 1 ) THEN               ! ice shelf point only  
    1574                e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
    1575                e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
    1576                ! test  
    1577                zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
    1578                IF( zdiff <= 0. .AND. lwp ) THEN   
    1579                   it = it + 1  
    1580                   WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
    1581                   WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
    1582                   WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
    1583                   WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
    1584                ENDIF  
    1585             ENDIF  
    1586          END DO  
    1587       END DO  
    1588       ! END (ISF) 
    1589  
    1590       ! Scale factors and depth at U-, V-, UW and VW-points 
    1591       DO jk = 1, jpk                        ! initialisation to z-scale factors 
    1592          e3u_0 (:,:,jk) = e3t_1d(jk) 
    1593          e3v_0 (:,:,jk) = e3t_1d(jk) 
    1594          e3uw_0(:,:,jk) = e3w_1d(jk) 
    1595          e3vw_0(:,:,jk) = e3w_1d(jk) 
    1596       END DO 
    1597       DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
    1598          DO jj = 1, jpjm1 
    1599             DO ji = 1, fs_jpim1   ! vector opt. 
    1600                e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
    1601                e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
    1602                e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
    1603                e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
    1604             END DO 
    1605          END DO 
    1606       END DO 
    1607       ! (ISF) define e3uw 
    1608       DO jk = 2,jpk                           
    1609          DO jj = 1, jpjm1  
    1610             DO ji = 1, fs_jpim1   ! vector opt.  
    1611                e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj  ,jk) ) & 
    1612                  &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj  ,jk-1) ) 
    1613                e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji  ,jj+1,jk) ) & 
    1614                  &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji  ,jj+1,jk-1) ) 
    1615             END DO  
    1616          END DO  
    1617       END DO 
    1618       !End (ISF) 
    1619        
    1620       CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
    1621       CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    1622       ! 
    1623       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1624          WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
    1625          WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk) 
    1626          WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk) 
    1627          WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk) 
    1628       END DO 
    1629        
    1630       ! Scale factor at F-point 
    1631       DO jk = 1, jpk                        ! initialisation to z-scale factors 
    1632          e3f_0(:,:,jk) = e3t_1d(jk) 
    1633       END DO 
    1634       DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
    1635          DO jj = 1, jpjm1 
    1636             DO ji = 1, fs_jpim1   ! vector opt. 
    1637                e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
    1638             END DO 
    1639          END DO 
    1640       END DO 
    1641       CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions 
    1642       ! 
    1643       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1644          WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk) 
    1645       END DO 
    1646 !!gm  bug ? :  must be a do loop with mj0,mj1 
    1647       !  
    1648       e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
    1649       e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)  
    1650       e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)  
    1651       e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)  
    1652       e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)  
    1653  
    1654       ! Control of the sign 
    1655       IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
    1656       IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
    1657       IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
    1658       IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
    1659       
    1660       ! Compute gdep3w_0 (vertical sum of e3w) 
    1661       WHERE (misfdep == 0) misfdep = 1 
    1662       DO jj = 1,jpj 
    1663          DO ji = 1,jpi 
    1664             gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
    1665             DO jk = 2, misfdep(ji,jj) 
    1666                gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    1667             END DO 
    1668             IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 
    1669             DO jk = misfdep(ji,jj) + 1, jpk 
    1670                gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    1671             END DO 
    1672         END DO 
    1673       END DO 
    1674       !                                               ! ================= ! 
    1675       IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
    1676          !                                            ! ================= ! 
    1677          DO jj = 1,jpj 
    1678             DO ji = 1, jpi 
    1679                ik = MAX( mbathy(ji,jj), 1 ) 
    1680                zprt(ji,jj,1) = e3t_0   (ji,jj,ik) 
    1681                zprt(ji,jj,2) = e3w_0   (ji,jj,ik) 
    1682                zprt(ji,jj,3) = e3u_0   (ji,jj,ik) 
    1683                zprt(ji,jj,4) = e3v_0   (ji,jj,ik) 
    1684                zprt(ji,jj,5) = e3f_0   (ji,jj,ik) 
    1685                zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 
    1686             END DO 
    1687          END DO 
    1688          WRITE(numout,*) 
    1689          WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1690          WRITE(numout,*) 
    1691          WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1692          WRITE(numout,*) 
    1693          WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1694          WRITE(numout,*) 
    1695          WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1696          WRITE(numout,*) 
    1697          WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1698          WRITE(numout,*) 
    1699          WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1700       ENDIF   
    1701       ! 
    1702       CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 
    17031741      CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 
    17041742      CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 
    1705       ! 
    1706       IF( nn_timing == 1 )  CALL timing_stop('zgr_zps') 
    1707       ! 
    1708    END SUBROUTINE zgr_zps 
     1743 
     1744      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf') 
     1745       
     1746   END SUBROUTINE 
    17091747 
    17101748   SUBROUTINE zgr_sco 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90

    • Property svn:keywords set to Id
    r5038 r5620  
    3939   !!---------------------------------------------------------------------- 
    4040   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
    41    !! $Id: dtatem.F90 2392 2010-11-15 21:20:05Z gm $  
     41   !! $Id$  
    4242   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4343   !!---------------------------------------------------------------------- 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r5038 r5620  
    6969      !! ** Purpose :   Initialization of the dynamics and tracer fields. 
    7070      !!---------------------------------------------------------------------- 
    71       ! - ML - needed for initialization of e3t_b 
    72       INTEGER  ::  ji,jj,jk     ! dummy loop indices 
    73       REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::  zuvd    ! U & V data workspace 
     71      INTEGER ::   ji, jj, jk   ! dummy loop indices 
     72      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zuvd    ! U & V data workspace 
    7473      !!---------------------------------------------------------------------- 
    7574      ! 
     
    8483      IF( lk_c1d ) CALL dta_uvd_init          ! Initialization of U & V input data 
    8584 
    86       rhd  (:,:,:  ) = 0._wp 
    87       rhop (:,:,:  ) = 0._wp 
    88       rn2  (:,:,:  ) = 0._wp 
    89       tsa  (:,:,:,:) = 0._wp    
    90       rab_b(:,:,:,:) = 0._wp 
    91       rab_n(:,:,:,:) = 0._wp 
     85      rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk 
     86      rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk 
     87      tsa  (:,:,:,:) = 0._wp                                   ! set one for all to 0 at level jpk 
     88      rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk 
    9289 
    9390      IF( ln_rstart ) THEN                    ! Restart from a file 
     
    113110         ELSEIF( cp_cfg == 'gyre' ) THEN          
    114111            CALL istate_gyre                     ! GYRE  configuration : start from pre-defined T-S fields 
    115         ELSEIF( cp_cfg == 'isomip' .OR. cp_cfg == 'isomip2') THEN 
    116             IF(lwp) WRITE(numout,*) 'Initialization of T+S for ISOMIP domain'  
    117             tsn(:,:,:,jp_tem)=-1.9*tmask(:,:,:)          ! ISOMIP configuration : start from constant T+S fields  
    118             tsn(:,:,:,jp_sal)=34.4*tmask(:,:,:) 
    119             tsb(:,:,:,:)=tsn(:,:,:,:)   
    120112         ELSE                                    ! Initial T-S, U-V fields read in files 
    121113            IF ( ln_tsd_init ) THEN              ! read 3D T and S data at nit000 
     
    137129         CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )        ! before potential and in situ densities 
    138130#if ! defined key_c1d 
    139          IF( ln_zps )    CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
    140             &                                      rhd, gru , grv, aru, arv, gzu, gzv, ge3ru, ge3rv,  &             ! 
    141             &                                      gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
     131         IF( ln_zps .AND. .NOT. ln_isfcav)                                 & 
     132            &            CALL zps_hde    ( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
     133            &                                            rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
     134         IF( ln_zps .AND.       ln_isfcav)                                 & 
     135            &            CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
     136            &                                            rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     137            &                                     gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
    142138#endif 
    143139         !    
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90

    r5038 r5620  
    4141   REAL(wp), PUBLIC ::   rt0      = 273.15_wp        !: freezing point of fresh water [Kelvin] 
    4242#if defined key_lim3 
    43    REAL(wp), PUBLIC ::   rt0_snow = 273.16_wp        !: melting point of snow         [Kelvin] 
    44    REAL(wp), PUBLIC ::   rt0_ice  = 273.16_wp        !: melting point of ice          [Kelvin] 
     43   REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp        !: melting point of snow         [Kelvin] 
     44   REAL(wp), PUBLIC ::   rt0_ice  = 273.15_wp        !: melting point of ice          [Kelvin] 
    4545#else 
    4646   REAL(wp), PUBLIC ::   rt0_snow = 273.15_wp        !: melting point of snow         [Kelvin] 
     
    5151   REAL(wp), PUBLIC ::   rcp                         !: ocean specific heat           [J/Kelvin] 
    5252   REAL(wp), PUBLIC ::   r1_rcp                      !: = 1. / rcp                    [Kelvin/J] 
     53   REAL(wp), PUBLIC ::   rau0_rcp                    !: = rau0 * rcp  
    5354   REAL(wp), PUBLIC ::   r1_rau0_rcp                 !: = 1. / ( rau0 * rcp ) 
    5455 
     
    8283   REAL(wp), PUBLIC ::   xlic     =  300.33e+6_wp    !: volumetric latent heat fusion of ice                  [J/m3] 
    8384   REAL(wp), PUBLIC ::   xsn      =    2.8e+6_wp     !: volumetric latent heat of sublimation of snow         [J/m3] 
     85#endif 
     86#if defined key_lim3 
     87   REAL(wp), PUBLIC ::   r1_rhoic                    !: 1 / rhoic 
     88   REAL(wp), PUBLIC ::   r1_rhosn                    !: 1 / rhosn 
    8489#endif 
    8590   !!---------------------------------------------------------------------- 
     
    166171      lfus = xlsn / rhosn        ! latent heat of fusion of fresh ice 
    167172#endif 
    168  
     173#if defined key_lim3 
     174      r1_rhoic = 1._wp / rhoic 
     175      r1_rhosn = 1._wp / rhosn 
     176#endif 
    169177      IF(lwp) THEN 
    170178         WRITE(numout,*) 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90

    r5038 r5620  
    1717   !!            3.3  ! 2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module 
    1818   !!             -   ! 2010-10  (R. Furner, G. Madec) runoff and cla added directly here 
     19   !!            3.6  ! 2014-11  (P. Mathiot)          isf            added directly here 
    1920   !!---------------------------------------------------------------------- 
    2021 
     
    9798      ! 
    9899      CALL wrk_alloc( jpi  , jpj+2, zwu               ) 
    99       CALL wrk_alloc( jpi+4, jpj  , zwv, kjstart = -1 ) 
     100      CALL wrk_alloc( jpi+4, jpj  , zwv, kistart = -1 ) 
    100101      ! 
    101102      IF( kt == nit000 ) THEN 
     
    236237      ! 
    237238      CALL wrk_dealloc( jpi  , jpj+2, zwu               ) 
    238       CALL wrk_dealloc( jpi+4, jpj  , zwv, kjstart = -1 ) 
     239      CALL wrk_dealloc( jpi+4, jpj  , zwv, kistart = -1 ) 
    239240      ! 
    240241      IF( nn_timing == 1 )  CALL timing_stop('div_cur') 
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90

    r5038 r5620  
    55   !!============================================================================== 
    66   !! History :  1.0  !  2006-11  (G. Madec)  Original code 
    7    !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
     7   !!            3.3  !  2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase 
     8   !!            3.6  !  2015-05  (N. Ducousso, G. Madec)  add Hollingsworth scheme as an option  
    89   !!---------------------------------------------------------------------- 
    910 
     
    1718   USE dynkeg          ! kinetic energy gradient          (dyn_keg      routine) 
    1819   USE dynzad          ! vertical advection               (dyn_zad      routine) 
     20   ! 
    1921   USE in_out_manager  ! I/O manager 
    2022   USE lib_mpp         ! MPP library 
     
    2527 
    2628   PUBLIC dyn_adv       ! routine called by step module 
    27    PUBLIC dyn_adv_init  ! routine called by opa module 
     29   PUBLIC dyn_adv_init  ! routine called by opa  module 
    2830  
     31   !                                    !* namdyn_adv namelist * 
    2932   LOGICAL, PUBLIC ::   ln_dynadv_vec   !: vector form flag 
     33   INTEGER, PUBLIC ::   nn_dynkeg       !: scheme of kinetic energy gradient: =0 C2 ; =1 Hollingsworth 
    3034   LOGICAL, PUBLIC ::   ln_dynadv_cen2  !: flux form - 2nd order centered scheme flag 
    3135   LOGICAL, PUBLIC ::   ln_dynadv_ubs   !: flux form - 3rd order UBS scheme flag 
     
    3842#  include "vectopt_loop_substitute.h90" 
    3943   !!---------------------------------------------------------------------- 
    40    !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     44   !! NEMO/OPA 3.6 , NEMO Consortium (2015) 
    4145   !! $Id$ 
    4246   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    6367      SELECT CASE ( nadv )                  ! compute advection trend and add it to general trend 
    6468      CASE ( 0 )      
    65                       CALL dyn_keg     ( kt )    ! vector form : horizontal gradient of kinetic energy 
    66                       CALL dyn_zad     ( kt )    ! vector form : vertical advection 
     69                      CALL dyn_keg     ( kt, nn_dynkeg )    ! vector form : horizontal gradient of kinetic energy 
     70                      CALL dyn_zad     ( kt )               ! vector form : vertical advection 
    6771      CASE ( 1 )      
    68                       CALL dyn_keg     ( kt )    ! vector form : horizontal gradient of kinetic energy 
    69                       CALL dyn_zad_zts ( kt )    ! vector form : vertical advection with sub-timestepping 
     72                      CALL dyn_keg     ( kt, nn_dynkeg )    ! vector form : horizontal gradient of kinetic energy 
     73                      CALL dyn_zad_zts ( kt )               ! vector form : vertical advection with sub-timestepping 
    7074      CASE ( 2 )  
    71                       CALL dyn_adv_cen2( kt )    ! 2nd order centered scheme 
     75                      CALL dyn_adv_cen2( kt )               ! 2nd order centered scheme 
    7276      CASE ( 3 )    
    73                       CALL dyn_adv_ubs ( kt )    ! 3rd order UBS      scheme 
     77                      CALL dyn_adv_ubs ( kt )               ! 3rd order UBS      scheme 
    7478      ! 
    75       CASE (-1 )                                 ! esopa: test all possibility with control print 
    76                       CALL dyn_keg     ( kt ) 
     79      CASE (-1 )                                            ! esopa: test all possibility with control print 
     80                      CALL dyn_keg     ( kt, nn_dynkeg ) 
    7781                      CALL dyn_zad     ( kt ) 
    7882                      CALL dyn_adv_cen2( kt ) 
     
    9296      !!              momentum advection formulation & scheme and set nadv 
    9397      !!---------------------------------------------------------------------- 
    94       INTEGER ::   ioptio 
    95       INTEGER  ::   ios                 ! Local integer output status for namelist read 
    96       !! 
    97       NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts 
     98      INTEGER ::   ioptio, ios   ! Local integer 
     99      ! 
     100      NAMELIST/namdyn_adv/ ln_dynadv_vec, nn_dynkeg, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts 
    98101      !!---------------------------------------------------------------------- 
    99  
     102      ! 
    100103      REWIND( numnam_ref )              ! Namelist namdyn_adv in reference namelist : Momentum advection scheme 
    101104      READ  ( numnam_ref, namdyn_adv, IOSTAT = ios, ERR = 901) 
     
    112115         WRITE(numout,*) '~~~~~~~~~~~' 
    113116         WRITE(numout,*) '       Namelist namdyn_adv : chose a advection formulation & scheme for momentum' 
    114          WRITE(numout,*) '          Vector/flux form (T/F)             ln_dynadv_vec  = ', ln_dynadv_vec 
    115          WRITE(numout,*) '          2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2 
    116          WRITE(numout,*) '          3rd order UBS advection scheme     ln_dynadv_ubs  = ', ln_dynadv_ubs 
    117          WRITE(numout,*) '      Sub timestepping of vertical advection ln_dynzad_zts  = ', ln_dynzad_zts 
     117         WRITE(numout,*) '          Vector/flux form (T/F)                           ln_dynadv_vec  = ', ln_dynadv_vec 
     118         WRITE(numout,*) '          = 0 standard scheme  ; =1 Hollingsworth scheme   nn_dynkeg      = ', nn_dynkeg 
     119         WRITE(numout,*) '          2nd order centred advection scheme               ln_dynadv_cen2 = ', ln_dynadv_cen2 
     120         WRITE(numout,*) '          3rd order UBS advection scheme                   ln_dynadv_ubs  = ', ln_dynadv_ubs 
     121         WRITE(numout,*) '          Sub timestepping of vertical advection           ln_dynzad_zts  = ', ln_dynzad_zts 
    118122      ENDIF 
    119123 
     
    126130      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' ) 
    127131      IF( ln_dynzad_zts .AND. .NOT. ln_dynadv_vec )   & 
    128           CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' ) 
     132         CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' ) 
     133      IF( nn_dynkeg /= nkeg_C2 .AND. nn_dynkeg /= nkeg_HW )   &   
     134         CALL ctl_stop( 'KEG scheme wrong value of nn_dynkeg' ) 
    129135 
    130136      !                               ! Set nadv 
     
    137143      IF(lwp) THEN                    ! Print the choice 
    138