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 5111 for branches/2015 – NEMO

Changeset 5111 for branches/2015


Ignore:
Timestamp:
2015-03-02T16:30:57+01:00 (9 years ago)
Author:
mathiot
Message:

add some missing if ln_isfcav, test of option compatibility with ln_isfcav + small documentation

Location:
branches/2015/dev_r5094_UKMO_ISFCLEAN
Files:
1 added
29 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Biblio/Biblio.bib

    r4560 r5111  
    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 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Chapters/Chap_DOM.tex

    r4147 r5111  
    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  
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Chapters/Chap_DYN.tex

    r4759 r5111  
    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) 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Chapters/Chap_SBC.tex

    r4661 r5111  
    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 (parametrisation)  
     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. 
     64An other exemple of modification is the ice shelf melting/freezing (see \S\ref{SBC_isf}),  
     65which act as sources of fresh water due to ice shelf melting/freezing. 
    6266 
    6367 
     
    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 parametrisation 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} parametrisation of isf melting.  
     971The effective melting length (\np{sn\_Leff\_isf}) is read in a file. 
     972 
     973\item[\np{nn\_isf}~=~3] 
     974A simple parametrisation 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 in 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 in 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 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Chapters/Chap_ZDF.tex

    r4147 r5111  
    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 defined 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 same, only the bottom friction is described in details. 
     842 
    840843Both the surface momentum flux (wind stress) and the bottom momentum  
    841844flux (bottom friction) enter the equations as a condition on the vertical  
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/nambfr

    r4147 r5111  
    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/ 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/namdyn_hpg

    r4147 r5111  
    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) 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/namsbc

    r4230 r5111  
    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 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/DOC/TexFiles/Namelist/namzgr

    r3294 r5111  
    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/ 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/CONFIG/SHARED/namelist_ref

    r5099 r5111  
    8383   ln_zps      = .true.    !  z-coordinate - partial steps   (T/F) 
    8484   ln_sco      = .false.   !  s- or hybrid z-s-coordinate    (T/F) 
    85    ln_isfcav   = .false.   !  ice shelf cavity 
     85   ln_isfcav   = .false.   !  ice shelf cavity               (T/F) 
    8686/ 
    8787!----------------------------------------------------------------------- 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

    r4990 r5111  
    105105      END DO 
    106106      IF( .NOT.lk_vvl ) THEN 
    107          DO ji=1,jpi 
    108             DO jj=1,jpj 
    109                zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
    110             END DO 
    111          END DO 
     107         IF ( ln_isfcav ) THEN 
     108            DO ji=1,jpi 
     109               DO jj=1,jpj 
     110                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
     111               END DO 
     112            END DO 
     113         ELSE 
     114            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     115         END IF 
    112116      END IF 
    113117      !                                          
     
    127131      END DO 
    128132      IF( .NOT.lk_vvl ) THEN 
    129          DO ji=1,jpi 
    130             DO jj=1,jpj 
    131                zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
    132             END DO 
    133          END DO 
     133         IF ( ln_isfcav ) THEN 
     134            DO ji=1,jpi 
     135               DO jj=1,jpj 
     136                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj) 
     137               END DO 
     138            END DO 
     139         ELSE 
     140            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1) 
     141         END IF 
    134142      END IF 
    135143      !     
     
    157165      END DO 
    158166      IF( .NOT.lk_vvl ) THEN 
    159          DO ji=1,jpi 
    160             DO jj=1,jpj 
    161                ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem)  
    162                zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal)  
    163             END DO 
    164          END DO 
     167         IF ( ln_isfcav ) THEN 
     168            DO ji=1,jpi 
     169               DO jj=1,jpj 
     170                  ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem)  
     171                  zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal)  
     172               END DO 
     173            END DO 
     174         ELSE 
     175            ztemp = ztemp + zarea_ssh(:,:) * tsn(:,:,1,jp_tem)  
     176            zsal  = zsal  + zarea_ssh(:,:) * tsn(:,:,1,jp_sal)  
     177         END IF 
    165178      ENDIF 
    166179      IF( lk_mpp ) THEN   
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90

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

    r5098 r5111  
    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! 
     
    473475         ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 
    474476         IF( cp_cfg == "isomip" .AND. ln_isfcav ) 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  
     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         !  
    484485         ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN 
    485486         !  
     
    536537            CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy ) 
    537538            CALL iom_close( inum ) 
    538             !   
     539            !                                                 
    539540            risfdep(:,:)=0._wp          
    540541            misfdep(:,:)=1              
     
    963964      !!---------------------------------------------------------------------- 
    964965      !! 
    965       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
     966      INTEGER  ::   ji, jj, jk       ! dummy loop indices 
    966967      INTEGER  ::   ik, it           ! temporary integers 
    967       INTEGER  ::   id, jd, nprocd 
    968968      LOGICAL  ::   ll_print         ! Allow  control print for debugging 
    969969      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
    970970      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
    971       REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth 
     971      REAL(wp) ::   zmax             ! Maximum depth 
    972972      REAL(wp) ::   zdiff            ! temporary scalar 
    973973      REAL(wp) ::   zrefdep          ! temporary scalar 
     
    10461046                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
    10471047                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
    1048                   e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
    1049                      &                          / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
     1048                  e3t_0  (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
     1049                     &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
    10501050                  e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
    10511051                     &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
     
    12021202      
    12031203      ! Compute gdep3w_0 (vertical sum of e3w) 
    1204       IF ( ln_isfcav ) THEN 
     1204      IF ( ln_isfcav ) THEN ! if cavity 
    12051205         WHERE (misfdep == 0) misfdep = 1 
    12061206         DO jj = 1,jpj 
     
    13351335            misfdep(jpi,:) = misfdep(  2  ,:)  
    13361336         ENDIF 
    1337   
     1337 
    13381338         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
    13391339            mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west 
    13401340            mbathy(jpi,:) = mbathy(  2  ,:) 
    13411341         ENDIF 
    1342   
     1342 
    13431343         ! split last cell if possible (only where water column is 2 cell or less) 
    13441344         DO jk = jpkm1, 1, -1 
     
    13581358            END WHERE 
    13591359         END DO 
    1360   
     1360 
    13611361  
    13621362 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 
     
    16391639               IF( zmbathy(ji,jj) .LT. misfdep(ji  ,jj+1) ) ibtestjp1 = 0 
    16401640               ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 
    1641                IF( ibtest == 0 ) THEN 
     1641               IF( ibtest == 0 .AND. misfdep(ji,jj) .GE. 2) THEN 
    16421642                  mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 
    16431643               END IF 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r5104 r5111  
    140140            &            CALL zps_hde    ( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
    141141            &                                            rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
    142          IF( ln_zps .AND. ln_isfcav)                                       & 
     142         IF( ln_zps .AND.       ln_isfcav)                                 & 
    143143            &            CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
    144144            &                                            rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90

    r4990 r5111  
    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 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r5098 r5111  
    167167                           & the standard jacobian formulation hpg_sco or & 
    168168                           & the pressure jacobian formulation hpg_prj') 
     169 
     170      IF( ln_hpg_isf .AND. .NOT. ln_isfcav )   & 
     171           &   CALL ctl_stop( ' hpg_isf not available is ln_isfcav = false ' ) 
    169172      ! 
    170173      !                               ! Set nhpg from ln_hpg_... flags 
     
    187190      IF( ( .NOT. ln_hpg_isf ) .AND. ln_isfcav )   & 
    188191          &  CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' ) 
     192      !  
     193      ! initialisation of ice load 
     194      riceload(:,:)=0.0 
    189195      ! 
    190196   END SUBROUTINE dyn_hpg_init 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

    r5098 r5111  
    9090         DO jj = 2, jpjm1                 ! vertical momentum advection at w-point 
    9191            DO ji = fs_2, fs_jpim1        ! vector opt. 
    92                zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1)-un(ji,jj,jk) ) !* wumask(ji,jj,jk) 
    93                zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1)-vn(ji,jj,jk) ) !* wvmask(ji,jj,jk) 
     92               zwuw(ji,jj,jk) = ( zww(ji+1,jj  ) + zww(ji,jj) ) * ( un(ji,jj,jk-1)-un(ji,jj,jk) ) 
     93               zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1)-vn(ji,jj,jk) ) 
    9494            END DO   
    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 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

    r5104 r5111  
    157157         ! 
    158158         !==   Local vertical density gradient at T-point   == !   (evaluated from N^2) 
    159          ! surface initialisation  
    160          zdzr(:,:,1) = 0._wp  
    161159         ! interior value 
    162160         DO jk = 2, jpkm1 
     
    169167               &                 * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) ) 
    170168         END DO 
     169         ! surface initialisation  
     170         zdzr(:,:,1) = 0._wp  
    171171         IF ( ln_isfcav ) THEN 
    172          ! if isf need to overwrite the interior value at at the first ocean point 
     172            ! if isf need to overwrite the interior value at at the first ocean point 
    173173            DO jj = 1, jpjm1 
    174174               DO ji = 1, jpim1 
     
    185185         ! ===========================      | vslp = d/dj( prd ) / d/dz( prd ) 
    186186         ! 
    187          DO jj = 2, jpjm1 
    188             DO ji = fs_2, fs_jpim1   ! vector opt. 
    189                IF (miku(ji,jj) .GT. miku(ji+1,jj)) zhmlpu(ji,jj) = hmlpt(ji  ,jj) 
    190                IF (miku(ji,jj) .LT. miku(ji+1,jj)) zhmlpu(ji,jj) = hmlpt(ji+1,jj) 
    191                IF (miku(ji,jj) .EQ. miku(ji+1,jj)) zhmlpu(ji,jj) = MAX(hmlpt(ji  ,jj), hmlpt(ji+1,jj)) 
    192                IF (mikv(ji,jj) .GT. miku(ji,jj+1)) zhmlpv(ji,jj) = hmlpt(ji  ,jj) 
    193                IF (mikv(ji,jj) .LT. miku(ji,jj+1)) zhmlpv(ji,jj) = hmlpt(ji,jj+1) 
    194                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 
    195197            ENDDO 
    196          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 
    197206         DO jk = 2, jpkm1                            !* Slopes at u and v points 
    198207            DO jj = 2, jpjm1 
     
    208217                  zbv = MIN(  zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav )  ) 
    209218                  !                                      ! uslp and vslp output in zwz and zww, resp. 
    210                   zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )  
    211                   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) ) 
    212221                  ! thickness of water column between surface and level k at u/v point 
    213                   zdepu = 0.5_wp * (( fsdept(ji,jj,jk) + fsdept(ji+1,jj  ,jk) )                   & 
    214                              - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj  ) )  & 
    215                              - fse3u(ji,jj,miku(ji,jj))                                         ) 
    216                   zdepv = 0.5_wp * (( fsdept(ji,jj,jk) + fsdept(ji  ,jj+1,jk) )                   & 
    217                              - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) & 
    218                              - fse3v(ji,jj,mikv(ji,jj))                                         ) 
    219                   zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps )                                              & 
    220                      &                 + zfi  * uslpml(ji,jj)                                                     & 
    221                      &                        * zdepu / MAX( zhmlpu(ji,jj), 5._wp ) 
     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) 
    222229                  zwz(ji,jj,jk) = zwz(ji,jj,jk) * wumask(ji,jj,jk) 
    223                   zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps )                                              & 
    224                      &                 + zfj  * vslpml(ji,jj)                                                     & 
    225                      &                        * zdepv / MAX( zhmlpv(ji,jj), 5._wp )  
     230                  zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps )                                          & 
     231                     &                 + zfj  * vslpml(ji,jj) * zdepv / zhmlpv(ji,jj)  
    226232                  zww(ji,jj,jk) = zww(ji,jj,jk) * wvmask(ji,jj,jk) 
    227233                   
     
    276282                  uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk  ) ) * 0.5_wp   & 
    277283                     &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5_wp   & 
    278                      &                            *   umask(ji,jj,jk-1) !* umask(ji,jj,jk) * umask(ji,jj,jk+1) 
     284                     &                            *   umask(ji,jj,jk-1) 
    279285                  vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk  ) ) * 0.5_wp   & 
    280286                     &                            * ( vmask(ji  ,jj,jk) + vmask(ji  ,jj,jk+1) ) * 0.5_wp   & 
    281                      &                            *   vmask(ji,jj,jk-1) !* vmask(ji,jj,jk) * vmask(ji,jj,jk+1) 
     287                     &                            *   vmask(ji,jj,jk-1) 
    282288               END DO 
    283289            END DO 
     
    292298               DO ji = fs_2, fs_jpim1   ! vector opt. 
    293299                  !                                  !* Local vertical density gradient evaluated from N^2 
    294                   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) 
    295301                  !                                  !* Slopes at w point 
    296302                  !                                        ! i- & j-gradient of density at w-points 
     
    308314                  zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj )  ) 
    309315                  !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 
    310                   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 
    311317                  zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj), 10._wp ) 
    312318                  zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) & 
     
    746752            DO ji = 1, jpi 
    747753               ik = nmln(ji,jj) - 1 
    748                IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN   ;   omlmask(ji,jj,jk) = 1._wp 
    749                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 
    750758               ENDIF 
    751759            END DO 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90

    r5098 r5111  
    8989         ! 
    9090         IF( kn_fwb == 3 .AND. nn_sssr /= 2 )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 requires nn_sssr = 2, we stop ' ) 
    91          ! 
    92          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. 
    9396         ! 
    9497#if ! defined key_lim2 &&  ! defined key_lim3 && ! defined key_cice 
     
    107110            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain 
    108111            zcoef = z_fwf * rcp 
    109             emp(:,:) = emp(:,:) - z_fwf  
    110             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 
    111114         ENDIF 
    112115         ! 
     
    139142         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN         ! correct the freshwater fluxes 
    140143            zcoef = fwfold * rcp 
    141             emp(:,:) = emp(:,:) + fwfold 
    142             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 
    143146         ENDIF 
    144147         ! 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90

    r4990 r5111  
    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      !                                                   ! ---------------------------------------- ! 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r4990 r5111  
    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 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90

    r5098 r5111  
    147147         ! Surface value 
    148148         IF( lk_vvl ) THEN    
    149             DO jj = 1, jpj 
    150                DO ji = 1, jpi 
    151                   zwz(ji,jj, mikt(ji,jj) ) = 0.e0          ! volume variable 
    152                END DO 
    153             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 
    154158         ELSE                 
    155             DO jj = 1, jpj 
    156                DO ji = 1, jpi 
    157                   zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
    158                END DO 
    159             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 
    160168         ENDIF 
    161169 
     
    212220         END DO 
    213221         ! surface value 
    214          DO jj = 1, jpj 
    215             DO ji = 1, jpi 
    216                zwz(ji,jj,mikt(ji,jj)) = 0.e0 
    217             END DO 
    218          END DO 
     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 
    219231         CALL lbc_lnk( zwx, 'U', -1. )   ;   CALL lbc_lnk( zwy, 'V', -1. )         ! Lateral bondary conditions 
    220232         CALL lbc_lnk( zwz, 'W',  1. ) 
     
    360372 
    361373         ! upstream tracer flux in the k direction 
    362          ! Surface value 
    363          IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0._wp                        ! volume variable 
    364          ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface  
    365          ENDIF 
    366374         ! Interior value 
    367375         DO jk = 2, jpkm1 
     
    374382            END DO 
    375383         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 
    376406 
    377407         ! total advective trend 
     
    583613 
    584614      DO jk = 1, jpkm1 
     615         ikm1 = MAX(jk-1,1) 
     616         z2dtt = p2dt(jk) 
    585617         DO jj = 2, jpjm1 
    586618            DO ji = fs_2, fs_jpim1   ! vector opt. 
    587                ikm1 = MAX(jk-1,1) 
    588                z2dtt = p2dt(jk) 
    589                 
     619 
    590620               ! search maximum in neighbourhood 
    591621               zup = MAX(  zbup(ji  ,jj  ,jk  ),   & 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r4990 r5111  
    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 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90

    r5098 r5111  
    187187            END DO 
    188188         END DO 
    189             ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 
     189         ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 
    190190         zdk1t(:,:,1) = ( ptb(:,:,1,jn  ) - ptb(:,:,2,jn) ) * wmask(:,:,2) 
    191191         zdkt (:,:,1) = zdk1t(:,:,1) 
     
    200200         END IF 
    201201 
    202             ! 2. Horizontal fluxes 
    203             ! --------------------    
     202         ! 2. Horizontal fluxes 
     203         ! --------------------    
    204204         DO jk = 1, jpkm1 
    205205            DO jj = 1 , jpjm1 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r5098 r5111  
    8282      !!          di(rho) = rd~ - rd(i,j,k)   or   rd(i+1,j,k) - rd~ 
    8383      !! 
    84       !! ** Action  : compute for top and bottom interfaces 
     84      !! ** Action  : compute for top interfaces 
    8585      !!              - pgtu, pgtv: horizontal gradient of tracer at u- & v-points 
    8686      !!              - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points 
     
    9191      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts  
    9292      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) 
     93      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad of prd at u- & v-pts (bottom) 
    9494      ! 
    9595      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
    96       INTEGER  ::   iku, ikv, ikum1, ikvm1,ikup1, ikvp1   ! partial step level (ocean bottom level) at u- and v-points 
    97       REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv, zdzwu, zdzwv, zdzwuip1, zdzwvjp1  ! temporary scalars 
     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 
    9898      REAL(wp), DIMENSION(jpi,jpj)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos 
    9999      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::  zti, ztj             !  
     
    172172            END DO 
    173173         END DO 
    174           
     174 
    175175         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
    176176         ! step and store it in  zri, zrj for each  case 
     
    193193            END DO 
    194194         END DO 
    195          CALL lbc_lnk( pgru   , 'U', -1. )   ;   CALL lbc_lnk( pgrv  , 'V', -1. )   ! Lateral boundary conditions 
     195         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions 
    196196         ! 
    197197      END IF 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r4990 r5111  
    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 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90

    r4990 r5111  
    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      ! 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r5104 r5111  
    336336               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   & 
    337337                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   &  
    338                   &           / (  fse3uw_n(ji,jj,jk)         & 
    339                   &              * fse3uw_b(ji,jj,jk) ) 
     338                  &                            / (  fse3uw_n(ji,jj,jk)               & 
     339                  &                              *  fse3uw_b(ji,jj,jk) ) 
    340340               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   & 
    341341                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   & 
  • branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/TOP_SRC/trcini.F90

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