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

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

Changeset 5120


Ignore:
Timestamp:
2015-03-03T17:11:55+01:00 (10 years ago)
Author:
acc
Message:

#1473. Re-organisation and optimisation of ice shelf cavity option. This commit merges changes from the dev_r5094_UKMO_ISFCLEAN branch onto the trunk. Results will change, even with ln_isfcav=F, due to a return to original definitions of the vertical metrics. All changes have been reviewed and SETTE tested.

Location:
trunk
Files:
56 edited
1 copied

Legend:

Unmodified
Added
Removed

  • trunk/DOC/TexFiles/Biblio/Biblio.bib

    r4560 r5120  
    271271  volume = {326}, 
    272272  pages = {677--684} 
     273} 
     274 
     275@ARTICLE{Beckmann2003, 
     276  author = {A. Beckmann and H. Goosse}, 
     277  title = {A parameterization of ice shelf-ocean interaction for climate models}, 
     278  journal = OM 
     279  year = {2003} 
     280  volume = {5} 
     281  pages = {157--170} 
    273282} 
    274283 

  • trunk/DOC/TexFiles/Chapters/Chap_DOM.tex

    r4147 r5120  
    493493$z(i,j,k,t)$ (Fig.~\ref{Fig_z_zps_s_sps}f). This option can be used with full step  
    494494bathymetry or $s$-coordinate (hybrid and partial step coordinates have not  
    495 yet been tested in NEMO v2.3).  
     495yet been tested in NEMO v2.3). If using $z$-coordinate with partial step bathymetry 
     496(\np{ln\_zps}~=~true), ocean cavity beneath ice shelves can be open (\np{ln\_isfcav}~=~true). 
    496497 
    497498Contrary to the horizontal grid, the vertical grid is computed in the code and no  

  • trunk/DOC/TexFiles/Chapters/Chap_DYN.tex

    r4759 r5120  
    627627\eqref{Eq_dynhpg_zco_surf} - \eqref{Eq_dynhpg_zco}, and $z_T$ is the depth of  
    628628the $T$-point evaluated from the sum of the vertical scale factors at the $w$-point  
    629 ($e_{3w}$).  
     629($e_{3w}$). 
     630  
     631$\bullet$ Traditional coding with adaptation for ice shelf cavities (\np{ln\_dynhpg\_isf}=true). 
     632This scheme need the activation of ice shelf cavities (\np{ln\_isfcav}=true). 
    630633 
    631634$\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}=true) 

  • trunk/DOC/TexFiles/Chapters/Chap_SBC.tex

    r4661 r5120  
    11% ================================================================ 
    2 % Chapter � Surface Boundary Condition (SBC, ICB)  
    3 % ================================================================ 
    4 \chapter{Surface Boundary Condition (SBC, ICB) } 
     2% Chapter � Surface Boundary Condition (SBC, ISF, ICB)  
     3% ================================================================ 
     4\chapter{Surface Boundary Condition (SBC, ISF, ICB) } 
    55\label{SBC} 
    66\minitoc 
    … …  
    4848below ice-covered areas (using observed ice-cover or a sea-ice model)  
    4949(\np{nn\_ice}~=~0,1, 2 or 3); the addition of river runoffs as surface freshwater  
    50 fluxes or lateral inflow (\np{ln\_rnf}~=~true); the addition of a freshwater flux adjustment  
    51 in order to avoid a mean sea-level drift (\np{nn\_fwb}~=~0,~1~or~2); the  
     50fluxes or lateral inflow (\np{ln\_rnf}~=~true); the addition of isf melting as lateral inflow (parameterisation)  
     51(\np{nn\_isf}~=~2 or 3 and \np{ln\_isfcav}~=~false) or as surface flux at the land-ice ocean interface 
     52(\np{nn\_isf}~=~1 or 4 and \np{ln\_isfcav}~=~true);  
     53the addition of a freshwater flux adjustment in order to avoid a mean sea-level drift (\np{nn\_fwb}~=~0,~1~or~2); the  
    5254transformation of the solar radiation (if provided as daily mean) into a diurnal  
    5355cycle (\np{ln\_dm2dc}~=~true); and a neutral drag coefficient can be read from an external wave  
    … …  
    6062Finally, the different options that further modify the fluxes applied to the ocean are discussed. 
    6163One of these is modification by icebergs (see \S\ref{ICB_icebergs}), which act as drifting sources of fresh water. 
     64Another example of modification is that due to the ice shelf melting/freezing (see \S\ref{SBC_isf}),  
     65which provides additional sources of fresh water. 
    6266 
    6367 
    … …  
    686690air temperature, sea-surface temperature, cloud cover and relative humidity. 
    687691Sensible heat and latent heat fluxes are computed by classical 
    688 bulk formulae parameterized according to \citet{Kondo1975}. 
     692bulk formulae parameterised according to \citet{Kondo1975}. 
    689693Details on the bulk formulae used can be found in \citet{Maggiore_al_PCE98} and \citet{Castellari_al_JMS1998}. 
    690694 
    … …  
    826830\Pi-g\delta = (1+k-h) \Pi_{A}(\lambda,\phi) 
    827831\end{equation} 
    828 with $k$ a number of Love estimated to 0.6 which parametrized the astronomical tidal land, 
    829 and $h$ a number of Love to 0.3 which parametrized the parametrization due to the astronomical tidal land. 
     832with $k$ a number of Love estimated to 0.6 which parameterised the astronomical tidal land, 
     833and $h$ a number of Love to 0.3 which parameterised the parameterisation due to the astronomical tidal land. 
    830834 
    831835% ================================================================ 
    … …  
    945949 
    946950%} 
    947  
    948  
     951% ================================================================ 
     952%        Ice shelf melting 
     953% ================================================================ 
     954\section   [Ice shelf melting (\textit{sbcisf})] 
     955                        {Ice shelf melting (\mdl{sbcisf})} 
     956\label{SBC_isf} 
     957%------------------------------------------namsbc_isf---------------------------------------------------- 
     958\namdisplay{namsbc_isf} 
     959%-------------------------------------------------------------------------------------------------------- 
     960Namelist variable in \ngn{namsbc}, \np{nn\_isf},  control the kind of ice shelf representation used.  
     961\begin{description} 
     962\item[\np{nn\_isf}~=~1] 
     963The ice shelf cavity is represented. The fwf and heat flux are computed.  
     964Full description, sensitivity and validation in preparation.  
     965 
     966\item[\np{nn\_isf}~=~2] 
     967A parameterisation of isf is used. The ice shelf cavity is not represented.  
     968The fwf is distributed along the ice shelf edge between the depth of the average grounding line (GL) 
     969(\np{sn\_depmax\_isf}) and the base of the ice shelf along the calving front (\np{sn\_depmin\_isf}) as in (\np{nn\_isf}~=~3).  
     970Furthermore the fwf is computed using the \citet{Beckmann2003} parameterisation of isf melting.  
     971The effective melting length (\np{sn\_Leff\_isf}) is read from a file. 
     972 
     973\item[\np{nn\_isf}~=~3] 
     974A simple parameterisation of isf is used. The ice shelf cavity is not represented.  
     975The fwf (\np{sn\_rnfisf}) is distributed along the ice shelf edge between the depth of the average grounding line (GL) 
     976(\np{sn\_depmax\_isf}) and the base of the ice shelf along the calving front (\np{sn\_depmin\_isf}). 
     977Full description, sensitivity and validation in preparation. 
     978 
     979\item[\np{nn\_isf}~=~4] 
     980The ice shelf cavity is represented. However, the fwf (\np{sn\_fwfisf}) and heat flux (\np{sn\_qisf}) are  
     981not computed but specified from file.  
     982\end{description} 
     983 
     984\np{nn\_isf}~=~1 and \np{nn\_isf}~=~2 compute a melt rate based on the water masse properties, ocean velocities and depth. 
     985 This flux is thus highly dependent of the model resolution (horizontal and vertical), realism of the water masse onto the shelf ... 
     986 
     987\np{nn\_isf}~=~3 and \np{nn\_isf}~=~4 read the melt rate and heat flux from a file. You have total control of the fwf scenario. 
     988 
     989 This can be usefull if the water masses on the shelf are not realistic or the resolution (horizontal/vertical) are too  
     990coarse to have realistic melting or for sensitivity studies where you want to control your input.  
     991Full description, sensitivity and validation in preparation.  
     992 
     993There is 2 ways to apply the fwf to NEMO. The first possibility (\np{ln\_divisf}~=~false) applied the fwf 
     994 and heat flux directly on the salinity and temperature tendancy. The second possibility (\np{ln\_divisf}~=~true) 
     995 apply the fwf as for the runoff fwf (see \S\ref{SBC_rnf}). The mass/volume addition due to the ice shelf melting is, 
     996 at each relevant depth level, added to the horizontal divergence (\textit{hdivn}) in the subroutine \rou{sbc\_isf\_div}  
     997(called from \mdl{divcur}).  
     998% 
    949999% ================================================================ 
    9501000%        Handling of icebergs 

  • trunk/DOC/TexFiles/Chapters/Chap_ZDF.tex

    r4147 r5120  
    830830% Bottom Friction 
    831831% ================================================================ 
    832 \section  [Bottom Friction (\textit{zdfbfr})]   {Bottom Friction (\mdl{zdfbfr} module)} 
     832\section  [Bottom and top Friction (\textit{zdfbfr})]   {Bottom Friction (\mdl{zdfbfr} module)} 
    833833\label{ZDF_bfr} 
    834834 
    … …  
    837837%-------------------------------------------------------------------------------------------------------------- 
    838838 
    839 Options are defined through the  \ngn{nambfr} namelist variables. 
     839Options to define the top and bottom friction are defined through the  \ngn{nambfr} namelist variables. 
     840The top friction is activated only if the ice shelf cavities are opened (\np{ln\_isfcav}~=~true). 
     841As the friction processes at the top and bottom are the represented similarly, only the bottom friction is described in detail. 
     842 
    840843Both the surface momentum flux (wind stress) and the bottom momentum  
    841844flux (bottom friction) enter the equations as a condition on the vertical  

  • trunk/DOC/TexFiles/Namelist/nambfr

    r4147 r5120  
    55                           !                              = 2 : nonlinear friction 
    66   rn_bfri1    =    4.e-4  !  bottom drag coefficient (linear case) 
    7    rn_bfri2    =    1.e-3  !  bottom drag coefficient (non linear case) 
     7   rn_bfri2    =    1.e-3  !  bottom drag coefficient (non linear case). Minimum coeft if ln_loglayer=T 
     8   rn_bfri2_max =   1.e-1  !  max. bottom drag coefficient (non linear case and ln_loglayer=T) 
    89   rn_bfeb2    =    2.5e-3 !  bottom turbulent kinetic energy background  (m2/s2) 
    9    rn_bfrz0    =    3.e-3  ! bottom roughness for loglayer bfr coeff  
     10   rn_bfrz0    =    3.e-3  !  bottom roughness [m] if ln_loglayer=T 
    1011   ln_bfr2d    = .false.   !  horizontal variation of the bottom friction coef (read a 2D mask file ) 
    1112   rn_bfrien   =    50.    !  local multiplying factor of bfr (ln_bfr2d=T) 
     13   rn_tfri1    =    4.e-4  !  top drag coefficient (linear case) 
     14   rn_tfri2    =    2.5e-3 !  top drag coefficient (non linear case). Minimum coeft if ln_loglayer=T 
     15   rn_tfri2_max =   1.e-1  !  max. top drag coefficient (non linear case and ln_loglayer=T) 
     16   rn_tfeb2    =    0.0    !  top turbulent kinetic energy background  (m2/s2) 
     17   rn_tfrz0    =    3.e-3  !  top roughness [m] if ln_loglayer=T 
     18   ln_tfr2d    = .false.   !  horizontal variation of the top friction coef (read a 2D mask file ) 
     19   rn_tfrien   =    50.    !  local multiplying factor of tfr (ln_tfr2d=T) 
     20 
    1221   ln_bfrimp   = .true.    !  implicit bottom friction (requires ln_zdfexp = .false. if true) 
     22   ln_loglayer = .false.   !  logarithmic formulation (non linear case) 
    1323/ 

  • trunk/DOC/TexFiles/Namelist/namdyn_hpg

    r4147 r5120  
    55   ln_hpg_zps  = .true.    !  z-coordinate - partial steps (interpolation) 
    66   ln_hpg_sco  = .false.   !  s-coordinate (standard jacobian formulation) 
     7   ln_hpg_isf  = .false.   !  s-coordinate (sco ) adapted to ice shelf cavity 
    78   ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial) 
    89   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme) 

  • trunk/DOC/TexFiles/Namelist/namsbc

    r4230 r5120  
    1919   ln_dm2dc    = .false.   !  daily mean to diurnal cycle on short wave 
    2020   ln_rnf      = .true.    !  runoffs                                   (T => fill namsbc_rnf) 
     21   nn_isf      = 0         !  ice shelf melting/freezing                (/=0 => fill namsbc_isf) 
     22                           !  0 =no isf                  1 = presence of ISF 
     23                           !  2 = bg03 parametrisation   3 = rnf file for isf 
     24                           !  4 = ISF fwf specified 
     25                           !  option 1 and 4 need ln_isfcav = .true. (domzgr) 
    2126   ln_ssr      = .true.    !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr) 
    2227   nn_fwb      = 3         !  FreshWater Budget: =0 unchecked 

  • trunk/DOC/TexFiles/Namelist/namzgr

    r3294 r5120  
    55   ln_zps      = .true.    !  z-coordinate - partial steps   (T/F) 
    66   ln_sco      = .false.   !  s- or hybrid z-s-coordinate    (T/F) 
     7   ln_isfcav   = .false.   !  ice shelf cavity               (T/F) 
    78/ 

  • trunk/NEMOGCM/CONFIG/ISOMIP/EXP00/namelist_cfg

    r4924 r5120  
    457457!----------------------------------------------------------------------- 
    458458   ln_hpg_zps  = .false.   !  z-coordinate - partial steps (interpolation) 
    459    ln_hpg_sco  = .true.    !  s-coordinate (standard jacobian formulation) 
     459   ln_hpg_isf  = .true.    !  s-coordinate adapted for isf (standard jacobian formulation) 
    460460   ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T) 
    461461                                 !           centered      time scheme  (F) 

  • trunk/NEMOGCM/CONFIG/SHARED/namelist_ref

    r5118 r5120  
    8585   ln_zps      = .true.    !  z-coordinate - partial steps   (T/F) 
    8686   ln_sco      = .false.   !  s- or hybrid z-s-coordinate    (T/F) 
    87    ln_isfcav   = .false.   !  ice shelf cavity 
     87   ln_isfcav   = .false.   !  ice shelf cavity               (T/F) 
    8888/ 
    8989!----------------------------------------------------------------------- 
    … …  
    846846   ln_hpg_zps  = .true.    !  z-coordinate - partial steps (interpolation) 
    847847   ln_hpg_sco  = .false.   !  s-coordinate (standard jacobian formulation) 
     848   ln_hpg_isf  = .false.   !  s-coordinate (sco ) adapted to isf 
    848849   ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial) 
    849850   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme) 

  • trunk/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90

    r4990 r5120  
    537537      CALL eos_rab( pts, rab_n )       ! now local thermal/haline expension ratio at T-points 
    538538      CALL bn2    ( pts, rab_n, rn2  ) ! now    Brunt-Vaisala 
    539       IF( ln_zps )   & ! Partial steps: before Horizontal DErivative 
    540         &    CALL zps_hde( kt, jpts, pts, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
    541         &                                        rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             ! 
    542         &                                 gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    )  ! of t, s, rd at the last ocean level 
     539        ! Partial steps: before Horizontal DErivative 
    543540        ! only gtsu, gtsv, rhd, gru , grv are used  
     541      IF( ln_zps  .AND. .NOT. ln_isfcav)                            & 
     542         &            CALL zps_hde    ( kt, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
     543         &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
     544      IF( ln_zps .AND.        ln_isfcav)                            & 
     545         &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
     546         &                                        rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     547         &                                 gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level 
     548 
    544549 
    545550 

  • trunk/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90

    r5118 r5120  
    233233         WRITE(numout,*) '                       NEMO team' 
    234234         WRITE(numout,*) '            Ocean General Circulation Model' 
    235          WRITE(numout,*) '                  version 3.5  (2012) ' 
     235         WRITE(numout,*) '                  version 3.6  (2015) ' 
    236236         WRITE(numout,*) 
    237237         WRITE(numout,*) 

  • trunk/NEMOGCM/NEMO/OOO_SRC/nemogcm.F90

    r5118 r5120  
    233233         WRITE(numout,*) '                       NEMO team' 
    234234         WRITE(numout,*) '            Ocean General Circulation Model' 
    235          WRITE(numout,*) '                  version 3.4  (2011) ' 
     235         WRITE(numout,*) '                  version 3.6  (2015) ' 
    236236         WRITE(numout,*) 
    237237         WRITE(numout,*) 

  • trunk/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r4998 r5120  
    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 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

    r5107 r5120  
    103103      END DO 
    104104      IF( .NOT.lk_vvl ) THEN 
    105          DO ji=1,jpi 
    106             DO jj=1,jpj 
    107                zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
    108             END DO 
    109          END DO 
     105         IF ( ln_isfcav ) THEN 
     106            DO ji=1,jpi 
     107               DO jj=1,jpj 
     108                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
     109               END DO 
     110            END DO 
     111         ELSE 
     112            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     113         END IF 
    110114      END IF 
    111115      !                                          
    … …  
    125129      END DO 
    126130      IF( .NOT.lk_vvl ) THEN 
    127          DO ji=1,jpi 
    128             DO jj=1,jpj 
    129                zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
    130             END DO 
    131          END DO 
     131         IF ( ln_isfcav ) THEN 
     132            DO ji=1,jpi 
     133               DO jj=1,jpj 
     134                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
     135               END DO 
     136            END DO 
     137         ELSE 
     138            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     139         END IF 
    132140      END IF 
    133141      !     
    … …  
    155163      END DO 
    156164      IF( .NOT.lk_vvl ) THEN 
    157          DO ji=1,jpi 
    158             DO jj=1,jpj 
    159                ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem)  
    160                zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal)  
    161             END DO 
    162          END DO 
     165         IF ( ln_isfcav ) THEN 
     166            DO ji=1,jpi 
     167               DO jj=1,jpj 
     168                  ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem)  
     169                  zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal)  
     170               END DO 
     171            END DO 
     172         ELSE 
     173            ztemp = ztemp + zarea_ssh(:,:) * tsn(:,:,1,jp_tem)  
     174            zsal  = zsal  + zarea_ssh(:,:) * tsn(:,:,1,jp_sal)  
     175         END IF 
    163176      ENDIF 
    164177      IF( lk_mpp ) THEN   

  • trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90

    r4990 r5120  
    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 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r4990 r5120  
    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 
    … …  
    400401         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk), STAT=ierr(11) ) 
    401402 
     403      ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) ) 
     404 
    402405#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) ) 
     406      ALLOCATE( npcoa(4,jpk), nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), STAT=ierr(12) ) 
    404407#endif 
    405408      ! 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r4990 r5120  
    281281      CALL lbc_lnk( fmask_i, 'F', 1._wp ) 
    282282 
     283      ! 3. Ocean/land mask at wu-, wv- and w points  
     284      !---------------------------------------------- 
     285      wmask (:,:,1) = tmask(:,:,1) ! ???????? 
     286      wumask(:,:,1) = umask(:,:,1) ! ???????? 
     287      wvmask(:,:,1) = vmask(:,:,1) ! ???????? 
     288      DO jk=2,jpk 
     289         wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1) 
     290         wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1)    
     291         wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1) 
     292      END DO 
    283293 
    284294      ! 4. ocean/land mask for the elliptic equation 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r5107 r5120  
    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 
    … …  
    589584      !! * Local declarations 
    590585      INTEGER                             :: ji,jj,jk       ! dummy loop indices 
     586      REAL(wp)                            :: zcoef 
    591587      !!---------------------------------------------------------------------- 
    592588 
    … …  
    635631      ! t- and w- points depth 
    636632      ! ---------------------- 
     633      ! set the isf depth as it is in the initial step 
    637634      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 
    638635      fsdepw_n(:,:,1) = 0.0_wp 
    639636      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 
    640       DO jj = 1,jpj 
    641          DO ji = 1,jpi 
    642             DO jk = 2,mikt(ji,jj)-1 
    643                fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 
    644                fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    645                fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 
    646             END DO 
    647             IF (mikt(ji,jj) .GT. 1) THEN 
    648                jk = mikt(ji,jj) 
    649                fsdept_n(ji,jj,jk) = gdepw_0(ji,jj,jk) + 0.5_wp * fse3w_n(ji,jj,jk) 
    650                fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
    651                fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
    652             END IF 
    653             DO jk = mikt(ji,jj)+1, jpk 
    654                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)) 
    655644               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 
    656                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) 
    657648            END DO 
    658649         END DO 
    659650      END DO 
     651 
    660652      ! Local depth and Inverse of the local depth of the water column at u- and v- points 
    661653      ! ---------------------------------------------------------------------------------- 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r5118 r5120  
    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! 
    … …  
    472474         ! 
    473475         ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 
    474          IF( cp_cfg == "isomip" ) THEN  
    475            !  
    476            risfdep(:,:)=200.e0  
    477            misfdep(:,:)=1  
    478            ij0 = 1 ; ij1 = 40  
    479            DO jj = mj0(ij0), mj1(ij1)  
    480               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp  
    481                 END DO  
     476         IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN  
     477            risfdep(:,:)=200.e0  
     478            misfdep(:,:)=1  
     479            ij0 = 1 ; ij1 = 40  
     480            DO jj = mj0(ij0), mj1(ij1)  
     481               risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp  
     482            END DO  
    482483            WHERE( bathy(:,:) <= 0._wp )  risfdep(:,:) = 0._wp  
    483            !  
    484          ELSEIF ( cp_cfg == "isomip2" ) THEN 
     484         !  
     485         ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN 
    485486         !  
    486487            risfdep(:,:)=0.e0 
    … …  
    540541            END IF 
    541542            CALL iom_close( inum ) 
    542             !   
     543            !                                                 
    543544            risfdep(:,:)=0._wp          
    544545            misfdep(:,:)=1              
    … …  
    588589      IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==! 
    589590         ! patch to avoid case bathy = ice shelf draft and bathy between 0 and zhmin 
    590          WHERE (bathy == risfdep) 
    591             bathy   = 0.0_wp ; risfdep = 0.0_wp 
    592          END WHERE 
     591         IF ( ln_isfcav ) THEN 
     592            WHERE (bathy == risfdep) 
     593               bathy   = 0.0_wp ; risfdep = 0.0_wp 
     594            END WHERE 
     595         END IF 
    593596         ! end patch 
    594597         IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level 
    … …  
    965968      !!---------------------------------------------------------------------- 
    966969      !! 
     970      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
     971      INTEGER  ::   ik, it           ! temporary integers 
     972      LOGICAL  ::   ll_print         ! Allow  control print for debugging 
     973      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
     974      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
     975      REAL(wp) ::   zmax             ! Maximum depth 
     976      REAL(wp) ::   zdiff            ! temporary scalar 
     977      REAL(wp) ::   zrefdep          ! temporary scalar 
     978      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt 
     979      !!--------------------------------------------------------------------- 
     980      ! 
     981      IF( nn_timing == 1 )  CALL timing_start('zgr_zps') 
     982      ! 
     983      CALL wrk_alloc( jpi, jpj, jpk, zprt ) 
     984      ! 
     985      IF(lwp) WRITE(numout,*) 
     986      IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps' 
     987      IF(lwp) WRITE(numout,*) '    ~~~~~~~ ' 
     988      IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used' 
     989 
     990      ll_print = .FALSE.                   ! Local variable for debugging 
     991       
     992      IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth 
     993         WRITE(numout,*) 
     994         WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)' 
     995         CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 
     996      ENDIF 
     997 
     998 
     999      ! bathymetry in level (from bathy_meter) 
     1000      ! =================== 
     1001      zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 
     1002      bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat) 
     1003      WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0 
     1004      ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level 
     1005      END WHERE 
     1006 
     1007      ! Compute mbathy for ocean points (i.e. the number of ocean levels) 
     1008      ! find the number of ocean levels such that the last level thickness 
     1009      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 
     1010      ! e3t_1d is the reference level thickness 
     1011      DO jk = jpkm1, 1, -1 
     1012         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
     1013         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
     1014      END DO 
     1015 
     1016      IF ( ln_isfcav ) CALL zgr_isf 
     1017 
     1018      ! Scale factors and depth at T- and W-points 
     1019      DO jk = 1, jpk                        ! intitialization to the reference z-coordinate 
     1020         gdept_0(:,:,jk) = gdept_1d(jk) 
     1021         gdepw_0(:,:,jk) = gdepw_1d(jk) 
     1022         e3t_0  (:,:,jk) = e3t_1d  (jk) 
     1023         e3w_0  (:,:,jk) = e3w_1d  (jk) 
     1024      END DO 
     1025      !  
     1026      DO jj = 1, jpj 
     1027         DO ji = 1, jpi 
     1028            ik = mbathy(ji,jj) 
     1029            IF( ik > 0 ) THEN               ! ocean point only 
     1030               ! max ocean level case 
     1031               IF( ik == jpkm1 ) THEN 
     1032                  zdepwp = bathy(ji,jj) 
     1033                  ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
     1034                  ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
     1035                  e3t_0(ji,jj,ik  ) = ze3tp 
     1036                  e3t_0(ji,jj,ik+1) = ze3tp 
     1037                  e3w_0(ji,jj,ik  ) = ze3wp 
     1038                  e3w_0(ji,jj,ik+1) = ze3tp 
     1039                  gdepw_0(ji,jj,ik+1) = zdepwp 
     1040                  gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
     1041                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
     1042                  ! 
     1043               ELSE                         ! standard case 
     1044                  IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
     1045                  ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
     1046                  ENDIF 
     1047!gm Bug?  check the gdepw_1d 
     1048                  !       ... on ik 
     1049                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
     1050                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
     1051                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
     1052                  e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
     1053                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
     1054                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
     1055                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
     1056                  !       ... on ik+1 
     1057                  e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1058                  e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     1059                  gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
     1060               ENDIF 
     1061            ENDIF 
     1062         END DO 
     1063      END DO 
     1064      ! 
     1065      it = 0 
     1066      DO jj = 1, jpj 
     1067         DO ji = 1, jpi 
     1068            ik = mbathy(ji,jj) 
     1069            IF( ik > 0 ) THEN               ! ocean point only 
     1070               e3tp (ji,jj) = e3t_0(ji,jj,ik) 
     1071               e3wp (ji,jj) = e3w_0(ji,jj,ik) 
     1072               ! test 
     1073               zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  ) 
     1074               IF( zdiff <= 0._wp .AND. lwp ) THEN  
     1075                  it = it + 1 
     1076                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
     1077                  WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
     1078                  WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
     1079                  WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  ) 
     1080               ENDIF 
     1081            ENDIF 
     1082         END DO 
     1083      END DO 
     1084      ! 
     1085      IF ( ln_isfcav ) THEN 
     1086      ! (ISF) Definition of e3t, u, v, w for ISF case 
     1087         DO jj = 1, jpj  
     1088            DO ji = 1, jpi  
     1089               ik = misfdep(ji,jj)  
     1090               IF( ik > 1 ) THEN               ! ice shelf point only  
     1091                  IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
     1092                  gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
     1093!gm Bug?  check the gdepw_0  
     1094               !       ... on ik  
     1095                  gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
     1096                     &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
     1097                     &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
     1098                  e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
     1099                  e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
     1100 
     1101                  IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
     1102                     e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
     1103                  ENDIF  
     1104               !       ... on ik / ik-1  
     1105                  e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
     1106                  e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
     1107! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
     1108                  gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
     1109               ENDIF  
     1110            END DO  
     1111         END DO  
     1112      !  
     1113         it = 0  
     1114         DO jj = 1, jpj  
     1115            DO ji = 1, jpi  
     1116               ik = misfdep(ji,jj)  
     1117               IF( ik > 1 ) THEN               ! ice shelf point only  
     1118                  e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
     1119                  e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
     1120               ! test  
     1121                  zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
     1122                  IF( zdiff <= 0. .AND. lwp ) THEN   
     1123                     it = it + 1  
     1124                     WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
     1125                     WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
     1126                     WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
     1127                     WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
     1128                  ENDIF  
     1129               ENDIF  
     1130            END DO  
     1131         END DO  
     1132      END IF 
     1133      ! END (ISF) 
     1134 
     1135      ! Scale factors and depth at U-, V-, UW and VW-points 
     1136      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     1137         e3u_0 (:,:,jk) = e3t_1d(jk) 
     1138         e3v_0 (:,:,jk) = e3t_1d(jk) 
     1139         e3uw_0(:,:,jk) = e3w_1d(jk) 
     1140         e3vw_0(:,:,jk) = e3w_1d(jk) 
     1141      END DO 
     1142      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
     1143         DO jj = 1, jpjm1 
     1144            DO ji = 1, fs_jpim1   ! vector opt. 
     1145               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
     1146               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
     1147               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
     1148               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
     1149            END DO 
     1150         END DO 
     1151      END DO 
     1152      IF ( ln_isfcav ) THEN 
     1153      ! (ISF) define e3uw (adapted for 2 cells in the water column) 
     1154      ! Need to test if the modification of only mikt and mbkt levels is enough 
     1155         DO jk = 2,jpk                           
     1156            DO jj = 1, jpjm1  
     1157               DO ji = 1, fs_jpim1   ! vector opt.  
     1158                  e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj  ,jk) ) & 
     1159                    &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj  ,jk-1) ) 
     1160                  e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji  ,jj+1,jk) ) & 
     1161                    &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji  ,jj+1,jk-1) ) 
     1162               END DO  
     1163            END DO  
     1164         END DO 
     1165      END IF 
     1166       
     1167      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
     1168      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
     1169      ! 
     1170      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1171         WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
     1172         WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk) 
     1173         WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk) 
     1174         WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk) 
     1175      END DO 
     1176       
     1177      ! Scale factor at F-point 
     1178      DO jk = 1, jpk                        ! initialisation to z-scale factors 
     1179         e3f_0(:,:,jk) = e3t_1d(jk) 
     1180      END DO 
     1181      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
     1182         DO jj = 1, jpjm1 
     1183            DO ji = 1, fs_jpim1   ! vector opt. 
     1184               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
     1185            END DO 
     1186         END DO 
     1187      END DO 
     1188      CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions 
     1189      ! 
     1190      DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
     1191         WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk) 
     1192      END DO 
     1193!!gm  bug ? :  must be a do loop with mj0,mj1 
     1194      !  
     1195      e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
     1196      e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)  
     1197      e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)  
     1198      e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)  
     1199      e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)  
     1200 
     1201      ! Control of the sign 
     1202      IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
     1203      IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
     1204      IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
     1205      IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
     1206      
     1207      ! Compute gdep3w_0 (vertical sum of e3w) 
     1208      IF ( ln_isfcav ) THEN ! if cavity 
     1209         WHERE (misfdep == 0) misfdep = 1 
     1210         DO jj = 1,jpj 
     1211            DO ji = 1,jpi 
     1212               gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
     1213               DO jk = 2, misfdep(ji,jj) 
     1214                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1215               END DO 
     1216               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)) 
     1217               DO jk = misfdep(ji,jj) + 1, jpk 
     1218                  gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1219               END DO 
     1220            END DO 
     1221         END DO 
     1222      ELSE ! no cavity 
     1223         gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 
     1224         DO jk = 2, jpk 
     1225            gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
     1226         END DO 
     1227      END IF 
     1228      !                                               ! ================= ! 
     1229      IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
     1230         !                                            ! ================= ! 
     1231         DO jj = 1,jpj 
     1232            DO ji = 1, jpi 
     1233               ik = MAX( mbathy(ji,jj), 1 ) 
     1234               zprt(ji,jj,1) = e3t_0   (ji,jj,ik) 
     1235               zprt(ji,jj,2) = e3w_0   (ji,jj,ik) 
     1236               zprt(ji,jj,3) = e3u_0   (ji,jj,ik) 
     1237               zprt(ji,jj,4) = e3v_0   (ji,jj,ik) 
     1238               zprt(ji,jj,5) = e3f_0   (ji,jj,ik) 
     1239               zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 
     1240            END DO 
     1241         END DO 
     1242         WRITE(numout,*) 
     1243         WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1244         WRITE(numout,*) 
     1245         WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1246         WRITE(numout,*) 
     1247         WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1248         WRITE(numout,*) 
     1249         WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1250         WRITE(numout,*) 
     1251         WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1252         WRITE(numout,*) 
     1253         WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
     1254      ENDIF   
     1255      ! 
     1256      CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 
     1257      ! 
     1258      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps') 
     1259      ! 
     1260   END SUBROUTINE zgr_zps 
     1261 
     1262   SUBROUTINE zgr_isf 
     1263      !!---------------------------------------------------------------------- 
     1264      !!                    ***  ROUTINE zgr_isf  *** 
     1265      !!    
     1266      !! ** Purpose :   check the bathymetry in levels 
     1267      !!    
     1268      !! ** Method  :   THe water column have to contained at least 2 cells 
     1269      !!                Bathymetry and isfdraft are modified (dig/close) to respect 
     1270      !!                this criterion. 
     1271      !!                  
     1272      !!    
     1273      !! ** Action  : - test compatibility between isfdraft and bathy  
     1274      !!              - bathy and isfdraft are modified 
     1275      !!---------------------------------------------------------------------- 
     1276      !!    
    9671277      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    9681278      INTEGER  ::   ik, it           ! temporary integers 
    … …  
    9751285      REAL(wp) ::   zdiff            ! temporary scalar 
    9761286      REAL(wp) ::   zrefdep          ! temporary scalar 
    977       REAL(wp) ::   zbathydiff, zrisfdepdiff  
    978       REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 3D workspace (ISH) 
    979       INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep   ! 3D workspace (ISH) 
    980       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt 
     1287      REAL(wp) ::   zbathydiff, zrisfdepdiff  ! isf temporary scalar 
     1288      REAL(wp), POINTER, DIMENSION(:,:)   ::   zrisfdep, zbathy, zmask   ! 2D workspace (ISH) 
     1289      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmisfdep         ! 2D workspace (ISH) 
    9811290      !!--------------------------------------------------------------------- 
    9821291      ! 
    983       IF( nn_timing == 1 )  CALL timing_start('zgr_zps') 
    984       ! 
    985       CALL wrk_alloc( jpi, jpj, jpk, zprt ) 
     1292      IF( nn_timing == 1 )  CALL timing_start('zgr_isf') 
     1293      ! 
    9861294      CALL wrk_alloc( jpi, jpj, zbathy, zmask, zrisfdep) 
    987       CALL wrk_alloc( jpi, jpj, zmbathy, zmisfdep) 
    988       ! 
    989       IF(lwp) WRITE(numout,*) 
    990       IF(lwp) WRITE(numout,*) '    zgr_zps : z-coordinate with partial steps' 
    991       IF(lwp) WRITE(numout,*) '    ~~~~~~~ ' 
    992       IF(lwp) WRITE(numout,*) '              mbathy is recomputed : bathy_level file is NOT used' 
    993  
    994       ll_print = .FALSE.                   ! Local variable for debugging 
    995        
    996       IF(lwp .AND. ll_print) THEN          ! control print of the ocean depth 
    997          WRITE(numout,*) 
    998          WRITE(numout,*) 'dom_zgr_zps:  bathy (in hundred of meters)' 
    999          CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 
    1000       ENDIF 
    1001  
    1002       ! bathymetry in level (from bathy_meter) 
    1003       ! =================== 
    1004       zmax = gdepw_1d(jpk) + e3t_1d(jpk)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 
    1005       bathy(:,:) = MIN( zmax ,  bathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat) 
    1006       WHERE( bathy(:,:) == 0._wp )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0 
    1007       ELSE WHERE                     ;   mbathy(:,:) = jpkm1   ! ocean : initialize mbathy to the max ocean level 
    1008       END WHERE 
    1009  
    1010       ! Compute mbathy for ocean points (i.e. the number of ocean levels) 
    1011       ! find the number of ocean levels such that the last level thickness 
    1012       ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 
    1013       ! e3t_1d is the reference level thickness 
    1014       DO jk = jpkm1, 1, -1 
    1015          zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 
    1016          WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
    1017       END DO 
     1295      CALL wrk_alloc( jpi, jpj, zmisfdep, zmbathy ) 
     1296 
     1297 
    10181298      ! (ISF) compute misfdep 
    10191299      WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) .NE. 0) ;   misfdep(:,:) = 1   ! open water : set misfdep to 1   
    … …  
    10591339            misfdep(jpi,:) = misfdep(  2  ,:)  
    10601340         ENDIF 
    1061   
     1341 
    10621342         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
    10631343            mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west 
    10641344            mbathy(jpi,:) = mbathy(  2  ,:) 
    10651345         ENDIF 
    1066   
     1346 
    10671347         ! split last cell if possible (only where water column is 2 cell or less) 
    10681348         DO jk = jpkm1, 1, -1 
    … …  
    10821362            END WHERE 
    10831363         END DO 
    1084   
     1364 
    10851365  
    10861366 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 
    … …  
    13631643               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0 
    13641644               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    1365                IF( ibtest == 0 ) THEN 
     1645               IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN 
    13661646                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 
    13671647               END IF 
    … …  
    14791759      ENDIF  
    14801760 
    1481       ! Scale factors and depth at T- and W-points 
    1482       DO jk = 1, jpk                        ! intitialization to the reference z-coordinate 
    1483          gdept_0(:,:,jk) = gdept_1d(jk) 
    1484          gdepw_0(:,:,jk) = gdepw_1d(jk) 
    1485          e3t_0  (:,:,jk) = e3t_1d  (jk) 
    1486          e3w_0  (:,:,jk) = e3w_1d  (jk) 
    1487       END DO 
    1488       !  
    1489       DO jj = 1, jpj 
    1490          DO ji = 1, jpi 
    1491             ik = mbathy(ji,jj) 
    1492             IF( ik > 0 ) THEN               ! ocean point only 
    1493                ! max ocean level case 
    1494                IF( ik == jpkm1 ) THEN 
    1495                   zdepwp = bathy(ji,jj) 
    1496                   ze3tp  = bathy(ji,jj) - gdepw_1d(ik) 
    1497                   ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
    1498                   e3t_0(ji,jj,ik  ) = ze3tp 
    1499                   e3t_0(ji,jj,ik+1) = ze3tp 
    1500                   e3w_0(ji,jj,ik  ) = ze3wp 
    1501                   e3w_0(ji,jj,ik+1) = ze3tp 
    1502                   gdepw_0(ji,jj,ik+1) = zdepwp 
    1503                   gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
    1504                   gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
    1505                   ! 
    1506                ELSE                         ! standard case 
    1507                   IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 
    1508                   ELSE                                       ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
    1509                   ENDIF 
    1510 !gm Bug?  check the gdepw_1d 
    1511                   !       ... on ik 
    1512                   gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
    1513                      &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
    1514                      &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
    1515                   e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
    1516                      &                          / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
    1517                   e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
    1518                      &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
    1519                   !       ... on ik+1 
    1520                   e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1521                   e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
    1522                   gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
    1523                ENDIF 
    1524             ENDIF 
    1525          END DO 
    1526       END DO 
    1527       ! 
    1528       it = 0 
    1529       DO jj = 1, jpj 
    1530          DO ji = 1, jpi 
    1531             ik = mbathy(ji,jj) 
    1532             IF( ik > 0 ) THEN               ! ocean point only 
    1533                e3tp (ji,jj) = e3t_0(ji,jj,ik) 
    1534                e3wp (ji,jj) = e3w_0(ji,jj,ik) 
    1535                ! test 
    1536                zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik  ) 
    1537                IF( zdiff <= 0._wp .AND. lwp ) THEN  
    1538                   it = it + 1 
    1539                   WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj 
    1540                   WRITE(numout,*) ' bathy = ', bathy(ji,jj) 
    1541                   WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 
    1542                   WRITE(numout,*) ' e3tp    = ', e3t_0  (ji,jj,ik), ' e3wp    = ', e3w_0  (ji,jj,ik  ) 
    1543                ENDIF 
    1544             ENDIF 
    1545          END DO 
    1546       END DO 
    1547       ! 
    1548       ! (ISF) Definition of e3t, u, v, w for ISF case 
    1549       DO jj = 1, jpj  
    1550          DO ji = 1, jpi  
    1551             ik = misfdep(ji,jj)  
    1552             IF( ik > 1 ) THEN               ! ice shelf point only  
    1553                IF( risfdep(ji,jj) < gdepw_1d(ik) )  risfdep(ji,jj)= gdepw_1d(ik)  
    1554                gdepw_0(ji,jj,ik) = risfdep(ji,jj)  
    1555 !gm Bug?  check the gdepw_0  
    1556                !       ... on ik  
    1557                gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
    1558                   &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
    1559                   &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
    1560                e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
    1561                e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
    1562  
    1563                IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only (2 cell water column)  
    1564                   e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
    1565                 ENDIF  
    1566                !       ... on ik / ik-1  
    1567                e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
    1568                e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
    1569 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
    1570                gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
    1571             ENDIF  
    1572          END DO  
    1573       END DO  
    1574       !  
    1575       it = 0  
    1576       DO jj = 1, jpj  
    1577          DO ji = 1, jpi  
    1578             ik = misfdep(ji,jj)  
    1579             IF( ik > 1 ) THEN               ! ice shelf point only  
    1580                e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
    1581                e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
    1582                ! test  
    1583                zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
    1584                IF( zdiff <= 0. .AND. lwp ) THEN   
    1585                   it = it + 1  
    1586                   WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
    1587                   WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)  
    1588                   WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
    1589                   WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
    1590                ENDIF  
    1591             ENDIF  
    1592          END DO  
    1593       END DO  
    1594       ! END (ISF) 
    1595  
    1596       ! Scale factors and depth at U-, V-, UW and VW-points 
    1597       DO jk = 1, jpk                        ! initialisation to z-scale factors 
    1598          e3u_0 (:,:,jk) = e3t_1d(jk) 
    1599          e3v_0 (:,:,jk) = e3t_1d(jk) 
    1600          e3uw_0(:,:,jk) = e3w_1d(jk) 
    1601          e3vw_0(:,:,jk) = e3w_1d(jk) 
    1602       END DO 
    1603       DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors 
    1604          DO jj = 1, jpjm1 
    1605             DO ji = 1, fs_jpim1   ! vector opt. 
    1606                e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 
    1607                e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 
    1608                e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 
    1609                e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 
    1610             END DO 
    1611          END DO 
    1612       END DO 
    1613       ! (ISF) define e3uw 
    1614       DO jk = 2,jpk                           
    1615          DO jj = 1, jpjm1  
    1616             DO ji = 1, fs_jpim1   ! vector opt.  
    1617                e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj  ,jk) ) & 
    1618                  &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj  ,jk-1) ) 
    1619                e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji  ,jj+1,jk) ) & 
    1620                  &   - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji  ,jj+1,jk-1) ) 
    1621             END DO  
    1622          END DO  
    1623       END DO 
    1624       !End (ISF) 
    1625        
    1626       CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
    1627       CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    1628       ! 
    1629       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1630          WHERE( e3u_0 (:,:,jk) == 0._wp )   e3u_0 (:,:,jk) = e3t_1d(jk) 
    1631          WHERE( e3v_0 (:,:,jk) == 0._wp )   e3v_0 (:,:,jk) = e3t_1d(jk) 
    1632          WHERE( e3uw_0(:,:,jk) == 0._wp )   e3uw_0(:,:,jk) = e3w_1d(jk) 
    1633          WHERE( e3vw_0(:,:,jk) == 0._wp )   e3vw_0(:,:,jk) = e3w_1d(jk) 
    1634       END DO 
    1635        
    1636       ! Scale factor at F-point 
    1637       DO jk = 1, jpk                        ! initialisation to z-scale factors 
    1638          e3f_0(:,:,jk) = e3t_1d(jk) 
    1639       END DO 
    1640       DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors 
    1641          DO jj = 1, jpjm1 
    1642             DO ji = 1, fs_jpim1   ! vector opt. 
    1643                e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 
    1644             END DO 
    1645          END DO 
    1646       END DO 
    1647       CALL lbc_lnk( e3f_0, 'F', 1._wp )       ! Lateral boundary conditions 
    1648       ! 
    1649       DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries) 
    1650          WHERE( e3f_0(:,:,jk) == 0._wp )   e3f_0(:,:,jk) = e3t_1d(jk) 
    1651       END DO 
    1652 !!gm  bug ? :  must be a do loop with mj0,mj1 
    1653       !  
    1654       e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:)     ! we duplicate factor scales for jj = 1 and jj = 2 
    1655       e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)  
    1656       e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)  
    1657       e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)  
    1658       e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)  
    1659  
    1660       ! Control of the sign 
    1661       IF( MINVAL( e3t_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3t_0 <= 0' ) 
    1662       IF( MINVAL( e3w_0  (:,:,:) ) <= 0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   e3w_0 <= 0' ) 
    1663       IF( MINVAL( gdept_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdept_0 <  0' ) 
    1664       IF( MINVAL( gdepw_0(:,:,:) ) <  0._wp )   CALL ctl_stop( '    zgr_zps :   e r r o r   gdepw_0 <  0' ) 
    1665       
    1666       ! Compute gdep3w_0 (vertical sum of e3w) 
    1667       WHERE (misfdep == 0) misfdep = 1 
    1668       DO jj = 1,jpj 
    1669          DO ji = 1,jpi 
    1670             gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
    1671             DO jk = 2, misfdep(ji,jj) 
    1672                gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    1673             END DO 
    1674             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)) 
    1675             DO jk = misfdep(ji,jj) + 1, jpk 
    1676                gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
    1677             END DO 
    1678         END DO 
    1679       END DO 
    1680       !                                               ! ================= ! 
    1681       IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
    1682          !                                            ! ================= ! 
    1683          DO jj = 1,jpj 
    1684             DO ji = 1, jpi 
    1685                ik = MAX( mbathy(ji,jj), 1 ) 
    1686                zprt(ji,jj,1) = e3t_0   (ji,jj,ik) 
    1687                zprt(ji,jj,2) = e3w_0   (ji,jj,ik) 
    1688                zprt(ji,jj,3) = e3u_0   (ji,jj,ik) 
    1689                zprt(ji,jj,4) = e3v_0   (ji,jj,ik) 
    1690                zprt(ji,jj,5) = e3f_0   (ji,jj,ik) 
    1691                zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 
    1692             END DO 
    1693          END DO 
    1694          WRITE(numout,*) 
    1695          WRITE(numout,*) 'domzgr e3t(mbathy)'      ;   CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1696          WRITE(numout,*) 
    1697          WRITE(numout,*) 'domzgr e3w(mbathy)'      ;   CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1698          WRITE(numout,*) 
    1699          WRITE(numout,*) 'domzgr e3u(mbathy)'      ;   CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1700          WRITE(numout,*) 
    1701          WRITE(numout,*) 'domzgr e3v(mbathy)'      ;   CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1702          WRITE(numout,*) 
    1703          WRITE(numout,*) 'domzgr e3f(mbathy)'      ;   CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1704          WRITE(numout,*) 
    1705          WRITE(numout,*) 'domzgr gdep3w(mbathy)'   ;   CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 
    1706       ENDIF   
    1707       ! 
    1708       CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 
    17091761      CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 
    17101762      CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 
    1711       ! 
    1712       IF( nn_timing == 1 )  CALL timing_stop('zgr_zps') 
    1713       ! 
    1714    END SUBROUTINE zgr_zps 
     1763 
     1764      IF( nn_timing == 1 )  CALL timing_stop('zgr_isf') 
     1765       
     1766   END SUBROUTINE 
    17151767 
    17161768   SUBROUTINE zgr_sco 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r4990 r5120  
    137137         CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )        ! before potential and in situ densities 
    138138#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 
     139         IF( ln_zps .AND. .NOT. ln_isfcav)                                 & 
     140            &            CALL zps_hde    ( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
     141            &                                            rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
     142         IF( ln_zps .AND.       ln_isfcav)                                 & 
     143            &            CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
     144            &                                            rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     145            &                                     gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
    142146#endif 
    143147         !    

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90

    r4990 r5120  
    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 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90

    r4990 r5120  
    127127      IF( ln_dynzad_zts .AND. .NOT. ln_dynadv_vec )   & 
    128128          CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' ) 
     129      IF( ln_dynzad_zts .AND. ln_isfcav )   & 
     130          CALL ctl_stop( 'Sub timestepping of vertical advection does not work with ln_isfcav = .TRUE.' ) 
    129131 
    130132      !                               ! Set nadv 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90

    r4990 r5120  
    8080              ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) 
    8181              va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) 
    82                
    83               ! (ISF) stability criteria for top friction 
    84               ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels 
    85               ikbv = mikv(ji,jj) 
    86               ! 
    87               ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
    88               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  tfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) & 
    89                  &             * (1.-umask(ji,jj,1)) 
    90               va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  tfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) & 
    91                  &             * (1.-vmask(ji,jj,1)) 
    92               ! (ISF) 
    93                
    9482           END DO 
    9583        END DO 
     84         
     85        IF ( ln_isfcav ) THEN 
     86           DO jj = 2, jpjm1 
     87              DO ji = 2, jpim1 
     88                 ! (ISF) stability criteria for top friction 
     89                 ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels 
     90                 ikbv = mikv(ji,jj) 
     91                 ! 
     92                 ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
     93                 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  tfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) & 
     94                    &             * (1.-umask(ji,jj,1)) 
     95                 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  tfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) & 
     96                    &             * (1.-vmask(ji,jj,1)) 
     97                 ! (ISF) 
     98              END DO 
     99           END DO 
     100        END IF 
    96101 
    97102        ! 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r4990 r5120  
    1616   !!            3.4  !  2011-11  (H. Liu) hpg_prj: Original code for s-coordinates 
    1717   !!                 !           (A. Coward) suppression of hel, wdj and rot options 
     18   !!            3.6  !  2014-11  (P. Mathiot) hpg_isf: original code for ice shelf cavity 
    1819   !!---------------------------------------------------------------------- 
    1920 
    … …  
    2526   !!       hpg_zps  : z-coordinate plus partial steps (interpolation) 
    2627   !!       hpg_sco  : s-coordinate (standard jacobian formulation) 
     28   !!       hpg_isf  : s-coordinate (sco formulation) adapted to ice shelf 
    2729   !!       hpg_djc  : s-coordinate (Density Jacobian with Cubic polynomial) 
    2830   !!       hpg_prj  : s-coordinate (Pressure Jacobian with Cubic polynomial) 
    … …  
    5557   LOGICAL , PUBLIC ::   ln_hpg_djc      !: s-coordinate (Density Jacobian with Cubic polynomial) 
    5658   LOGICAL , PUBLIC ::   ln_hpg_prj      !: s-coordinate (Pressure Jacobian scheme) 
     59   LOGICAL , PUBLIC ::   ln_hpg_isf      !: s-coordinate similar to sco modify for isf 
    5760   LOGICAL , PUBLIC ::   ln_dynhpg_imp   !: semi-implicite hpg flag 
    5861 
    … …  
    97100      CASE (  3 )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial) 
    98101      CASE (  4 )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme) 
     102      CASE (  5 )   ;   CALL hpg_isf    ( kt )      ! s-coordinate similar to sco modify for ice shelf 
    99103      END SELECT 
    100104      ! 
    … …  
    128132      !! 
    129133      NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     & 
    130          &                 ln_hpg_djc, ln_hpg_prj, ln_dynhpg_imp 
     134         &                 ln_hpg_djc, ln_hpg_prj, ln_hpg_isf, ln_dynhpg_imp 
    131135      !!---------------------------------------------------------------------- 
    132136      ! 
    … …  
    148152         WRITE(numout,*) '      z-coord. - partial steps (interpolation)          ln_hpg_zps    = ', ln_hpg_zps 
    149153         WRITE(numout,*) '      s-coord. (standard jacobian formulation)          ln_hpg_sco    = ', ln_hpg_sco 
     154         WRITE(numout,*) '      s-coord. (standard jacobian formulation) for isf  ln_hpg_isf    = ', ln_hpg_isf 
    150155         WRITE(numout,*) '      s-coord. (Density Jacobian: Cubic polynomial)     ln_hpg_djc    = ', ln_hpg_djc 
    151156         WRITE(numout,*) '      s-coord. (Pressure Jacobian: Cubic polynomial)    ln_hpg_prj    = ', ln_hpg_prj 
    … …  
    158163                           & either  ln_hpg_sco or  ln_hpg_prj instead') 
    159164      ! 
    160       IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj) )   & 
     165      IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) )   & 
    161166         &   CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:& 
    162167                           & the standard jacobian formulation hpg_sco or & 
    163168                           & the pressure jacobian formulation hpg_prj') 
     169 
     170      IF(       ln_hpg_isf .AND. .NOT. ln_isfcav )   & 
     171         &   CALL ctl_stop( ' hpg_isf not available if ln_isfcav = false ' ) 
     172      IF( .NOT. ln_hpg_isf .AND.       ln_isfcav )   & 
     173         &   CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' ) 
    164174      ! 
    165175      !                               ! Set nhpg from ln_hpg_... flags 
    … …  
    169179      IF( ln_hpg_djc )   nhpg = 3 
    170180      IF( ln_hpg_prj )   nhpg = 4 
     181      IF( ln_hpg_isf )   nhpg = 5 
    171182      ! 
    172183      !                               ! Consistency check 
    … …  
    177188      IF( ln_hpg_djc )   ioptio = ioptio + 1 
    178189      IF( ln_hpg_prj )   ioptio = ioptio + 1 
     190      IF( ln_hpg_isf )   ioptio = ioptio + 1 
    179191      IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 
    180       IF( (ln_hpg_zco .OR. ln_hpg_zps .OR. ln_hpg_djc .OR. ln_hpg_prj ) .AND. nn_isf .NE. 0 )   & 
    181           &  CALL ctl_stop( 'Only hpg_sco has been corrected to work with ice shelf cavity.' ) 
     192      !  
     193      ! initialisation of ice load 
     194      riceload(:,:)=0.0 
    182195      ! 
    183196   END SUBROUTINE dyn_hpg_init 
    … …  
    345358   END SUBROUTINE hpg_zps 
    346359 
    347  
    348360   SUBROUTINE hpg_sco( kt ) 
    349361      !!--------------------------------------------------------------------- 
    … …  
    366378      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    367379      !! 
     380      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
     381      REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars 
     382      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
     383      !!---------------------------------------------------------------------- 
     384      ! 
     385      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 
     386      ! 
     387      IF( kt == nit000 ) THEN 
     388         IF(lwp) WRITE(numout,*) 
     389         IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' 
     390         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used' 
     391      ENDIF 
     392 
     393      ! Local constant initialization 
     394      zcoef0 = - grav * 0.5_wp 
     395      ! To use density and not density anomaly 
     396      IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume 
     397      ELSE                 ;     znad = 0._wp         ! Fixed volume 
     398      ENDIF 
     399 
     400      ! Surface value 
     401      DO jj = 2, jpjm1 
     402         DO ji = fs_2, fs_jpim1   ! vector opt. 
     403            ! hydrostatic pressure gradient along s-surfaces 
     404            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   & 
     405               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) ) 
     406            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   & 
     407               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) ) 
     408            ! s-coordinate pressure gradient correction 
     409            zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
     410               &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 
     411            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
     412               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
     413            ! add to the general momentum trend 
     414            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
     415            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 
     416         END DO 
     417      END DO 
     418 
     419      ! interior value (2=<jk=<jpkm1) 
     420      DO jk = 2, jpkm1 
     421         DO jj = 2, jpjm1 
     422            DO ji = fs_2, fs_jpim1   ! vector opt. 
     423               ! hydrostatic pressure gradient along s-surfaces 
     424               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   & 
     425                  &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   & 
     426                  &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  ) 
     427               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   & 
     428                  &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   & 
     429                  &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  ) 
     430               ! s-coordinate pressure gradient correction 
     431               zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
     432                  &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 
     433               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
     434                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
     435               ! add to the general momentum trend 
     436               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
     437               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 
     438            END DO 
     439         END DO 
     440      END DO 
     441      ! 
     442      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
     443      ! 
     444   END SUBROUTINE hpg_sco 
     445 
     446   SUBROUTINE hpg_isf( kt ) 
     447      !!--------------------------------------------------------------------- 
     448      !!                  ***  ROUTINE hpg_sco  *** 
     449      !! 
     450      !! ** Method  :   s-coordinate case. Jacobian scheme. 
     451      !!      The now hydrostatic pressure gradient at a given level, jk, 
     452      !!      is computed by taking the vertical integral of the in-situ 
     453      !!      density gradient along the model level from the suface to that 
     454      !!      level. s-coordinates (ln_sco): a corrective term is added 
     455      !!      to the horizontal pressure gradient : 
     456      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ] 
     457      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ] 
     458      !!      add it to the general momentum trend (ua,va). 
     459      !!         ua = ua - 1/e1u * zhpi 
     460      !!         va = va - 1/e2v * zhpj 
     461      !!      iceload is added and partial cell case are added to the top and bottom 
     462      !!       
     463      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
     464      !!---------------------------------------------------------------------- 
     465      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     466      !! 
    368467      INTEGER  ::   ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j                 ! dummy loop indices 
    369468      REAL(wp) ::   zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1             ! temporary scalars 
    … …  
    379478     IF( kt == nit000 ) THEN 
    380479         IF(lwp) WRITE(numout,*) 
    381          IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' 
     480         IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 
    382481         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used' 
    383482      ENDIF 
    … …  
    565664!================================================================================== 
    566665 
    567 # if defined key_vectopt_loop 
    568          jj = 1 
    569          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    570 # else 
    571666      DO jj = 2, jpjm1 
    572667         DO ji = 2, jpim1 
    573 # endif 
    574668            iku = mbku(ji,jj) 
    575669            ikv = mbkv(ji,jj) 
    … …  
    598692               va(ji,jj,ikv)   =  va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 
    599693            END IF 
    600 # if ! defined key_vectopt_loop 
    601          END DO 
    602 # endif 
     694         END DO 
    603695      END DO 
    604696      
    … …  
    610702      CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 
    611703      ! 
    612    END SUBROUTINE hpg_sco 
     704   END SUBROUTINE hpg_isf 
    613705 
    614706 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

    r4990 r5120  
    250250      IF( ( ioptio > 1 .AND. .NOT. lk_esopa ) .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) )   & 
    251251           &   CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' ) 
    252       IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. nn_isf .NE. 0 )   & 
     252      IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. ln_isfcav )   & 
    253253           &   CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' ) 
    254254      ! 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r5032 r5120  
    2222   USE dom_oce         ! ocean space and time domain 
    2323   USE sbc_oce         ! surface boundary condition: ocean 
     24   USE sbcisf          ! ice shelf variable (fwfisf) 
    2425   USE dynspg_oce      ! surface pressure gradient variables 
    2526   USE phycst          ! physical constants 
    … …  
    453454      !                                         ! Surface net water flux and rivers 
    454455      IF (ln_bt_fw) THEN 
    455          zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) ) 
     456         zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) ) 
    456457      ELSE 
    457          zssh_frc(:,:) = zraur * z1_2 * (emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)) 
     458         zssh_frc(:,:) = zraur * z1_2 * (  emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)   & 
     459                &                        + rdivisf * ( fwfisf(:,:) + fwfisf_b(:,:) )       ) 
    458460      ENDIF 
    459461#if defined key_asminc 

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

    r4990 r5120  
    9595         END DO    
    9696      END DO 
    97       DO jj = 2, jpjm1              ! Surface and bottom values set to zero 
    98          DO ji = fs_2, fs_jpim1           ! vector opt. 
    99             zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp 
    100             zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp 
    101             zwuw(ji,jj,jpk) = 0._wp 
    102             zwvw(ji,jj,jpk) = 0._wp 
    103          END DO   
    104       END DO 
     97      ! 
     98      ! Surface and bottom advective fluxes set to zero 
     99      IF ( ln_isfcav ) THEN 
     100         DO jj = 2, jpjm1 
     101            DO ji = fs_2, fs_jpim1           ! vector opt. 
     102               zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp 
     103               zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp 
     104               zwuw(ji,jj,jpk) = 0._wp 
     105               zwvw(ji,jj,jpk) = 0._wp 
     106            END DO 
     107         END DO 
     108      ELSE 
     109         DO jj = 2, jpjm1         
     110            DO ji = fs_2, fs_jpim1           ! vector opt. 
     111               zwuw(ji,jj, 1 ) = 0._wp 
     112               zwvw(ji,jj, 1 ) = 0._wp 
     113               zwuw(ji,jj,jpk) = 0._wp 
     114               zwvw(ji,jj,jpk) = 0._wp 
     115            END DO   
     116         END DO 
     117      END IF 
    105118 
    106119      DO jk = 1, jpkm1              ! Vertical momentum advection at u- and v-points 
    … …  
    196209         END DO 
    197210      END DO 
    198  
    199       DO jj = 2, jpjm1                    ! Surface and bottom advective fluxes set to zero 
     211      ! 
     212      ! Surface and bottom advective fluxes set to zero 
     213      DO jj = 2, jpjm1         
    200214         DO ji = fs_2, fs_jpim1           ! vector opt. 
    201             zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp 
    202             zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp 
     215            zwuw(ji,jj, 1 ) = 0._wp 
     216            zwvw(ji,jj, 1 ) = 0._wp 
    203217            zwuw(ji,jj,jpk) = 0._wp 
    204218            zwvw(ji,jj,jpk) = 0._wp 
    … …  
    228242            DO jj = 2, jpjm1                 ! vertical momentum advection at w-point 
    229243               DO ji = fs_2, fs_jpim1        ! vector opt. 
    230                   zwuw(ji,jj,jk) = ( zww(ji+1,jj  ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) 
    231                   zwvw(ji,jj,jk) = ( zww(ji  ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) 
     244                  zwuw(ji,jj,jk) = ( zww(ji+1,jj  ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) !* wumask(ji,jj,jk) 
     245                  zwvw(ji,jj,jk) = ( zww(ji  ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) !* wvmask(ji,jj,jk) 
    232246               END DO   
    233247            END DO    

  • trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r4990 r5120  
    105105               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 
    106106               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 
    107                ikbu = miku(ji,jj)       ! ocean top level at u- and v-points  
    108                ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
    109                IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu) 
    110                IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv) 
    111             END DO 
    112          END DO 
     107            END DO 
     108         END DO 
     109         IF ( ln_isfcav ) THEN 
     110            DO jj = 2, jpjm1 
     111               DO ji = 2, jpim1 
     112                  ikbu = miku(ji,jj)       ! ocean top level at u- and v-points  
     113                  ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
     114                  IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu) 
     115                  IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv) 
     116               END DO 
     117            END DO 
     118         END IF 
    113119      ENDIF 
    114120 
    … …  
    145151               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua 
    146152               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va 
    147                ikbu = miku(ji,jj)         ! top ocean level at u- and v-points  
    148                ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
    149                ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu) 
    150                ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv) 
    151                ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua 
    152                va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va 
    153             END DO 
    154          END DO 
     153            END DO 
     154         END DO 
     155         IF ( ln_isfcav ) THEN 
     156            DO jj = 2, jpjm1         
     157               DO ji = fs_2, fs_jpim1   ! vector opt. 
     158                  ikbu = miku(ji,jj)         ! top ocean level at u- and v-points  
     159                  ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
     160                  ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu) 
     161                  ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv) 
     162                  ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua 
     163                  va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va 
     164               END DO 
     165            END DO 
     166         END IF 
    155167      ENDIF 
    156168#endif 
    … …  
    167179               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl   * fse3u_a(ji,jj,jk)   ! after scale factor at T-point 
    168180               zcoef = - p2dt / ze3ua       
    169                zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
    170                zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk) 
    171                zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)  
    172                zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1) 
    173                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 
     181               zzwi          = zcoef * avmu  (ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
     182               zwi(ji,jj,jk) = zzwi  * wumask(ji,jj,jk  ) 
     183               zzws          = zcoef * avmu  (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)  
     184               zws(ji,jj,jk) = zzws  * wumask(ji,jj,jk+1) 
     185               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    174186            END DO 
    175187         END DO 
    … …  
    198210      ! 
    199211      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    200       DO jj = 2, jpjm1    
    201          DO ji = fs_2, fs_jpim1   ! vector opt. 
    202             DO jk = miku(ji,jj)+1, jpkm1 
     212      DO jk = 2, jpkm1 
     213         DO jj = 2, jpjm1    
     214            DO ji = fs_2, fs_jpim1   ! vector opt. 
    203215               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    204216            END DO 
    … …  
    208220      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    209221         DO ji = fs_2, fs_jpim1   ! vector opt. 
    210             ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,miku(ji,jj)) + r_vvl   * fse3u_a(ji,jj,miku(ji,jj))  
    211222#if defined key_dynspg_ts 
    212             ua(ji,jj,miku(ji,jj)) = ua(ji,jj,miku(ji,jj)) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    213                &                                      / ( ze3ua * rau0 )  
     223            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1)  
     224            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     225               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1)  
    214226#else 
    215             ua(ji,jj,miku(ji,jj)) = ub(ji,jj,miku(ji,jj)) & 
    216                &                   + p2dt *(ua(ji,jj,miku(ji,jj)) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    217                &                                  / ( fse3u(ji,jj,miku(ji,jj)) * rau0     ) )  
    218 #endif 
    219             DO jk = miku(ji,jj)+1, jpkm1 
     227            ua(ji,jj,1) = ub(ji,jj,1) & 
     228               &                   + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     229               &                                      / ( fse3u(ji,jj,1) * rau0     ) * umask(ji,jj,1) )  
     230#endif 
     231         END DO 
     232      END DO 
     233      DO jk = 2, jpkm1 
     234         DO jj = 2, jpjm1 
     235            DO ji = fs_2, fs_jpim1 
    220236#if defined key_dynspg_ts 
    221237               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side 
    … …  
    231247         DO ji = fs_2, fs_jpim1   ! vector opt. 
    232248            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    233             DO jk = jpk-2, miku(ji,jj), -1 
     249         END DO 
     250      END DO 
     251      DO jk = jpk-2, 1, -1 
     252         DO jj = 2, jpjm1 
     253            DO ji = fs_2, fs_jpim1 
    234254               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    235255            END DO 
    … …  
    260280               zcoef = - p2dt / ze3va 
    261281               zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
    262                zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk) 
     282               zwi(ji,jj,jk) =  zzwi * wvmask(ji,jj,jk) 
    263283               zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 
    264                zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1) 
    265                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 
     284               zws(ji,jj,jk) =  zzws * wvmask(ji,jj,jk+1) 
     285               zwd(ji,jj,jk) = 1._wp - zzwi - zzws 
    266286            END DO 
    267287         END DO 
    … …  
    290310      ! 
    291311      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    292       DO jj = 2, jpjm1    
    293          DO ji = fs_2, fs_jpim1   ! vector opt. 
    294             DO jk = mikv(ji,jj)+1, jpkm1         
     312      DO jk = 2, jpkm1         
     313         DO jj = 2, jpjm1    
     314            DO ji = fs_2, fs_jpim1   ! vector opt. 
    295315               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    296316            END DO 
    … …  
    300320      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    301321         DO ji = fs_2, fs_jpim1   ! vector opt. 
    302             ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,mikv(ji,jj)) + r_vvl   * fse3v_a(ji,jj,mikv(ji,jj))  
    303322#if defined key_dynspg_ts             
    304             va(ji,jj,mikv(ji,jj)) = va(ji,jj,mikv(ji,jj)) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     323            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1)  
     324            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    305325               &                                      / ( ze3va * rau0 )  
    306326#else 
    307             va(ji,jj,mikv(ji,jj)) = vb(ji,jj,mikv(ji,jj)) & 
    308                &                   + p2dt *(va(ji,jj,mikv(ji,jj)) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    309                &                                                       / ( fse3v(ji,jj,mikv(ji,jj)) * rau0     )  ) 
    310 #endif 
    311             DO jk = mikv(ji,jj)+1, jpkm1 
     327            va(ji,jj,1) = vb(ji,jj,1) & 
     328               &                   + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     329               &                                                       / ( fse3v(ji,jj,1) * rau0     )  ) 
     330#endif 
     331         END DO 
     332      END DO 
     333      DO jk = 2, jpkm1 
     334         DO jj = 2, jpjm1 
     335            DO ji = fs_2, fs_jpim1   ! vector opt. 
    312336#if defined key_dynspg_ts 
    313337               zrhs = va(ji,jj,jk)   ! zrhs=right hand side 
    … …  
    323347         DO ji = fs_2, fs_jpim1   ! vector opt. 
    324348            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    325             DO jk = jpk-2, mikv(ji,jj), -1 
     349         END DO 
     350      END DO 
     351      DO jk = jpk-2, 1, -1 
     352         DO jj = 2, jpjm1 
     353            DO ji = fs_2, fs_jpim1 
    326354               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    327355            END DO 
    … …  
    349377              avmu(ji,jj,ikbu+1) = 0.e0 
    350378              avmv(ji,jj,ikbv+1) = 0.e0 
    351               ikbu = miku(ji,jj)         ! ocean top level at u- and v-points  
    352               ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
    353               IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0 
    354               IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0 
    355379           END DO 
    356380        END DO 
     381        IF (ln_isfcav) THEN 
     382           DO jj = 2, jpjm1 
     383              DO ji = 2, jpim1 
     384                 ikbu = miku(ji,jj)         ! ocean top level at u- and v-points  
     385                 ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
     386                 IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0 
     387                 IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0 
     388              END DO 
     389           END DO 
     390        END IF 
    357391      ENDIF 
    358392      ! 

  • trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

    r5016 r5120  
    142142            DO jj = 1, jpjm1 
    143143               DO ji = 1, jpim1 
    144 ! IF should be useless check zpshde (PM) 
    145                IF ( mbku(ji,jj) > 1 ) zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 
    146                IF ( mbkv(ji,jj) > 1 ) zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 
     144                  zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 
     145                  zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 
     146               END DO 
     147            END DO 
     148         ENDIF 
     149         IF( ln_zps .AND. ln_isfcav ) THEN           ! partial steps correction at the bottom ocean level 
     150            DO jj = 1, jpjm1 
     151               DO ji = 1, jpim1 
    147152               IF ( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj)  
    148153               IF ( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj) 
    … …  
    151156         ENDIF 
    152157         ! 
    153          zdzr(:,:,1) = 0._wp        !==   Local vertical density gradient at T-point   == !   (evaluated from N^2) 
    154          DO jk = 1, jpkm1 
     158         !==   Local vertical density gradient at T-point   == !   (evaluated from N^2) 
     159         ! interior value 
     160         DO jk = 2, jpkm1 
    155161            !                                ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 
    156162            !                                !   trick: tmask(ik  )  = 0   =>   all pn2   = 0   =>   zdzr = 0 
    … …  
    162168         END DO 
    163169         ! surface initialisation  
    164          DO jj = 1, jpjm1 
    165             DO ji = 1, jpim1 
    166               zdzr(ji,jj,1:mikt(ji,jj)) = 0._wp 
    167             END DO 
    168          END DO 
     170         zdzr(:,:,1) = 0._wp  
     171         IF ( ln_isfcav ) THEN 
     172            ! if isf need to overwrite the interior value at at the first ocean point 
     173            DO jj = 1, jpjm1 
     174               DO ji = 1, jpim1 
     175                  zdzr(ji,jj,1:mikt(ji,jj)) = 0._wp 
     176               END DO 
     177            END DO 
     178         END IF 
    169179         ! 
    170180         !                          !==   Slopes just below the mixed layer   ==! 
    … …  
    175185         ! ===========================      | vslp = d/dj( prd ) / d/dz( prd ) 
    176186         ! 
    177          DO jj = 2, jpjm1 
    178             DO ji = fs_2, fs_jpim1   ! vector opt. 
    179                IF (miku(ji,jj) .GT. miku(ji+1,jj)) zhmlpu(ji,jj) = hmlpt(ji  ,jj) 
    180                IF (miku(ji,jj) .LT. miku(ji+1,jj)) zhmlpu(ji,jj) = hmlpt(ji+1,jj) 
    181                IF (miku(ji,jj) .EQ. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj), hmlpt(ji+1,jj)) 
    182                IF (mikv(ji,jj) .GT. miku(ji,jj+1)) zhmlpv(ji,jj) = hmlpt(ji  ,jj) 
    183                IF (mikv(ji,jj) .LT. miku(ji,jj+1)) zhmlpv(ji,jj) = hmlpt(ji,jj+1) 
    184                IF (mikv(ji,jj) .EQ. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji,jj+1)) 
     187         IF ( ln_isfcav ) THEN 
     188            DO jj = 2, jpjm1 
     189               DO ji = fs_2, fs_jpim1   ! vector opt. 
     190                  IF (miku(ji,jj) .GT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp) 
     191                  IF (miku(ji,jj) .LT. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji+1,jj  ),                   5._wp) 
     192                  IF (miku(ji,jj) .EQ. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji+1,jj  ), 5._wp) 
     193                  IF (mikv(ji,jj) .GT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ),                   5._wp) 
     194                  IF (mikv(ji,jj) .LT. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj+1),                   5._wp) 
     195                  IF (mikv(ji,jj) .EQ. miku(ji,jj+1)) zhmlpv(ji,jj) = MAX(hmlpt(ji  ,jj  ), hmlpt(ji  ,jj+1), 5._wp) 
     196               ENDDO 
    185197            ENDDO 
    186          ENDDO 
     198         ELSE 
     199            DO jj = 2, jpjm1 
     200               DO ji = fs_2, fs_jpim1   ! vector opt. 
     201                  zhmlpu(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp) 
     202                  zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp) 
     203               ENDDO 
     204            ENDDO 
     205         END IF 
    187206         DO jk = 2, jpkm1                            !* Slopes at u and v points 
    188207            DO jj = 2, jpjm1 
    … …  
    198217                  zbv = MIN(  zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav )  ) 
    199218                  !                                      ! uslp and vslp output in zwz and zww, resp. 
    200                   zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )  
    201                   zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )  
     219                  zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj  ,jk) ) 
     220                  zfj = MAX( omlmask(ji,jj,jk), omlmask(ji  ,jj+1,jk) ) 
    202221                  ! thickness of water column between surface and level k at u/v point 
    203                   zdepu = 0.5_wp * (( fsdept(ji,jj,jk) + fsdept(ji+1,jj  ,jk) )                   & 
    204                              - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj  ) )  & 
    205                              - fse3u(ji,jj,miku(ji,jj))                                         ) 
    206                   zdepv = 0.5_wp * (( fsdept(ji,jj,jk) + fsdept(ji  ,jj+1,jk) )                   & 
    207                              - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) & 
    208                              - fse3v(ji,jj,mikv(ji,jj))                                         ) 
    209                   zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps )                                              & 
    210                      &                 + zfi  * uslpml(ji,jj)                                                     & 
    211                      &                        * zdepu / MAX( zhmlpu(ji,jj), 5._wp ) 
    212                   zwz(ji,jj,jk) = zwz(ji,jj,jk) * umask(ji,jj,jk) * umask(ji,jj,jk-1) 
    213                   zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps )                                              & 
    214                      &                 + zfj  * vslpml(ji,jj)                                                     & 
    215                      &                        * zdepv / MAX( zhmlpv(ji,jj), 5._wp )  
    216                   zww(ji,jj,jk) = zww(ji,jj,jk) * vmask(ji,jj,jk) * vmask(ji,jj,jk-1) 
     222                  zdepu = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji+1,jj  ,jk) )                              & 
     223                                   - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj  ) ) - fse3u(ji,jj,miku(ji,jj)) ) 
     224                  zdepv = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji  ,jj+1,jk) )                              & 
     225                                   - 2 * MAX( risfdep(ji,jj), risfdep(ji  ,jj+1) ) - fse3v(ji,jj,mikv(ji,jj)) ) 
     226                  ! 
     227                  zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps )                                          & 
     228                     &                 + zfi  * uslpml(ji,jj) * zdepu / zhmlpu(ji,jj) 
     229                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * wumask(ji,jj,jk) 
     230                  zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps )                                          & 
     231                     &                 + zfj  * vslpml(ji,jj) * zdepv / zhmlpv(ji,jj)  
     232                  zww(ji,jj,jk) = zww(ji,jj,jk) * wvmask(ji,jj,jk) 
    217233                   
    218234                  
    … …  
    266282                  uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk  ) ) * 0.5_wp   & 
    267283                     &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5_wp   & 
    268                      &                            *   umask(ji,jj,jk-1) !* umask(ji,jj,jk) * umask(ji,jj,jk+1) 
     284                     &                            *   umask(ji,jj,jk-1) 
    269285                  vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk  ) ) * 0.5_wp   & 
    270286                     &                            * ( vmask(ji  ,jj,jk) + vmask(ji  ,jj,jk+1) ) * 0.5_wp   & 
    271                      &                            *   vmask(ji,jj,jk-1) !* vmask(ji,jj,jk) * vmask(ji,jj,jk+1) 
     287                     &                            *   vmask(ji,jj,jk-1) 
    272288               END DO 
    273289            END DO 
    … …  
    282298               DO ji = fs_2, fs_jpim1   ! vector opt. 
    283299                  !                                  !* Local vertical density gradient evaluated from N^2 
    284                   zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
     300                  zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) * wmask(ji,jj,jk) 
    285301                  !                                  !* Slopes at w point 
    286302                  !                                        ! i- & j-gradient of density at w-points 
    … …  
    298314                  zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj )  ) 
    299315                  !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 
    300                   zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )    ! zfk=1 in the ML otherwise zfk=0 
     316                  zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0 
    301317                  zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj), 10._wp ) 
    302318                  zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) & 
    303                      &            + zck * wslpiml(ji,jj) * zfk  ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
     319                     &            + zck * wslpiml(ji,jj) * zfk  ) * wmask(ji,jj,jk) 
    304320                  zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) & 
    305                      &            + zck * wslpjml(ji,jj) * zfk  ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
     321                     &            + zck * wslpjml(ji,jj) * zfk  ) * wmask(ji,jj,jk) 
    306322 
    307323!!gm  modif to suppress omlmask....  (as in Griffies operator) 
    … …  
    356372                  zck =   ( umask(ji,jj,jk) + umask(ji-1,jj,jk) )   & 
    357373                     &  * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 
    358                   wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck * tmask(ji,jj,jk-1) * tmask(ji,jj,jk) 
    359                   wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck * tmask(ji,jj,jk-1) * tmask(ji,jj,jk) 
     374                  wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck * wmask(ji,jj,jk) 
     375                  wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck * wmask(ji,jj,jk) 
    360376               END DO 
    361377            END DO 
    … …  
    423439                  vslp(ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)  
    424440                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) & 
    425                     &                              * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * 0.5  
     441                    &                              * wmask(ji,jj,jk) * 0.5  
    426442                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) & 
    427                     &                              * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * 0.5  
     443                    &                              * wmask(ji,jj,jk) * 0.5  
    428444               END DO  
    429445            END DO  
    … …  
    736752            DO ji = 1, jpi 
    737753               ik = nmln(ji,jj) - 1 
    738                IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN   ;   omlmask(ji,jj,jk) = 1._wp 
    739                ELSE                  ;   omlmask(ji,jj,jk) = 0._wp 
     754               IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN 
     755                  omlmask(ji,jj,jk) = 1._wp 
     756               ELSE 
     757                  omlmask(ji,jj,jk) = 0._wp 
    740758               ENDIF 
    741759            END DO 
    … …  
    794812            zbj = MIN(  zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj )  ) 
    795813            !                        !- i- & j-slope at w-points (wslpiml, wslpjml) 
    796             wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 
    797             wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 
     814            wslpiml(ji,jj) = zai / ( zbi - zeps ) * wmask (ji,jj,ik) 
     815            wslpjml(ji,jj) = zaj / ( zbj - zeps ) * wmask (ji,jj,ik) 
    798816         END DO 
    799817      END DO 

  • trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r4990 r5120  
    9898   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   emp_tot           !: total E-P over ocean and ice                 [Kg/m2/s] 
    9999   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fmmflx            !: freshwater budget: freezing/melting          [Kg/m2/s] 
    100    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rnf    , rnf_b    !: river runoff   [Kg/m2/s]   
     100   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rnf    , rnf_b    !: river runoff        [Kg/m2/s]   
     101   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fwfisf , fwfisf_b !: ice shelf melting   [Kg/m2/s]   
    101102   !! 
    102103   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  sbc_tsc, sbc_tsc_b  !: sbc content trend                      [K.m/s] jpi,jpj,jpts 
    … …  
    147148         &      sfx    (jpi,jpj) , sfx_b(jpi,jpj) , emp_tot(jpi,jpj), fmmflx(jpi,jpj), STAT=ierr(2) ) 
    148149         ! 
    149       ALLOCATE( rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,     & 
    150          &      rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) ) 
     150      ALLOCATE( fwfisf  (jpi,jpj), rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,     & 
     151         &      fwfisf_b(jpi,jpj), rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) ) 
    151152         ! 
    152153      ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     & 

  • trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90

    r4990 r5120  
    88   !!            3.0  ! 2006-08  (G. Madec)  Surface module 
    99   !!            3.2  ! 2009-07  (C. Talandier) emp mean s spread over erp area  
     10   !!            3.6  ! 2014-11  (P. Mathiot  ) add ice shelf melting 
    1011   !!---------------------------------------------------------------------- 
    1112 
    … …  
    8889         ! 
    8990         IF( kn_fwb == 3 .AND. nn_sssr /= 2 )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 requires nn_sssr = 2, we stop ' ) 
    90          ! 
    91          area = glob_sum( e1e2t(:,:) )           ! interior global domain surface 
     91         IF( kn_fwb == 3 .AND. ln_isfcav    )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 with ln_isfcav = .TRUE. not working, we stop ' ) 
     92         ! 
     93         area = glob_sum( e1e2t(:,:) * tmask(:,:,1))           ! interior global domain surface 
     94         ! isf cavities are excluded because it can feedback to the melting with generation of inhibition of plumes 
     95         ! and in case of no melt, it can generate HSSW. 
    9296         ! 
    9397#if ! defined key_lim2 &&  ! defined key_lim3 && ! defined key_cice 
    … …  
    106110            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain 
    107111            zcoef = z_fwf * rcp 
    108             emp(:,:) = emp(:,:) - z_fwf  
    109             qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) ! account for change to the heat budget due to fw correction 
     112            emp(:,:) = emp(:,:) - z_fwf              * tmask(:,:,1) 
     113            qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction 
    110114         ENDIF 
    111115         ! 
    … …  
    138142         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! correct the freshwater fluxes 
    139143            zcoef = fwfold * rcp 
    140             emp(:,:) = emp(:,:) + fwfold 
    141             qns(:,:) = qns(:,:) - zcoef * sst_m(:,:) ! account for change to the heat budget due to fw correction 
     144            emp(:,:) = emp(:,:) + fwfold             * tmask(:,:,1) 
     145            qns(:,:) = qns(:,:) - zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction 
    142146         ENDIF 
    143147         ! 
    … …  
    158162            zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) ) 
    159163            !                                                  ! fwf global mean (excluding ocean to ice/snow exchanges)  
    160             z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - snwice_fmass(:,:) ) ) / area 
     164            z_fwf     = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) - snwice_fmass(:,:) ) ) / area 
    161165            !             
    162166            IF( z_fwf < 0._wp ) THEN         ! spread out over >0 erp area to increase evaporation 

  • trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r4990 r5120  
    77   !! History :  3.2   !  2011-02  (C.Harris  ) Original code isf cav 
    88   !!            X.X   !  2006-02  (C. Wang   ) Original code bg03 
    9    !!            3.4   !  2013-03  (P. Mathiot) Merging 
     9   !!            3.4   !  2013-03  (P. Mathiot) Merging + parametrization 
    1010   !!---------------------------------------------------------------------- 
    1111 
    … …  
    3737 
    3838   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc    
    39    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fwfisf_b, fwfisf  !: evaporation damping   [kg/m2/s] 
    40    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf            !: net heat flux from ice shelf 
     39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf              !: net heat flux from ice shelf 
    4140   REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m] 
    4241   LOGICAL , PUBLIC ::   ln_divisf                   !: flag to correct divergence  
    … …  
    309308      sbc_isf_alloc = 0       ! set to zero if no array to be allocated 
    310309      IF( .NOT. ALLOCATED( qisf ) ) THEN 
    311          ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts)              , & 
    312                &    qisf(jpi,jpj)     , fwfisf(jpi,jpj)     , fwfisf_b(jpi,jpj)   , & 
    313                &    rhisf_tbl(jpi,jpj), r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj)  , & 
    314                &    ttbl(jpi,jpj)     , stbl(jpi,jpj)       , utbl(jpi,jpj)       , & 
    315                &    vtbl(jpi, jpj)    , risfLeff(jpi,jpj)   , rhisf_tbl_0(jpi,jpj), & 
    316                &    ralpha(jpi,jpj)   , misfkt(jpi,jpj)     , misfkb(jpi,jpj)     , & 
     310         ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj)   , & 
     311               &    rhisf_tbl(jpi,jpj)    , r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj)  , & 
     312               &    ttbl(jpi,jpj)         , stbl(jpi,jpj)       , utbl(jpi,jpj)       , & 
     313               &    vtbl(jpi, jpj)        , risfLeff(jpi,jpj)   , rhisf_tbl_0(jpi,jpj), & 
     314               &    ralpha(jpi,jpj)       , misfkt(jpi,jpj)     , misfkb(jpi,jpj)     , & 
    317315               &    STAT= sbc_isf_alloc ) 
    318316         ! 

  • trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r4990 r5120  
    1313   !!            3.4  ! 2011-11  (C. Harris) CICE added as an option 
    1414   !!            3.5  ! 2012-11  (A. Coward, G. Madec) Rethink of heat, mass and salt surface fluxes 
     15   !!            3.6  ! 2014-11  (P. Mathiot, C. Harris) add ice shelves melting                     
    1516   !!---------------------------------------------------------------------- 
    1617 
    … …  
    179180         IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' ) 
    180181         fwfisf  (:,:) = 0.0_wp 
     182         fwfisf_b(:,:) = 0.0_wp 
    181183      END IF 
    182184      IF( nn_ice == 0  )   fr_i(:,:) = 0.e0        ! no ice in the domain, ice fraction is always zero 

  • trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90

    r4990 r5120  
    6161      !!--------------------------------------------------------------------- 
    6262       
    63       !                                        !* first wet T-, U-, V- ocean level (ISF) variables (T, S, depth, velocity) 
     63      !                                        !* surface T-, U-, V- ocean level variables (T, S, depth, velocity) 
    6464      DO jj = 1, jpj 
    6565         DO ji = 1, jpi 
    66             zub(ji,jj)        = ub (ji,jj,miku(ji,jj)) 
    67             zvb(ji,jj)        = vb (ji,jj,mikv(ji,jj)) 
    6866            zts(ji,jj,jp_tem) = tsn(ji,jj,mikt(ji,jj),jp_tem) 
    6967            zts(ji,jj,jp_sal) = tsn(ji,jj,mikt(ji,jj),jp_sal) 
    7068         END DO 
    7169      END DO 
     70      zub(:,:)        = ub (:,:,1       ) 
     71      zvb(:,:)        = vb (:,:,1       ) 
    7272      ! 
    7373      IF( lk_vvl ) THEN 
    74          DO jj = 1, jpj 
    75             DO ji = 1, jpi 
    76                zdep(ji,jj) = fse3t_n(ji,jj,mikt(ji,jj)) 
    77             END DO 
    78          END DO 
     74         zdep(:,:) = fse3t_n(:,:,1) 
    7975      ENDIF 
    8076      !                                                   ! ---------------------------------------- ! 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r4990 r5120  
    206206      IF( lk_esopa         )   ioptio =          1 
    207207 
    208       IF( ( ln_traadv_muscl .OR. ln_traadv_muscl2 .OR. ln_traadv_ubs .OR. ln_traadv_qck ) .AND. nn_isf .NE. 0 )  & 
     208      IF( ( ln_traadv_muscl .OR. ln_traadv_muscl2 .OR. ln_traadv_ubs .OR. ln_traadv_qck .OR. ln_traadv_tvd_zts ) .AND. ln_isfcav )  & 
    209209      &   CALL ctl_stop( 'Only traadv_cen2 and traadv_tvd is compatible with ice shelf cavity') 
    210210 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90

    r4990 r5120  
    106106      ENDIF 
    107107      ! 
    108       zwi(:,:,:) = 0.e0 ; zwz(:,:,:) = 0.e0 
     108      zwi(:,:,:) = 0.e0 ;  
    109109      ! 
    110110      !                                                          ! =========== 
    111111      DO jn = 1, kjpt                                            ! tracer loop 
    112112         !                                                       ! =========== 
    113          ! 1. Bottom value : flux set to zero 
     113         ! 1. Bottom and k=1 value : flux set to zero 
    114114         ! ---------------------------------- 
    115115         zwx(:,:,jpk) = 0.e0    ;    zwz(:,:,jpk) = 0.e0 
    116116         zwy(:,:,jpk) = 0.e0    ;    zwi(:,:,jpk) = 0.e0 
    117  
     117           
     118         zwz(:,:,1  ) = 0._wp 
    118119         ! 2. upstream advection with initial mass fluxes & intermediate update 
    119120         ! -------------------------------------------------------------------- 
    … …  
    134135 
    135136         ! upstream tracer flux in the k direction 
     137         ! Interior value 
     138         DO jk = 2, jpkm1 
     139            DO jj = 1, jpj 
     140               DO ji = 1, jpi 
     141                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
     142                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
     143                  zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 
     144               END DO 
     145            END DO 
     146         END DO 
    136147         ! Surface value 
    137148         IF( lk_vvl ) THEN    
    138             DO jj = 1, jpj 
    139                DO ji = 1, jpi 
    140                   zwz(ji,jj, mikt(ji,jj) ) = 0.e0                         ! volume variable 
    141                END DO 
    142             END DO 
     149            IF ( ln_isfcav ) THEN 
     150               DO jj = 1, jpj 
     151                  DO ji = 1, jpi 
     152                     zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable 
     153                  END DO 
     154               END DO 
     155            ELSE 
     156               zwz(:,:,1) = 0.e0          ! volume variable 
     157            END IF 
    143158         ELSE                 
    144             DO jj = 1, jpj 
    145                DO ji = 1, jpi 
    146                   zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
    147                END DO 
    148             END DO    
     159            IF ( ln_isfcav ) THEN 
     160               DO jj = 1, jpj 
     161                  DO ji = 1, jpi 
     162                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     163                  END DO 
     164               END DO    
     165            ELSE 
     166               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface 
     167            END IF 
    149168         ENDIF 
    150          ! Interior value 
    151          DO jj = 1, jpj 
    152             DO ji = 1, jpi 
    153                DO jk = mikt(ji,jj)+1, jpkm1 
    154                   zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
    155                   zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
    156                   zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) 
    157                END DO 
    158             END DO 
    159          END DO 
    160169 
    161170         ! total advective trend 
    … …  
    202211       
    203212         ! antidiffusive flux on k 
    204          zwz(:,:,1) = 0.e0         ! Surface value 
    205          ! 
    206          DO jj = 1, jpj 
    207             DO ji = 1, jpi 
    208                ik=mikt(ji,jj) 
    209                ! surface value 
    210                zwz(ji,jj,1:ik) = 0.e0 
    211                ! Interior value 
    212                DO jk = mikt(ji,jj)+1, jpkm1                     
     213         ! Interior value 
     214         DO jk = 2, jpkm1                     
     215            DO jj = 1, jpj 
     216               DO ji = 1, jpi 
    213217                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk) 
    214218               END DO 
    215219            END DO 
    216220         END DO 
     221         ! surface value 
     222         IF ( ln_isfcav ) THEN 
     223            DO jj = 1, jpj 
     224               DO ji = 1, jpi 
     225                  zwz(ji,jj,mikt(ji,jj)) = 0.e0 
     226               END DO 
     227            END DO 
     228         ELSE 
     229            zwz(:,:,1) = 0.e0 
     230         END IF 
    217231         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions 
    218232         CALL lbc_lnk( zwz, 'W',  1. ) 
    … …  
    358372 
    359373         ! upstream tracer flux in the k direction 
    360          ! Surface value 
    361          IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0._wp                        ! volume variable 
    362          ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface  
    363          ENDIF 
    364374         ! Interior value 
    365375         DO jk = 2, jpkm1 
    … …  
    372382            END DO 
    373383         END DO 
     384         ! Surface value 
     385         IF( lk_vvl ) THEN 
     386            IF ( ln_isfcav ) THEN 
     387               DO jj = 1, jpj 
     388                  DO ji = 1, jpi 
     389                     zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable +    isf 
     390                  END DO 
     391               END DO 
     392            ELSE 
     393               zwz(:,:,1) = 0.e0                              ! volume variable + no isf 
     394            END IF 
     395         ELSE 
     396            IF ( ln_isfcav ) THEN 
     397               DO jj = 1, jpj 
     398                  DO ji = 1, jpi 
     399                     zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface +    isf 
     400                  END DO 
     401               END DO 
     402            ELSE 
     403               zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)                                               ! linear free surface + no isf 
     404            END IF 
     405         ENDIF 
    374406 
    375407         ! total advective trend 
    … …  
    580612         &        paft * tmask + zbig * ( 1._wp - tmask )  ) 
    581613 
    582       DO jj = 2, jpjm1 
    583          DO ji = fs_2, fs_jpim1   ! vector opt. 
    584             DO jk = mikt(ji,jj), jpkm1 
    585                ikm1 = MAX(jk-1,mikt(ji,jj)) 
    586                z2dtt = p2dt(jk) 
    587                 
     614      DO jk = 1, jpkm1 
     615         ikm1 = MAX(jk-1,1) 
     616         z2dtt = p2dt(jk) 
     617         DO jj = 2, jpjm1 
     618            DO ji = fs_2, fs_jpim1   ! vector opt. 
     619 
    588620               ! search maximum in neighbourhood 
    589621               zup = MAX(  zbup(ji  ,jj  ,jk  ),   & 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r4990 r5120  
    290290      IF(lwp) WRITE(numout,*) '              homogeneous ocean T = ', zt0, ' S = ',zs0 
    291291 
     292      ! Initialisation of gtui/gtvi in case of no cavity 
     293      IF ( .NOT. ln_isfcav ) THEN 
     294         gtui(:,:,:) = 0.0_wp 
     295         gtvi(:,:,:) = 0.0_wp 
     296      END IF 
    292297      !                                        ! T & S profile (to be coded +namelist parameter 
    293298 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90

    r4990 r5120  
    116116               END DO 
    117117            END DO 
    118  
    119118            !                          !==  Laplacian  ==! 
    120119            ! 
    … …  
    125124               END DO 
    126125            END DO 
     126            ! 
    127127            IF( ln_zps ) THEN                ! set gradient at partial step level (last ocean level) 
    128128               DO jj = 1, jpjm1 
    … …  
    130130                     IF( mbku(ji,jj) == jk )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgu(ji,jj,jn) 
    131131                     IF( mbkv(ji,jj) == jk )  ztv(ji,jj,jk) = zeev(ji,jj) * pgv(ji,jj,jn) 
    132                      ! (ISH) 
    133                      IF( miku(ji,jj) == jk )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgui(ji,jj,jn) 
    134                      IF( mikv(ji,jj) == jk )  ztu(ji,jj,jk) = zeev(ji,jj) * pgvi(ji,jj,jn) 
    135132                  END DO 
    136133               END DO 
    137134            ENDIF 
     135            ! (ISH) 
     136            IF( ln_zps .AND. ln_isfcav ) THEN ! set gradient at partial step level (first ocean level in a cavity) 
     137               DO jj = 1, jpjm1 
     138                  DO ji = 1, jpim1 
     139                     IF( miku(ji,jj) == MAX(jk,2) )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgui(ji,jj,jn) 
     140                     IF( mikv(ji,jj) == MAX(jk,2) )  ztu(ji,jj,jk) = zeev(ji,jj) * pgvi(ji,jj,jn) 
     141                  END DO 
     142               END DO 
     143            ENDIF 
     144            ! 
    138145            DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient 
    139146               DO ji = fs_2, fs_jpim1   ! vector opt. 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90

    r4990 r5120  
    106106      ! 
    107107      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
     108      INTEGER  ::  ikt 
    108109      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars 
    109110      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      - 
    … …  
    149150            END DO 
    150151         END DO 
     152 
     153         ! partial cell correction 
    151154         IF( ln_zps ) THEN      ! partial steps correction at the last ocean level  
    152155            DO jj = 1, jpjm1 
    153156               DO ji = 1, fs_jpim1   ! vector opt. 
    154157! IF useless if zpshde defines pgu everywhere 
    155                   IF (mbku(ji,jj) > 1) zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
    156                   IF (mbkv(ji,jj) > 1) zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
    157                   ! (ISF) 
     158                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
     159                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
     160               END DO 
     161            END DO 
     162         ENDIF 
     163         IF( ln_zps .AND. ln_isfcav ) THEN      ! partial steps correction at the first wet level beneath a cavity 
     164            DO jj = 1, jpjm1 
     165               DO ji = 1, fs_jpim1   ! vector opt. 
    158166                  IF (miku(ji,jj) > 1) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)           
    159167                  IF (mikv(ji,jj) > 1) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)      
    160168               END DO 
    161169            END DO 
    162          ENDIF 
     170         END IF 
    163171 
    164172         !!---------------------------------------------------------------------- 
    165173         !!   II - horizontal trend  (full) 
    166174         !!---------------------------------------------------------------------- 
    167 !CDIR PARALLEL DO PRIVATE( zdk1t )  
    168          !                                                ! =============== 
    169          DO jj = 1, jpj                                 ! Horizontal slab 
    170             !                                             ! =============== 
    171             DO ji = 1, jpi   ! vector opt. 
    172                DO jk = mikt(ji,jj), jpkm1 
    173                ! 1. Vertical tracer gradient at level jk and jk+1 
    174                ! ------------------------------------------------ 
    175                ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 
    176                   zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn) - ptb(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1) 
    177                ! 
    178                   IF( jk == mikt(ji,jj) ) THEN  ;   zdkt(ji,jj,jk) = zdk1t(ji,jj,jk) 
    179                   ELSE                          ;   zdkt(ji,jj,jk) = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
    180                   ENDIF 
     175!!!!!!!!!!CDIR PARALLEL DO PRIVATE( zdk1t )  
     176            ! 1. Vertical tracer gradient at level jk and jk+1 
     177            ! ------------------------------------------------ 
     178         !  
     179         ! interior value  
     180         DO jk = 2, jpkm1                
     181            DO jj = 1, jpj 
     182               DO ji = 1, jpi   ! vector opt. 
     183                  zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn  ) - ptb(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1) 
     184                  ! 
     185                  zdkt(ji,jj,jk)  = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn  ) ) * wmask(ji,jj,jk) 
    181186               END DO 
    182187            END DO 
    183188         END DO 
    184  
    185             ! 2. Horizontal fluxes 
    186             ! --------------------    
    187          DO jj = 1 , jpjm1 
    188             DO ji = 1, fs_jpim1   ! vector opt. 
    189                DO jk = mikt(ji,jj), jpkm1 
     189         ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 
     190         zdk1t(:,:,1) = ( ptb(:,:,1,jn  ) - ptb(:,:,2,jn) ) * wmask(:,:,2) 
     191         zdkt (:,:,1) = zdk1t(:,:,1) 
     192         IF ( ln_isfcav ) THEN 
     193            DO jj = 1, jpj 
     194               DO ji = 1, jpi   ! vector opt. 
     195                  ikt = mikt(ji,jj) ! surface level 
     196                  zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn  ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1) 
     197                  zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt) 
     198               END DO 
     199            END DO 
     200         END IF 
     201 
     202         ! 2. Horizontal fluxes 
     203         ! --------------------    
     204         DO jk = 1, jpkm1 
     205            DO jj = 1 , jpjm1 
     206               DO ji = 1, fs_jpim1   ! vector opt. 
    190207                  zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) 
    191208                  zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) 
    … …  
    208225               END DO 
    209226            END DO 
    210          END DO 
    211227 
    212228            ! II.4 Second derivative (divergence) and add to the general trend 
    213229            ! ---------------------------------------------------------------- 
    214          DO jj = 2 , jpjm1 
    215             DO ji = fs_2, fs_jpim1   ! vector opt. 
    216                DO jk = mikt(ji,jj), jpkm1 
    217                   zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
     230            DO jj = 2 , jpjm1 
     231               DO ji = fs_2, fs_jpim1   ! vector opt. 
     232                  zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 
    218233                  ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  ) 
    219234                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
    … …  
    278293            DO jj = 2, jpjm1 
    279294               DO ji = fs_2, fs_jpim1   ! vector opt. 
    280                   zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
     295                  zcoef0 = - fsahtw(ji,jj,jk) * wmask(ji,jj,jk) 
    281296                  ! 
    282297                  zmsku = 1./MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)      & 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90

    r4990 r5120  
    102102               END DO 
    103103            END DO 
    104             IF( ln_zps ) THEN      ! set gradient at partial step level 
     104            IF( ln_zps ) THEN      ! set gradient at partial step level for the last ocean cell 
    105105               DO jj = 1, jpjm1 
    106106                  DO ji = 1, fs_jpim1   ! vector opt. 
    … …  
    116116                        ztv(ji,jj,jk) = zabe2 * pgv(ji,jj,jn) 
    117117                     ENDIF 
    118                       
    119                      ! (ISH) 
     118                  END DO 
     119               END DO 
     120            ENDIF 
     121            ! (ISH) 
     122            IF( ln_zps .AND. ln_isfcav ) THEN      ! set gradient at partial step level for the first ocean cell 
     123                                                   ! into a cavity 
     124               DO jj = 1, jpjm1 
     125                  DO ji = 1, fs_jpim1   ! vector opt. 
    120126                     ! ice shelf level level MAX(2,jk) => only where ice shelf 
    121127                     iku = miku(ji,jj)  

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r4990 r5120  
    99   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps 
    1010   !!             -   !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC 
     11   !!            3.6  !  2014-11  (P. Mathiot) isf melting forcing  
    1112   !!---------------------------------------------------------------------- 
    1213 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90

    r4990 r5120  
    122122            DO jj=1, jpj 
    123123               DO ji=1, jpi 
    124                   zwt(ji,jj,1:mikt(ji,jj)) = 0._wp 
     124                  zwt(ji,jj,1) = 0._wp 
    125125               END DO 
    126126            END DO 
    … …  
    184184            DO jj = 2, jpjm1 
    185185               DO ji = fs_2, fs_jpim1 
    186                   zwt(ji,jj,1:mikt(ji,jj)) = zwd(ji,jj,1:mikt(ji,jj)) 
    187                   DO jk = mikt(ji,jj)+1, jpkm1 
     186                  zwt(ji,jj,1) = zwd(ji,jj,1) 
     187               END DO 
     188            END DO 
     189            DO jk = 2, jpkm1 
     190               DO jj = 2, jpjm1 
     191                  DO ji = fs_2, fs_jpim1 
    188192                     zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 
    189193                  END DO 
    … …  
    196200         DO jj = 2, jpjm1 
    197201            DO ji = fs_2, fs_jpim1 
    198                ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,mikt(ji,jj)) 
    199                ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,mikt(ji,jj)) 
    200                pta(ji,jj,mikt(ji,jj),jn) = ze3tb * ptb(ji,jj,mikt(ji,jj),jn)                     & 
    201                   &                      + p2dt(mikt(ji,jj)) * ze3tn * pta(ji,jj,mikt(ji,jj),jn) 
    202                DO jk = mikt(ji,jj)+1, jpkm1 
     202               ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1) 
     203               ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1) 
     204               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn)                     & 
     205                  &                      + p2dt(1) * ze3tn * pta(ji,jj,1,jn) 
     206            END DO 
     207         END DO 
     208         DO jk = 2, jpkm1 
     209            DO jj = 2, jpjm1 
     210               DO ji = fs_2, fs_jpim1 
    203211                  ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,jk) 
    204212                  ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t  (ji,jj,jk) 
    … …  
    213221            DO ji = fs_2, fs_jpim1 
    214222               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
    215                DO jk = jpk-2, mikt(ji,jj), -1 
     223            END DO 
     224         END DO 
     225         DO jk = jpk-2, 1, -1 
     226            DO jj = 2, jpjm1 
     227               DO ji = fs_2, fs_jpim1 
    216228                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   & 
    217229                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk) 

  • trunk/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r4990 r5120  
    88   !!             -   !  2004-03  (C. Ethe)  adapted for passive tracers 
    99   !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA  
     10   !!            3.6  !  2014-11  (P. Mathiot) Add zps_hde_isf (needed to open a cavity) 
    1011   !!====================================================================== 
    1112    
    … …  
    2728   PRIVATE 
    2829 
    29    PUBLIC   zps_hde    ! routine called by step.F90 
     30   PUBLIC   zps_hde     ! routine called by step.F90 
     31   PUBLIC   zps_hde_isf ! routine called by step.F90 
    3032 
    3133   !! * Substitutions 
    … …  
    4042 
    4143   SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv,   & 
     44      &                          prd, pgru, pgrv    ) 
     45      !!---------------------------------------------------------------------- 
     46      !!                     ***  ROUTINE zps_hde  *** 
     47      !!                     
     48      !! ** Purpose :   Compute the horizontal derivative of T, S and rho 
     49      !!      at u- and v-points with a linear interpolation for z-coordinate 
     50      !!      with partial steps. 
     51      !! 
     52      !! ** Method  :   In z-coord with partial steps, scale factors on last  
     53      !!      levels are different for each grid point, so that T, S and rd  
     54      !!      points are not at the same depth as in z-coord. To have horizontal 
     55      !!      gradients again, we interpolate T and S at the good depth :  
     56      !!      Linear interpolation of T, S    
     57      !!         Computation of di(tb) and dj(tb) by vertical interpolation: 
     58      !!          di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~ 
     59      !!          dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~ 
     60      !!         This formulation computes the two cases: 
     61      !!                 CASE 1                   CASE 2   
     62      !!         k-1  ___ ___________   k-1   ___ ___________ 
     63      !!                    Ti  T~                  T~  Ti+1 
     64      !!                  _____                        _____ 
     65      !!         k        |   |Ti+1     k           Ti |   | 
     66      !!                  |   |____                ____|   | 
     67      !!              ___ |   |   |           ___  |   |   | 
     68      !!                   
     69      !!      case 1->   e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then 
     70      !!          t~ = t(i+1,j  ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1) 
     71      !!        ( t~ = t(i  ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1)  ) 
     72      !!          or 
     73      !!      case 2->   e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then 
     74      !!          t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i ) 
     75      !!        ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) ) 
     76      !!          Idem for di(s) and dj(s)           
     77      !! 
     78      !!      For rho, we call eos which will compute rd~(t~,s~) at the right 
     79      !!      depth zh from interpolated T and S for the different formulations 
     80      !!      of the equation of state (eos). 
     81      !!      Gradient formulation for rho : 
     82      !!          di(rho) = rd~ - rd(i,j,k)   or   rd(i+1,j,k) - rd~ 
     83      !! 
     84      !! ** Action  : compute for top interfaces 
     85      !!              - pgtu, pgtv: horizontal gradient of tracer at u- & v-points 
     86      !!              - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points 
     87      !!---------------------------------------------------------------------- 
     88      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index 
     89      INTEGER                              , INTENT(in   )           ::  kjpt        ! number of tracers 
     90      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields 
     91      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts  
     92      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
     93      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad of prd at u- & v-pts (bottom) 
     94      ! 
     95      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
     96      INTEGER  ::   iku, ikv, ikum1, ikvm1   ! partial step level (ocean bottom level) at u- and v-points 
     97      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars 
     98      REAL(wp), DIMENSION(jpi,jpj)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos 
     99      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::  zti, ztj             !  
     100      !!---------------------------------------------------------------------- 
     101      ! 
     102      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde') 
     103      ! 
     104      pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 
     105      zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 
     106      zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ; 
     107      ! 
     108      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
     109         ! 
     110         DO jj = 1, jpjm1 
     111            DO ji = 1, jpim1 
     112               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
     113               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
     114               ze3wu = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
     115               ze3wv = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
     116               ! 
     117               ! i- direction 
     118               IF( ze3wu >= 0._wp ) THEN      ! case 1 
     119                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku) 
     120                  ! interpolated values of tracers 
     121                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
     122                  ! gradient of  tracers 
     123                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     124               ELSE                           ! case 2 
     125                  zmaxu = -ze3wu / fse3w(ji,jj,iku) 
     126                  ! interpolated values of tracers 
     127                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
     128                  ! gradient of tracers 
     129                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     130               ENDIF 
     131               ! 
     132               ! j- direction 
     133               IF( ze3wv >= 0._wp ) THEN      ! case 1 
     134                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv) 
     135                  ! interpolated values of tracers 
     136                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
     137                  ! gradient of tracers 
     138                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     139               ELSE                           ! case 2 
     140                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv) 
     141                  ! interpolated values of tracers 
     142                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
     143                  ! gradient of tracers 
     144                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     145               ENDIF 
     146            END DO 
     147         END DO 
     148         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     149         ! 
     150      END DO 
     151 
     152      ! horizontal derivative of density anomalies (rd) 
     153      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
     154         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ;  
     155         DO jj = 1, jpjm1 
     156            DO ji = 1, jpim1 
     157               iku = mbku(ji,jj) 
     158               ikv = mbkv(ji,jj) 
     159               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
     160               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
     161               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1 
     162               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2 
     163               ENDIF 
     164               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1 
     165               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2 
     166               ENDIF 
     167            END DO 
     168         END DO 
     169 
     170         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
     171         ! step and store it in  zri, zrj for each  case 
     172         CALL eos( zti, zhi, zri )   
     173         CALL eos( ztj, zhj, zrj ) 
     174 
     175         ! Gradient of density at the last level  
     176         DO jj = 1, jpjm1 
     177            DO ji = 1, jpim1 
     178               iku = mbku(ji,jj) 
     179               ikv = mbkv(ji,jj) 
     180               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
     181               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
     182               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1 
     183               ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2 
     184               ENDIF 
     185               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1 
     186               ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2 
     187               ENDIF 
     188            END DO 
     189         END DO 
     190         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions 
     191         ! 
     192      END IF 
     193      ! 
     194      IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde') 
     195      ! 
     196   END SUBROUTINE zps_hde 
     197   ! 
     198   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv,   & 
    42199      &                          prd, pgru, pgrv, pmru, pmrv, pgzu, pgzv, pge3ru, pge3rv,  & 
    43       &                   sgtu, sgtv, sgru, sgrv, smru, smrv, sgzu, sgzv, sge3ru, sge3rv ) 
     200      &                   pgtui, pgtvi, pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi ) 
    44201      !!---------------------------------------------------------------------- 
    45202      !!                     ***  ROUTINE zps_hde  *** 
    … …  
    82239      !! 
    83240      !! ** Action  : compute for top and bottom interfaces 
    84       !!              - pgtu, pgtv, sgtu, sgtv: horizontal gradient of tracer at u- & v-points 
    85       !!              - pgru, pgrv, sgru, sgtv: horizontal gradient of rho (if present) at u- & v-points 
    86       !!              - pmru, pmrv, smru, smrv: horizontal sum of rho at u- & v- point (used in dynhpg with vvl) 
    87       !!              - pgzu, pgzv, sgzu, sgzv: horizontal gradient of z at u- and v- point (used in dynhpg with vvl) 
    88       !!              - pge3ru, pge3rv, sge3ru, sge3rv: horizontal gradient of rho weighted by local e3w at u- & v-points  
     241      !!              - pgtu, pgtv, pgtui, pgtvi: horizontal gradient of tracer at u- & v-points 
     242      !!              - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points 
     243      !!              - pmru, pmrv, pmrui, pmrvi: horizontal sum of rho at u- & v- point (used in dynhpg with vvl) 
     244      !!              - pgzu, pgzv, pgzui, pgzvi: horizontal gradient of z at u- and v- point (used in dynhpg with vvl) 
     245      !!              - pge3ru, pge3rv, pge3rui, pge3rvi: horizontal gradient of rho weighted by local e3w at u- & v-points  
    89246      !!---------------------------------------------------------------------- 
    90247      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index 
    … …  
    92249      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields 
    93250      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts  
    94       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  sgtu, sgtv  ! hor. grad. of stra at u- & v-pts (ISF) 
     251      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui, pgtvi  ! hor. grad. of stra at u- & v-pts (ISF) 
    95252      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
    96253      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv      ! hor. grad of prd at u- & v-pts (bottom) 
    … …  
    98255      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu, pgzv      ! hor. grad of z   at u- & v-pts (bottom) 
    99256      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru, pge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (bottom) 
    100       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  sgru, sgrv      ! hor. grad of prd at u- & v-pts (top) 
    101       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  smru, smrv      ! hor. sum  of prd at u- & v-pts (top) 
    102       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  sgzu, sgzv      ! hor. grad of z   at u- & v-pts (top) 
    103       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  sge3ru, sge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (top) 
     257      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui, pgrvi      ! hor. grad of prd at u- & v-pts (top) 
     258      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmrui, pmrvi      ! hor. sum  of prd at u- & v-pts (top) 
     259      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzui, pgzvi      ! hor. grad of z   at u- & v-pts (top) 
     260      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3rui, pge3rvi  ! hor. grad of prd weighted by local e3w at u- & v-pts (top) 
    104261      ! 
    105262      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
    … …  
    110267      !!---------------------------------------------------------------------- 
    111268      ! 
    112       IF( nn_timing == 1 )  CALL timing_start( 'zps_hde') 
     269      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_isf') 
    113270      ! 
    114271      pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 
    115       sgtu(:,:,:)=0.0_wp ; sgtv(:,:,:)=0.0_wp ; 
     272      pgtui(:,:,:)=0.0_wp ; pgtvi(:,:,:)=0.0_wp ; 
    116273      zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ; 
    117274      zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ; 
    … …  
    256413                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) ) 
    257414                  ! gradient of tracers 
    258                   sgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     415                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
    259416               ELSE                           ! case 2 
    260417                  zmaxu = - ze3wu / fse3w(ji,jj,iku+1) 
    … …  
    262419                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) ) 
    263420                  ! gradient of  tracers 
    264                   sgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     421                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
    265422               ENDIF 
    266423               ! 
    … …  
    271428                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) ) 
    272429                  ! gradient of tracers 
    273                   sgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     430                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
    274431               ELSE                           ! case 2 
    275432                  zmaxv =  - ze3wv / fse3w(ji,jj,ikv+1) 
    … …  
    277434                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) ) 
    278435                  ! gradient of tracers 
    279                   sgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     436                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
    280437               ENDIF 
    281438            END DO!! 
    282439         END DO!! 
    283          CALL lbc_lnk( sgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( sgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     440         CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
    284441         ! 
    285442      END DO 
    … …  
    287444      ! horizontal derivative of density anomalies (rd) 
    288445      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
    289          sgru(:,:)  =0.0_wp ; sgrv(:,:)  =0.0_wp ; 
    290          sgzu(:,:)  =0.0_wp ; sgzv(:,:)  =0.0_wp ; 
    291          smru(:,:)  =0.0_wp ; smru(:,:)  =0.0_wp ; 
    292          sge3ru(:,:)=0.0_wp ; sge3rv(:,:)=0.0_wp ; 
     446         pgrui(:,:)  =0.0_wp ; pgrvi(:,:)  =0.0_wp ; 
     447         pgzui(:,:)  =0.0_wp ; pgzvi(:,:)  =0.0_wp ; 
     448         pmrui(:,:)  =0.0_wp ; pmrui(:,:)  =0.0_wp ; 
     449         pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ; 
    293450 
    294451         DO jj = 1, jpjm1 
    … …  
    321478               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 
    322479               IF( ze3wu >= 0._wp ) THEN 
    323                  sgzu  (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku) 
    324                  sgru  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1 
    325                  smru  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1  
    326                  sge3ru(ji,jj) = umask(ji,jj,iku+1)                                                                  & 
     480                 pgzui  (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku) 
     481                 pgrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1 
     482                 pmrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1  
     483                 pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                  & 
    327484                                * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   & 
    328485                                   - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1 
    329486               ELSE 
    330                  sgzu  (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu) 
    331                  sgru  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) - zri(ji,jj) )      ! i: 2 
    332                  smru  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2 
    333                  sge3ru(ji,jj) = umask(ji,jj,iku+1)                                                                   & 
     487                 pgzui  (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu) 
     488                 pgrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) - zri(ji,jj) )      ! i: 2 
     489                 pmrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2 
     490                 pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                   & 
    334491                                * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  & 
    335492                                   -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2 
    336493               ENDIF 
    337494               IF( ze3wv >= 0._wp ) THEN 
    338                  sgzv  (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv)  
    339                  sgrv  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )        ! j: 1 
    340                  smrv  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1 
    341                  sge3rv(ji,jj) = vmask(ji,jj,ikv+1)                                                                  &  
     495                 pgzvi  (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv)  
     496                 pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )        ! j: 1 
     497                 pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1 
     498                 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                  &  
    342499                                * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  & 
    343500                                   - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1 
    344501                                  ! + 2 due to the formulation in density and not in anomalie in hpg sco 
    345502               ELSE 
    346                  sgzv  (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv) 
    347                  sgrv  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )     ! j: 2 
    348                  smrv  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2 
    349                  sge3rv(ji,jj) = vmask(ji,jj,ikv+1)                                                                   & 
     503                 pgzvi  (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv) 
     504                 pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )     ! j: 2 
     505                 pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2 
     506                 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                   & 
    350507                                * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) & 
    351508                                   -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2 
    … …  
    353510            END DO 
    354511         END DO 
    355          CALL lbc_lnk( sgru   , 'U', -1. )   ;   CALL lbc_lnk( sgrv   , 'V', -1. )   ! Lateral boundary conditions 
    356          CALL lbc_lnk( smru   , 'U',  1. )   ;   CALL lbc_lnk( smrv   , 'V',  1. )   ! Lateral boundary conditions 
    357          CALL lbc_lnk( sgzu   , 'U', -1. )   ;   CALL lbc_lnk( sgzv   , 'V', -1. )   ! Lateral boundary conditions 
    358          CALL lbc_lnk( sge3ru , 'U', -1. )   ;   CALL lbc_lnk( sge3rv , 'V', -1. )   ! Lateral boundary conditions 
     512         CALL lbc_lnk( pgrui   , 'U', -1. )   ;   CALL lbc_lnk( pgrvi   , 'V', -1. )   ! Lateral boundary conditions 
     513         CALL lbc_lnk( pmrui   , 'U',  1. )   ;   CALL lbc_lnk( pmrvi   , 'V',  1. )   ! Lateral boundary conditions 
     514         CALL lbc_lnk( pgzui   , 'U', -1. )   ;   CALL lbc_lnk( pgzvi   , 'V', -1. )   ! Lateral boundary conditions 
     515         CALL lbc_lnk( pge3rui , 'U', -1. )   ;   CALL lbc_lnk( pge3rvi , 'V', -1. )   ! Lateral boundary conditions 
    359516         ! 
    360517      END IF   
    361518      ! 
    362       IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde') 
    363       ! 
    364    END SUBROUTINE zps_hde 
    365  
     519      IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde_isf') 
     520      ! 
     521   END SUBROUTINE zps_hde_isf 
    366522   !!====================================================================== 
    367523END MODULE zpshde 

  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r4990 r5120  
    120120                  zbfrt(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp) 
    121121                  zbfrt(ji,jj) = MIN(zbfrt(ji,jj), rn_bfri2_max) 
    122 ! (ISF) 
    123                   ikbt = mikt(ji,jj) 
    124 ! JC: possible WAD implementation should modify line below if layers vanish 
    125                   ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp 
    126                   ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 
    127                   ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max) 
    128  
    129122               END DO 
    130123            END DO 
     124! (ISF) 
     125            IF ( ln_isfcav ) THEN 
     126               DO jj = 1, jpj 
     127                  DO ji = 1, jpi 
     128                     ikbt = mikt(ji,jj) 
     129! JC: possible WAD implementation should modify line below if layers vanish 
     130                     ztmp = (1-tmask(ji,jj,1)) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp 
     131                     ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 
     132                     ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max) 
     133                  END DO 
     134               END DO 
     135            END IF 
    131136         !    
    132137         ELSE 
    … …  
    152157               ! 
    153158               ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 
    154                IF ( miku(ji,jj) + 2 .GE. mbku(ji,jj) ) THEN 
    155                   bfrua(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) )   & 
    156                                &            + ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) ) & 
    157                                &          * zecu * (1._wp - umask(ji,jj,1)) 
    158                END IF 
    159                IF ( mikv(ji,jj) + 2 .GE. mbkv(ji,jj) ) THEN 
    160                   bfrva(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) )   & 
    161                                &            + ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) ) & 
    162                                &          * zecv * (1._wp - vmask(ji,jj,1)) 
    163                END IF 
    164                ! (ISF) ======================================================================== 
    165                ikbu = miku(ji,jj)         ! ocean bottom level at u- and v-points  
    166                ikbv = mikv(ji,jj)         ! (deepest ocean u- and v-points) 
    167                ! 
    168                zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     & 
    169                   &           + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  ) 
    170                zuv  = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     & 
    171                   &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  ) 
    172                ! 
    173                zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 ) 
    174                zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 ) 
    175                ! 
    176                tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) * zecu * (1._wp - umask(ji,jj,1)) 
    177                tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) * zecv * (1._wp - vmask(ji,jj,1)) 
    178                ! (ISF) END ==================================================================== 
    179                ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 
    180                IF ( miku(ji,jj) + 2 .GE. mbku(ji,jj) ) THEN 
    181                   tfrua(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) )   & 
    182                                &            + ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) ) & 
    183                                &          * zecu * (1._wp - umask(ji,jj,1)) 
    184                END IF 
    185                IF ( mikv(ji,jj) + 2 .GE. mbkv(ji,jj) ) THEN 
    186                   tfrva(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) )   & 
    187                                &            + ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) ) & 
    188                                &          * zecv * (1._wp - vmask(ji,jj,1)) 
     159               IF ( ln_isfcav ) THEN 
     160                  IF ( miku(ji,jj) + 1 .GE. mbku(ji,jj) ) THEN 
     161                     bfrua(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) )   & 
     162                                  &            + ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) ) & 
     163                                  &          * zecu * (1._wp - umask(ji,jj,1)) 
     164                  END IF 
     165                  IF ( mikv(ji,jj) + 1 .GE. mbkv(ji,jj) ) THEN 
     166                     bfrva(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) )   & 
     167                                  &            + ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) ) & 
     168                                  &          * zecv * (1._wp - vmask(ji,jj,1)) 
     169                  END IF 
    189170               END IF 
    190171            END DO 
    191172         END DO 
     173         IF ( ln_isfcav ) THEN 
     174            DO jj = 2, jpjm1 
     175               DO ji = 2, jpim1 
     176                  ! (ISF) ======================================================================== 
     177                  ikbu = miku(ji,jj)         ! ocean bottom level at u- and v-points  
     178                  ikbv = mikv(ji,jj)         ! (deepest ocean u- and v-points) 
     179                  ! 
     180                  zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     & 
     181                     &           + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  ) 
     182                  zuv  = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     & 
     183                     &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  ) 
     184              ! 
     185                  zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 ) 
     186                  zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 ) 
     187              ! 
     188                  tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) * zecu * (1._wp - umask(ji,jj,1)) 
     189                  tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) * zecv * (1._wp - vmask(ji,jj,1)) 
     190              ! (ISF) END ==================================================================== 
     191              ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 
     192                  IF ( miku(ji,jj) + 1 .GE. mbku(ji,jj) ) THEN 
     193                     tfrua(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) )   & 
     194                                  &            + ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) ) & 
     195                                  &          * zecu * (1._wp - umask(ji,jj,1)) 
     196                  END IF 
     197                  IF ( mikv(ji,jj) + 1 .GE. mbkv(ji,jj) ) THEN 
     198                     tfrva(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) )   & 
     199                                  &            + ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) ) & 
     200                                  &          * zecv * (1._wp - vmask(ji,jj,1)) 
     201                  END IF 
     202               END DO 
     203            END DO 
     204         END IF 
    192205         ! 
    193206         CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition 

  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

    r4990 r5120  
    156156         END DO 
    157157         ! mask zmsk in order to have avt and avs masked 
    158          zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk) 
     158         zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk) 
    159159 
    160160 
    … …  
    191191               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    & 
    192192                  &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   & 
    193                   &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * umask(ji,jj,jk) 
     193                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * wumask(ji,jj,jk) 
    194194               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    & 
    195195                  &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   & 
    196                   &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * vmask(ji,jj,jk) 
     196                  &                  avs(ji,jj,jk), avs(ji,jj+1,jk) )  * wvmask(ji,jj,jk) 
    197197            END DO 
    198198         END DO 
    … …  
    255255      IF( zdf_ddm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_ddm_init : unable to allocate arrays' ) 
    256256      !                               ! initialization to masked Kz 
    257       avs(:,:,:) = rn_avt0 * tmask(:,:,:)  
     257      avs(:,:,:) = rn_avt0 * wmask(:,:,:)  
    258258      ! 
    259259   END SUBROUTINE zdf_ddm_init 

  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90

    r4990 r5120  
    1414   !!---------------------------------------------------------------------- 
    1515   USE par_oce         ! mesh and scale factors 
    16    USE sbc_oce         ! surface module (only for nn_isf in the option compatibility test) 
    1716   USE ldftra_oce      ! ocean active tracers: lateral physics 
    1817   USE ldfdyn_oce      ! ocean dynamics lateral physics 
    … …  
    118117      IF( ioptio == 0 .OR. ioptio > 1 .AND. .NOT. lk_esopa )   & 
    119118         &   CALL ctl_stop( ' one and only one vertical diffusion option has to be defined ' ) 
    120       IF( ( lk_zdfric .OR. lk_zdfgls .OR. lk_zdfkpp ) .AND. nn_isf .NE. 0 )   & 
     119      IF( ( lk_zdfric .OR. lk_zdfgls .OR. lk_zdfkpp ) .AND. ln_isfcav )   & 
    121120         &   CALL ctl_stop( ' only zdfcst and zdftke were tested with ice shelves cavities ' ) 
    122121      ! 

  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r5112 r5120  
    2626   !!                 !                                + cleaning of the parameters + bugs correction 
    2727   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase 
     28   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability 
    2829   !!---------------------------------------------------------------------- 
    2930#if defined key_zdftke   ||   defined key_esopa 
    … …  
    236237      zfact3 = 0.5_wp       * rn_ediss 
    237238      ! 
     239      ! 
    238240      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    239241      !                     !  Surface boundary condition on tke 
    240242      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
     243      IF ( ln_isfcav ) THEN 
     244         DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin 
     245            DO ji = fs_2, fs_jpim1   ! vector opt. 
     246               en(ji,jj,mikt(ji,jj))=rn_emin * tmask(ji,jj,1) 
     247            END DO 
     248         END DO 
     249      END IF 
    241250      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
    242251         DO ji = fs_2, fs_jpim1   ! vector opt. 
    243             IF (mikt(ji,jj) .GT. 1) THEN 
    244                en(ji,jj,mikt(ji,jj))=rn_emin 
    245             ELSE 
    246                en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
    247             END IF 
     252            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
    248253         END DO 
    249254      END DO 
    … …  
    301306         END DO 
    302307         zcof = 0.016 / SQRT( zrhoa * zcdrag ) 
     308!CDIR NOVERRCHK 
    303309         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en 
    304             DO jj = 2, jpjm1 
     310!CDIR NOVERRCHK 
     311            DO jj = 2, jpjm1 
     312!CDIR NOVERRCHK 
    305313               DO ji = fs_2, fs_jpim1   ! vector opt. 
    306314                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift 
    … …  
    309317                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 
    310318                  !                                           ! TKE Langmuir circulation source term 
    311                   en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk) 
     319                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    312320               END DO 
    313321            END DO 
    … …  
    328336               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   & 
    329337                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   &  
    330                   &           / (  fse3uw_n(ji,jj,jk)         & 
    331                   &              * fse3uw_b(ji,jj,jk) ) 
     338                  &                            / (  fse3uw_n(ji,jj,jk)               & 
     339                  &                              *  fse3uw_b(ji,jj,jk) ) 
    332340               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   & 
    333341                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   & 
    … …  
    338346      END DO 
    339347      ! 
    340       DO jj = 2, jpjm1 
    341          DO ji = fs_2, fs_jpim1   ! vector opt. 
    342             DO jk = mikt(ji,jj)+1, jpkm1           !* Matrix and right hand side in en 
     348      DO jk = 2, jpkm1           !* Matrix and right hand side in en 
     349         DO jj = 2, jpjm1 
     350            DO ji = fs_2, fs_jpim1   ! vector opt. 
    343351               zcof   = zfact1 * tmask(ji,jj,jk) 
    344352               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal 
    … …  
    357365               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    & 
    358366                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) & 
    359                   &                                 * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
    360             END DO 
    361             !                          !* Matrix inversion from level 2 (tke prescribed at level 1) 
    362             DO jk = mikt(ji,jj)+2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
     367                  &                                 * wmask(ji,jj,jk) 
     368            END DO 
     369         END DO 
     370      END DO 
     371      !                          !* Matrix inversion from level 2 (tke prescribed at level 1) 
     372      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
     373         DO jj = 2, jpjm1 
     374            DO ji = fs_2, fs_jpim1    ! vector opt. 
    363375               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
    364376            END DO 
    365             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    366             zd_lw(ji,jj,mikt(ji,jj)+1) = en(ji,jj,mikt(ji,jj)+1) - zd_lw(ji,jj,mikt(ji,jj)+1) * en(ji,jj,mikt(ji,jj))    ! Surface boudary conditions on tke 
    367             ! 
    368             DO jk = mikt(ji,jj)+2, jpkm1 
     377         END DO 
     378      END DO 
     379      ! 
     380      ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
     381      DO jj = 2, jpjm1 
     382         DO ji = fs_2, fs_jpim1   ! vector opt. 
     383            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke 
     384         END DO 
     385      END DO 
     386      DO jk = 3, jpkm1 
     387         DO jj = 2, jpjm1 
     388            DO ji = fs_2, fs_jpim1    ! vector opt. 
    369389               zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1) 
    370390            END DO 
    371             ! 
    372             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
     391         END DO 
     392      END DO 
     393      ! 
     394      ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
     395      DO jj = 2, jpjm1 
     396         DO ji = fs_2, fs_jpim1   ! vector opt. 
    373397            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 
    374             ! 
    375             DO jk = jpk-2, mikt(ji,jj)+1, -1 
     398         END DO 
     399      END DO 
     400      DO jk = jpk-2, 2, -1 
     401         DO jj = 2, jpjm1 
     402            DO ji = fs_2, fs_jpim1    ! vector opt. 
    376403               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
    377404            END DO 
    378             ! 
    379             DO jk = mikt(ji,jj), jpkm1                             ! set the minimum value of tke 
    380                en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk) 
     405         END DO 
     406      END DO 
     407      DO jk = 2, jpkm1                             ! set the minimum value of tke 
     408         DO jj = 2, jpjm1 
     409            DO ji = fs_2, fs_jpim1   ! vector opt. 
     410               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 
    381411            END DO 
    382412         END DO 
    … …  
    391421               DO ji = fs_2, fs_jpim1   ! vector opt. 
    392422                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    393                      &                                 * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) * tmask(ji,jj,1) 
     423                     &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    394424               END DO 
    395425            END DO 
    … …  
    400430               jk = nmln(ji,jj) 
    401431               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    402                   &                                 * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) * tmask(ji,jj,1) 
     432                  &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    403433            END DO 
    404434         END DO 
    … …  
    416446                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    417447                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    418                      &                        * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * tmask(ji,jj,1) 
     448                     &                        * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    419449               END DO 
    420450            END DO 
    … …  
    484514      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2) 
    485515      ! 
     516      ! initialisation of interior minimum value (avoid a 2d loop with mikt) 
     517      zmxlm(:,:,:)  = rmxl_min     
     518      zmxld(:,:,:)  = rmxl_min 
     519      ! 
    486520      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 
    487521         DO jj = 2, jpjm1 
    488522            DO ji = fs_2, fs_jpim1 
    489                IF (mikt(ji,jj) .GT. 1) THEN 
    490                   zmxlm(ji,jj,mikt(ji,jj)) = rmxl_min 
    491                ELSE 
    492                   zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
    493                   zmxlm(ji,jj,mikt(ji,jj)) = MAX( rn_mxl0, zraug * taum(ji,jj) ) 
    494                END IF 
     523               zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
     524               zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 
    495525            END DO 
    496526         END DO 
    497527      ELSE  
    498          DO jj = 2, jpjm1 
    499             DO ji = fs_2, fs_jpim1                         ! surface set to the minimum value 
    500                zmxlm(ji,jj,mikt(ji,jj)) = MAX( tmask(ji,jj,1) * rn_mxl0, rmxl_min) 
    501             END DO 
    502          END DO 
     528         zmxlm(:,:,1) = rn_mxl0 
    503529      ENDIF 
    504       zmxlm(:,:,jpk)  = rmxl_min     ! last level set to the interior minium value 
    505       ! 
    506 !CDIR NOVERRCHK 
    507       DO jj = 2, jpjm1 
    508 !CDIR NOVERRCHK 
    509          DO ji = fs_2, fs_jpim1   ! vector opt. 
    510             !CDIR NOVERRCHK 
    511             DO jk = mikt(ji,jj)+1, jpkm1              ! interior value : l=sqrt(2*e/n^2) 
     530      ! 
     531!CDIR NOVERRCHK 
     532      DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2) 
     533!CDIR NOVERRCHK 
     534         DO jj = 2, jpjm1 
     535!CDIR NOVERRCHK 
     536            DO ji = fs_2, fs_jpim1   ! vector opt. 
    512537               zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    513                zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  ) 
    514             END DO 
    515             zmxld(ji,jj,mikt(ji,jj)) = zmxlm(ji,jj,mikt(ji,jj))   ! surface set to the minimum value  
     538               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 
     539            END DO 
    516540         END DO 
    517541      END DO 
    … …  
    519543      !                     !* Physical limits for the mixing length 
    520544      ! 
    521       zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the zmxlm   value 
     545      zmxld(:,:,1  ) = zmxlm(:,:,1)   ! surface set to the minimum value  
    522546      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value 
    523547      ! 
    524548      SELECT CASE ( nn_mxl ) 
    525549      ! 
     550      ! where wmask = 0 set zmxlm == fse3w 
    526551      CASE ( 0 )           ! bounded by the distance to surface and bottom 
    527          DO jj = 2, jpjm1 
    528             DO ji = fs_2, fs_jpim1   ! vector opt. 
    529                DO jk = mikt(ji,jj)+1, jpkm1 
     552         DO jk = 2, jpkm1 
     553            DO jj = 2, jpjm1 
     554               DO ji = fs_2, fs_jpim1   ! vector opt. 
    530555                  zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   & 
    531556                  &            fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 
    532                   zmxlm(ji,jj,jk) = zemxl 
    533                   zmxld(ji,jj,jk) = zemxl 
     557                  ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 
     558                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 
     559                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 
    534560               END DO 
    535561            END DO 
    … …  
    537563         ! 
    538564      CASE ( 1 )           ! bounded by the vertical scale factor 
    539          DO jj = 2, jpjm1 
    540             DO ji = fs_2, fs_jpim1   ! vector opt. 
    541                DO jk = mikt(ji,jj)+1, jpkm1 
     565         DO jk = 2, jpkm1 
     566            DO jj = 2, jpjm1 
     567               DO ji = fs_2, fs_jpim1   ! vector opt. 
    542568                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 
    543569                  zmxlm(ji,jj,jk) = zemxl 
    … …  
    548574         ! 
    549575      CASE ( 2 )           ! |dk[xml]| bounded by e3t : 
    550          DO jj = 2, jpjm1 
    551             DO ji = fs_2, fs_jpim1   ! vector opt. 
    552                DO jk = mikt(ji,jj)+1, jpkm1         ! from the surface to the bottom : 
     576         DO jk = 2, jpkm1         ! from the surface to the bottom : 
     577            DO jj = 2, jpjm1 
     578               DO ji = fs_2, fs_jpim1   ! vector opt. 
    553579                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    554580               END DO 
    555                DO jk = jpkm1, mikt(ji,jj)+1, -1     ! from the bottom to the surface : 
     581            END DO 
     582         END DO 
     583         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : 
     584            DO jj = 2, jpjm1 
     585               DO ji = fs_2, fs_jpim1   ! vector opt. 
    556586                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    557587                  zmxlm(ji,jj,jk) = zemxl 
    … …  
    562592         ! 
    563593      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t : 
    564          DO jj = 2, jpjm1 
    565             DO ji = fs_2, fs_jpim1   ! vector opt. 
    566                DO jk = mikt(ji,jj)+1, jpkm1         ! from the surface to the bottom : lup 
     594         DO jk = 2, jpkm1         ! from the surface to the bottom : lup 
     595            DO jj = 2, jpjm1 
     596               DO ji = fs_2, fs_jpim1   ! vector opt. 
    567597                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    568598               END DO 
    569                DO jk = jpkm1, mikt(ji,jj)+1, -1     ! from the bottom to the surface : ldown 
     599            END DO 
     600         END DO 
     601         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown 
     602            DO jj = 2, jpjm1 
     603               DO ji = fs_2, fs_jpim1   ! vector opt. 
    570604                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    571605               END DO 
    … …  
    604638               zsqen = SQRT( en(ji,jj,jk) ) 
    605639               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen 
    606                avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * tmask(ji,jj,jk) 
    607                avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
     640               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk) 
     641               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
    608642               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 
    609643            END DO 
    … …  
    612646      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
    613647      ! 
    614       DO jj = 2, jpjm1 
    615          DO ji = fs_2, fs_jpim1   ! vector opt. 
    616             DO jk = miku(ji,jj)+1, jpkm1            !* vertical eddy viscosity at u- and v-points 
    617                avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk) 
    618             END DO 
    619             DO jk = mikv(ji,jj)+1, jpkm1 
    620                avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk) 
     648      DO jk = 2, jpkm1            !* vertical eddy viscosity at wu- and wv-points 
     649         DO jj = 2, jpjm1 
     650            DO ji = fs_2, fs_jpim1   ! vector opt. 
     651               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk) 
     652               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk) 
    621653            END DO 
    622654         END DO 
    … …  
    625657      ! 
    626658      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt 
    627          DO jj = 2, jpjm1 
    628             DO ji = fs_2, fs_jpim1   ! vector opt. 
    629                DO jk = mikt(ji,jj)+1, jpkm1 
     659         DO jk = 2, jpkm1 
     660            DO jj = 2, jpjm1 
     661               DO ji = fs_2, fs_jpim1   ! vector opt. 
    630662                  zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 
    631663                  !                                          ! shear 
    … …  
    639671!!gm and even better with the use of the "true" ri_crit=0.22222...  (this change the results!) 
    640672!!gm              zpdlr = MAX(  0.1_wp,  ri_crit / MAX( ri_crit , zri )  ) 
    641                   avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
     673                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 
    642674# if defined key_c1d 
    643                   e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number 
    644                   e_ric(ji,jj,jk) = zri   * tmask(ji,jj,jk)  ! c1d config. : save Ri 
     675                  e_pdl(ji,jj,jk) = zpdlr * wmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number 
     676                  e_ric(ji,jj,jk) = zri   * wmask(ji,jj,jk)  ! c1d config. : save Ri 
    645677# endif 
    646678              END DO 
    … …  
    749781      !                               !* set vertical eddy coef. to the background value 
    750782      DO jk = 1, jpk 
    751          avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) 
    752          avm (:,:,jk) = avmb(jk) * tmask(:,:,jk) 
    753          avmu(:,:,jk) = avmb(jk) * umask(:,:,jk) 
    754          avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 
     783         avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 
     784         avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 
     785         avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 
     786         avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 
    755787      END DO 
    756788      dissl(:,:,:) = 1.e-12_wp 
    … …  
    814846           en(:,:,:) = rn_emin * tmask(:,:,:) 
    815847           DO jk = 1, jpk                           ! set the Kz to the background value 
    816               avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) 
    817               avm (:,:,jk) = avmb(jk) * tmask(:,:,jk) 
    818               avmu(:,:,jk) = avmb(jk) * umask(:,:,jk) 
    819               avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 
     848              avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 
     849              avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 
     850              avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 
     851              avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 
    820852           END DO 
    821853        ENDIF 

  • trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90

    r5021 r5120  
    126126      zkz(:,:) = 0.e0               !* Associated potential energy consummed over the whole water column 
    127127      DO jk = 2, jpkm1 
    128          zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk)* tmask(:,:,jk) * tmask(:,:,jk-1) 
     128         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk) * wmask(:,:,jk) 
    129129      END DO 
    130130 
    … …  
    135135      END DO 
    136136 
    137       DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 
    138          DO ji = 1, jpi 
    139             DO jk = mikt(ji,jj)+1, jpkm1     !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s 
    140                zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. )   !kz max = 300 cm2/s 
     137      DO jk = 2, jpkm1     !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s 
     138         DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 
     139            DO ji = 1, jpi 
     140               zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk)  !kz max = 300 cm2/s 
    141141            END DO 
    142142         END DO 
    … …  
    166166      !                          !   Update  mixing coefs  !                           
    167167      !                          ! ----------------------- ! 
    168       DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 
    169          DO ji = 1, jpi 
    170             DO jk = mikt(ji,jj)+1, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
    171                avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk) 
    172                avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk) 
     168      DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
     169         DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 
     170            DO ji = 1, jpi 
     171               avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk) 
     172               avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk) 
    173173            END DO 
    174174         END DO 
    175175      END DO 
    176176       
    177       DO jj = 2, jpjm1 
    178          DO ji = fs_2, fs_jpim1  ! vector opt. 
    179             DO jk = mikt(ji,jj)+1, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
    180                avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj  ,jk) ) * umask(ji,jj,jk) 
    181                avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji  ,jj+1,jk) ) * vmask(ji,jj,jk) 
     177      DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
     178         DO jj = 2, jpjm1 
     179            DO ji = fs_2, fs_jpim1  ! vector opt. 
     180               avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj  ,jk) ) * wumask(ji,jj,jk) 
     181               avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk) 
    182182            END DO 
    183183         END DO 
    … …  
    457457         ztpc = 0.e0 
    458458         zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:) 
    459          DO jj = 1, jpj 
    460             DO ji = 1, jpi 
    461                DO jk= mikt(ji,jj)+1, jpkm1 
    462                   ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 
     459         DO jk= 2, jpkm1 
     460            DO jj = 1, jpj 
     461               DO ji = 1, jpi 
     462                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
    463463               END DO 
    464464            END DO 
    … …  
    473473         zav_tide(:,:,:) = MIN( zav_tide(:,:,:), 60.e-4 )    
    474474         zkz(:,:) = 0.e0 
    475          DO jj = 1, jpj 
    476             DO ji = 1, jpi 
    477                DO jk = mikt(ji,jj)+1, jpkm1 
    478                zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zav_tide(ji,jj,jk)* tmask(ji,jj,jk) 
     475         DO jk = 2, jpkm1 
     476            DO jj = 1, jpj 
     477               DO ji = 1, jpi 
     478                  zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zav_tide(ji,jj,jk)* wmask(ji,jj,jk) 
    479479               END DO 
    480480            END DO 
    … …  
    498498         WRITE(numout,*) '          Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) ) 
    499499 
    500          DO jj = 1, jpj 
    501             DO ji = 1, jpi 
    502                DO jk = mikt(ji,jj)+1, jpkm1 
    503                   zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. )   !kz max = 300 cm2/s 
     500         DO jk = 2, jpkm1 
     501            DO jj = 1, jpj 
     502               DO ji = 1, jpi 
     503                  zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk)  !kz max = 300 cm2/s 
    504504               END DO 
    505505            END DO 
    … …  
    510510            DO jj = 1, jpj 
    511511               DO ji = 1, jpi 
    512                   ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 
     512                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
    513513               END DO 
    514514            END DO 
    … …  
    519519         DO jk = 1, jpk 
    520520            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zav_tide(:,:,jk)     * tmask_i(:,:) )   & 
    521                &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) ) 
     521               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) ) 
    522522            ztpc = 1.E50 
    523523            DO jj = 1, jpj 
    … …  
    540540            END DO 
    541541            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zkz(:,:)     * tmask_i(:,:) )   & 
    542                &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) ) 
     542               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) ) 
    543543            WRITE(numout,*) '                jk= ', jk,'   ', ze_z * 1.e4,' cm2/s' 
    544544         END DO 
    … …  
    546546            zkz(:,:) = az_tmx(:,:,jk) /rn_n2min 
    547547            ze_z =                  SUM( e1t(:,:) * e2t(:,:) * zkz(:,:)     * tmask_i(:,:) )   & 
    548                &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) ) 
     548               &     / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) ) 
    549549            WRITE(numout,*)  
    550550            WRITE(numout,*) '          N2 min - jk= ', jk,'   ', ze_z * 1.e4,' cm2/s min= ',MINVAL(zkz)*1.e4,   & 

  • trunk/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r5118 r5120  
    342342         WRITE(numout,*) '                       NEMO team' 
    343343         WRITE(numout,*) '            Ocean General Circulation Model' 
    344          WRITE(numout,*) '                  version 3.4  (2011) ' 
     344         WRITE(numout,*) '                  version 3.6  (2015) ' 
    345345         WRITE(numout,*) 
    346346         WRITE(numout,*) 

  • trunk/NEMOGCM/NEMO/OPA_SRC/step.F90

    r5012 r5120  
    122122      IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP closure scheme for Kz 
    123123      IF( lk_zdfcst  ) THEN                              ! Constant Kz (reset avt, avm[uv] to the background value) 
    124          avt (:,:,:) = rn_avt0 * tmask(:,:,:) 
    125          avmu(:,:,:) = rn_avm0 * umask(:,:,:) 
    126          avmv(:,:,:) = rn_avm0 * vmask(:,:,:) 
     124         avt (:,:,:) = rn_avt0 * wmask (:,:,:) 
     125         avmu(:,:,:) = rn_avm0 * wumask(:,:,:) 
     126         avmv(:,:,:) = rn_avm0 * wvmask(:,:,:) 
    127127      ENDIF 
    128128      IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths 
    … …  
    145145      ! 
    146146      IF( lk_ldfslp ) THEN                            ! slope of lateral mixing 
    147                          CALL eos( tsb, rhd, gdept_0(:,:,:) )             ! before in situ density 
    148          IF( ln_zps )    CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
    149             &                                      rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             ! 
    150             &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
     147                         CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density 
     148         IF( ln_zps .AND. .NOT. ln_isfcav)                               & 
     149            &            CALL zps_hde    ( kstp, jpts, tsb, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
     150            &                                          rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
     151         IF( ln_zps .AND.       ln_isfcav)                               & 
     152            &            CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
     153            &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     154            &                                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level 
    151155         IF( ln_traldf_grif ) THEN                           ! before slope for Griffies operator 
    152156                         CALL ldf_slp_grif( kstp ) 
    … …  
    177181          ! is necessary to compute momentum advection for the rhs of barotropic loop: 
    178182                            CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 
    179             IF( ln_zps )    CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
    180                &                                      rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             ! 
    181                &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
     183            IF( ln_zps .AND. .NOT. ln_isfcav)                               & 
     184               &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
     185               &                                          rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
     186            IF( ln_zps .AND.       ln_isfcav)                               & 
     187               &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
     188               &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     189               &                                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
    182190 
    183191                                  ua(:,:,:) = 0.e0             ! set dynamics trends to zero 
    … …  
    253261                             CALL tra_nxt( kstp )                ! tracer fields at next time step 
    254262                             CALL eos    ( tsa, rhd, rhop, fsdept_n(:,:,:) )  ! Time-filtered in situ density for hpg computation 
    255          IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
    256             &                                      rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             ! 
    257             &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
     263            IF( ln_zps .AND. .NOT. ln_isfcav)                                & 
     264               &             CALL zps_hde    ( kstp, jpts, tsa, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
     265               &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
     266            IF( ln_zps .AND.       ln_isfcav)                                & 
     267               &             CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
     268               &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     269               &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
    258270      ELSE                                                  ! centered hpg  (eos then time stepping) 
    259271         IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case 
    260272                             CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation 
    261          IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
    262          &                                      rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             ! 
    263          &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
     273         IF( ln_zps .AND. .NOT. ln_isfcav)                                   & 
     274               &             CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
     275               &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
     276         IF( ln_zps .AND.       ln_isfcav)                                   &  
     277               &             CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
     278               &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
     279               &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level 
    264280         ENDIF 
    265281         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection 

  • trunk/NEMOGCM/NEMO/OPA_SRC/timing.F90

    r3610 r5120  
    211211         WRITE(numtime,*) '                             NEMO team' 
    212212         WRITE(numtime,*) '                  Ocean General Circulation Model' 
    213          WRITE(numtime,*) '                        version 3.3  (2010) ' 
     213         WRITE(numtime,*) '                        version 3.6  (2015) ' 
    214214         WRITE(numtime,*) 
    215215         WRITE(numtime,*) '                        Timing Informations ' 

  • trunk/NEMOGCM/NEMO/SAS_SRC/nemogcm.F90

    r5118 r5120  
    250250         WRITE(numout,*) '                       NEMO team' 
    251251         WRITE(numout,*) '            Ocean General Circulation Model' 
    252          WRITE(numout,*) '                  version 3.4  (2011) ' 
     252         WRITE(numout,*) '                  version 3.6  (2015) ' 
    253253         WRITE(numout,*) '             StandAlone Surface version (SAS) ' 
    254254         WRITE(numout,*) 

  • trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90

    r4990 r5120  
    8282      IF( .NOT. Agrif_Root())   CALL Agrif_Update_Trc( kstp )   ! Update tracer at AGRIF zoom boundaries : children only 
    8383#endif 
    84          IF( ln_zps    )        CALL zps_hde( kstp, jptra, trn, pgtu=gtru, pgtv=gtrv, sgtu=gtrui, sgtv=gtrvi )  ! Partial steps: now horizontal gradient of passive 
     84 
     85         IF( ln_zps  .AND. .NOT. ln_isfcav)        & 
     86            &            CALL zps_hde    ( kstp, jptra, trn, gtru, gtrv )   ! Partial steps: now horizontal gradient of passive 
     87         IF( ln_zps .AND.        ln_isfcav)        & 
     88            &            CALL zps_hde_isf( kstp, jptra, trn, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! Partial steps: now horizontal gradient of passive 
    8589                                                                ! tracers at the bottom ocean level 
    8690         ! 

  • trunk/NEMOGCM/NEMO/TOP_SRC/trcini.F90

    r4990 r5120  
    143143  
    144144      tra(:,:,:,:) = 0._wp 
    145       IF( ln_zps .AND. .NOT. lk_c1d )   &              ! Partial steps: before horizontal gradient of passive 
    146         &    CALL zps_hde( nit000, jptra, trn, pgtu=gtru, pgtv=gtrv, sgtu=gtrui, sgtv=gtrvi )       ! tracers at the bottom ocean level 
     145      IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav )   &              ! Partial steps: before horizontal gradient of passive 
     146        &    CALL zps_hde    ( nit000, jptra, trn, gtru, gtrv  )  ! Partial steps: before horizontal gradient 
     147      IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav )   & 
     148        &    CALL zps_hde_isf( nit000, jptra, trn, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )       ! tracers at the bottom ocean level 
     149 
    147150 
    148151      ! 
Note: See TracChangeset for help on using the changeset viewer.