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 6012 – NEMO

Changeset 6012


Ignore:
Timestamp:
2015-12-07T16:11:45+01:00 (8 years ago)
Author:
mathiot
Message:

merge MetO branch with dev_r5151_UKMO_ISF

Location:
branches/2015/dev_MetOffice_merge_2015
Files:
37 edited
1 copied

Legend:

Unmodified
Added
Removed

  • branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Biblio/Biblio.bib

    r6005 r6012  
    289289  author = {A. Beckmann and H. Goosse}, 
    290290  title = {A parameterization of ice shelf-ocean interaction for climate models}, 
    291   journal = OM 
    292   year = {2003} 
    293   volume = {5} 
     291  journal = OM, 
     292  year = {2003}, 
     293  volume = {5}, 
    294294  pages = {157--170} 
    295295} 
    … …  
    13721372} 
    13731373 
     1374@ARTICLE{Holland1999,  
     1375  author = {D. Holland and A. Jenkins}, 
     1376  title = {Modeling Thermodynamic Ice-Ocean Interactions at the Base of an Ice Shelf}, 
     1377  journal = JPO, 
     1378  year = {1999}, 
     1379  volume = {29},   
     1380  pages = {1787--1800}, 
     1381} 
     1382 
    13741383@ARTICLE{HollowayOM86, 
    13751384  author = {Greg Holloway}, 
    … …  
    15161525  journal = GRL, 
    15171526  pages = {811--814} 
     1527} 
     1528 
     1529@ARTICLE{Jenkins1991, 
     1530  author = {A. Jenkins}, 
     1531  title = {A one-dimensional model of ice shelf-ocean interaction}, 
     1532  journal = JGR, 
     1533  year = {1991}, 
     1534  volume = {96},  number = {C11}, 
     1535  pages = {2298--2312} 
     1536} 
     1537 
     1538@ARTICLE{Jenkins2010, 
     1539  author = {A. Jenkins}, 
     1540  title = {observation and parameterization of ablation at the base of Ronne Ice Shelf, Antarctica}, 
     1541  journal = JPO, 
     1542  year = {2010}, 
     1543  volume = {40},  number = {10}, 
     1544  pages = {2298--2312} 
    15181545} 
    15191546 
    … …  
    19231950  volume = {51}, 
    19241951  pages = {737--769} 
     1952} 
     1953 
     1954@ARTICLE{Losch2008, 
     1955  author = {M. Losch}, 
     1956  title = {Modeling ice shelf cavities in a z coordinate ocean general circulation model}, 
     1957  journal = JGR, 
     1958  year = {2008}, 
     1959  volume = {113},  number = {C13}, 
    19251960} 
    19261961 

  • branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Chapters/Chap_SBC.tex

    r6006 r6012  
    4949(\np{nn\_ice}~=~0,1, 2 or 3); the addition of river runoffs as surface freshwater  
    5050fluxes 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);  
     51 or as surface flux at the land-ice ocean interface (\np{ln\_isf}=~true);  
    5352the addition of a freshwater flux adjustment in order to avoid a mean sea-level drift (\np{nn\_fwb}~=~0,~1~or~2); the  
    5453transformation of the solar radiation (if provided as daily mean) into a diurnal  
    … …  
    958957\namdisplay{namsbc_isf} 
    959958%-------------------------------------------------------------------------------------------------------- 
    960 Namelist variable in \ngn{namsbc}, \np{nn\_isf},  control the kind of ice shelf representation used.  
     959Namelist variable in \ngn{namsbc}, \np{nn\_isf}, control the kind of ice shelf representation used.  
    961960\begin{description} 
    962961\item[\np{nn\_isf}~=~1] 
    963 The ice shelf cavity is represented. The fwf and heat flux are computed.  
     962The ice shelf cavity is represented. The fwf and heat flux are computed. 2 bulk formulations are available: the ISOMIP one (\np{nn\_isfblk = 1}) described in (\np{nn\_isfblk = 2}), the 3 equation formulation described in \citet{Jenkins1991}. In addition to this,  
     9633 different way to compute the exchange coefficient are available. $\gamma\_{T/S}$ is constant (\np{nn\_gammablk = 0}), $\gamma\_{T/S}$ is velocity dependant \citep{Jenkins2010} (\np{nn\_gammablk = 1}) and $\gamma\_{T/S}$ is velocity dependant and stratification dependent \citep{Holland1999} (\np{nn\_gammablk = 2}). For each of them, the thermal/salt exchange coefficient (\np{rn\_gammat0} and \np{rn\_gammas0}) have to be specified (the default values are for the ISOMIP case).  
    964964Full description, sensitivity and validation in preparation.  
    965965 
    … …  
    969969(\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).  
    970970Furthermore the fwf is computed using the \citet{Beckmann2003} parameterisation of isf melting.  
    971 The effective melting length (\np{sn\_Leff\_isf}) is read from a file. 
     971The effective melting length (\np{sn\_Leff\_isf}) is read from a file and the exchange coefficients are set as (\np{rn\_gammat0}) and (\np{rn\_gammas0}). 
    972972 
    973973\item[\np{nn\_isf}~=~3] 
    … …  
    987987\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. 
    988988 
    989  This can be usefull if the water masses on the shelf are not realistic or the resolution (horizontal/vertical) are too  
     989This can be usefull if the water masses on the shelf are not realistic or the resolution (horizontal/vertical) are too  
    990990coarse to have realistic melting or for sensitivity studies where you want to control your input.  
    991991Full description, sensitivity and validation in preparation.  
    992992 
    993 There 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}).  
     993\np{rn\_hisf\_tbl} is the top boundary layer (tbl) thickness used by the Losch parametrisation \citep{Losch2008} to compute the melt. if 0, temperature/salt/velocity in the top cell is used to compute the melt. 
     994Otherwise, NEMO used the mean value into the tbl.  
    998995 
    999996\section{ Ice sheet coupling} 

  • branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namsbc

    r5120 r6012  
    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) 
     21   ln_isf      = .false.   !  ice shelf melting/freezing                (T => fill namsbc_isf) 
    2622   ln_ssr      = .true.    !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr) 
    2723   nn_fwb      = 3         !  FreshWater Budget: =0 unchecked 

  • branches/2015/dev_MetOffice_merge_2015/DOC/TexFiles/Namelist/namsbc_isf

    r5120 r6012  
    22&namsbc_isf    !  Top boundary layer (ISF) 
    33!----------------------------------------------------------------------- 
     4   nn_isf      = -99       !  ice shelf melting/freezing                (/=0 => fill namsbc_isf) 
     5                           !  1 = ISF explicit (cavity open) 
     6                           !  2 = bg03 parametrisation (cavity closed)   
     7                           !  3 = ISF  specified in depth along the calving front (cavity closed) 
     8                           !  4 = ISF melting specified (cavity open) 
     9                           !  option 1 and 4 need ln_isfcav = .true. (domzgr) 
    410!              ! file name ! frequency (hours) ! variable ! time interpol. !  clim   ! 'yearly'/ ! weights  ! rotation ! 
    511!              !           !  (if <0  months)  !   name   !    (logical)   !  (T/F)  ! 'monthly' ! filename ! pairing  ! 
    612! nn_isf == 4 
    7    sn_qisf      = 'rnfisf' ,         -12      ,'sohflisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   '' 
    813   sn_fwfisf    = 'rnfisf' ,         -12      ,'sowflisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   '' 
    914! nn_isf == 3 
    1015   sn_rnfisf    = 'runoffs' ,         -12      ,'sofwfisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   '' 
    11 ! nn_isf == 2 and 3 
     16! nn_isf == 2 or 3 
    1217   sn_depmax_isf = 'runoffs' ,       -12        ,'sozisfmax' ,   .false.  , .true.  , 'yearly'  ,  ''      ,   '' 
    1318   sn_depmin_isf = 'runoffs' ,       -12        ,'sozisfmin' ,   .false.  , .true.  , 'yearly'  ,  ''      ,   '' 
    1419! nn_isf == 2 
    1520   sn_Leff_isf = 'rnfisf' ,       0          ,'Leff'         ,   .false.  , .true.  , 'yearly'  ,  ''      ,   '' 
    16 ! for all case 
    17    ln_divisf   = .true.  ! apply isf melting as a mass flux or in the salinity trend. (maybe I should remove this option as for runoff?) 
    1821! only for nn_isf = 1 or 2 
    1922   rn_gammat0  = 1.0e-4   ! gammat coefficient used in blk formula 
    2023   rn_gammas0  = 1.0e-4   ! gammas coefficient used in blk formula 
     24! only for nn_isf = 1 or 4 
     25   nn_isfblk   =  1       ! 1 ISOMIP ; 2 conservative (3 equation formulation, Jenkins et al. 1991 ??) 
    2126! only for nn_isf = 1 
    22    nn_isfblk   =  1       ! 1 ISOMIP ; 2 conservative (3 equation formulation, Jenkins et al. 1991 ??) 
    23    rn_hisf_tbl =  30.      ! thickness of the top boundary layer           (Losh et al. 2008) 
     27   rn_hisf_tbl =  30.     ! thickness of the top boundary layer           (Losh et al. 2008) 
    2428                          ! 0 => thickness of the tbl = thickness of the first wet cell 
    25    ln_conserve = .true.   ! conservative case (take into account meltwater advection) 
    2629   nn_gammablk = 1        ! 0 = cst Gammat (= gammat/s) 
    2730                          ! 1 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010) 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/SHARED/namelist_ref

    r6008 r6012  
    245245   ln_dm2dc    = .false.   !  daily mean to diurnal cycle on short wave 
    246246   ln_rnf      = .true.    !  runoffs                                   (T   => fill namsbc_rnf) 
    247    nn_isf      = 0         !  ice shelf melting/freezing                (/=0 => fill namsbc_isf) 
    248                            !  0 =no isf                  1 = presence of ISF 
    249                            !  2 = bg03 parametrisation   3 = rnf file for isf 
    250                            !  4 = ISF fwf specified 
    251                            !  option 1 and 4 need ln_isfcav = .true. (domzgr) 
     247   ln_isf      = .false.   !  ice shelf                                 (T   => fill namsbc_isf) 
    252248   ln_ssr      = .true.    !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr) 
    253249   nn_fwb      = 2         !  FreshWater Budget: =0 unchecked 
    … …  
    431427&namsbc_isf    !  Top boundary layer (ISF) 
    432428!----------------------------------------------------------------------- 
    433 !              ! file name ! frequency (hours) ! variable ! time interpol. !  clim   ! 'yearly'/ ! weights  ! rotation ! 
    434 !              !           !  (if <0  months)  !   name   !    (logical)   !  (T/F)  ! 'monthly' ! filename ! pairing  ! 
     429!              ! file name ! frequency (hours) ! variable ! time interpol. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
     430!              !           !  (if <0  months)  !   name   !    (logical)   !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      ! 
    435431! nn_isf == 4 
    436    sn_qisf      = 'rnfisf' ,         -12      ,'sohflisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   '' 
    437    sn_fwfisf    = 'rnfisf' ,         -12      ,'sowflisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   '' 
     432   sn_fwfisf   = 'rnfisf'  ,         -12       ,'sowflisf',     .false.    , .true.  , 'yearly'  ,  ''      ,   ''     , '' 
    438433! nn_isf == 3 
    439    sn_rnfisf    = 'runoffs' ,         -12      ,'sofwfisf',    .false.      , .true.  , 'yearly'  ,  ''      ,  '' 
     434   sn_rnfisf   = 'rnfisf'  ,         -12       ,'sofwfisf',     .false.    , .true.  , 'yearly'  ,  ''      ,   ''     , '' 
    440435! nn_isf == 2 and 3 
    441    sn_depmax_isf = 'runoffs' ,       -12        ,'sozisfmax' ,   .false.  , .true.  , 'yearly'  ,  ''      ,  '' 
    442    sn_depmin_isf = 'runoffs' ,       -12        ,'sozisfmin' ,   .false.  , .true.  , 'yearly'  ,  ''      ,  '' 
     436   sn_depmax_isf = 'rnfisf' ,        -12       ,'sozisfmax' ,   .false.    , .true.  , 'yearly'  ,  ''      ,   ''     , '' 
     437   sn_depmin_isf = 'rnfisf' ,        -12       ,'sozisfmin' ,   .false.    , .true.  , 'yearly'  ,  ''      ,   ''     , '' 
    443438! nn_isf == 2 
    444    sn_Leff_isf = 'rnfisf' ,       0          ,'Leff'         ,   .false.  , .true.  , 'yearly'  ,  ''      ,  '' 
     439   sn_Leff_isf = 'rnfisf'  ,         -12       ,'Leff'    ,     .false.    , .true.  , 'yearly'  ,  ''      ,   ''     , '' 
    445440! for all case 
    446    ln_divisf   = .true.  ! apply isf melting as a mass flux or in the salinity trend. (maybe I should remove this option as for runoff?) 
     441   nn_isf      = 1         !  ice shelf melting/freezing 
     442                           !  1 = presence of ISF    2 = bg03 parametrisation  
     443                           !  3 = rnf file for isf   4 = ISF fwf specified 
     444                           !  option 1 and 4 need ln_isfcav = .true. (domzgr) 
    447445! only for nn_isf = 1 or 2 
    448446   rn_gammat0  = 1.0e-4   ! gammat coefficient used in blk formula 
    449447   rn_gammas0  = 1.0e-4   ! gammas coefficient used in blk formula 
     448! only for nn_isf = 1 or 4 
     449   rn_hisf_tbl =  30.      ! thickness of the top boundary layer    (Losh et al. 2008) 
     450                          ! 0 => thickness of the tbl = thickness of the first wet cell 
    450451! only for nn_isf = 1 
    451    nn_isfblk   =  1       ! 1 ISOMIP ; 2 conservative (3 equation formulation, Jenkins et al. 1991 ??) 
    452    rn_hisf_tbl =  30.      ! thickness of the top boundary layer           (Losh et al. 2008) 
    453                           ! 0 => thickness of the tbl = thickness of the first wet cell 
    454    ln_conserve = .true.   ! conservative case (take into account meltwater advection) 
     452   nn_isfblk   = 1        ! 1 ISOMIP  like: 2 equations formulation (Hunter et al., 2006) 
     453                          ! 2 ISOMIP+ like: 3 equations formulation (Asay-Davis et al., 2015) 
    455454   nn_gammablk = 1        ! 0 = cst Gammat (= gammat/s) 
    456455                          ! 1 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010) 
    457                           !     if you want to keep the cd as in global config, adjust rn_gammat0 to compensate 
    458                           ! 2 = velocity and stability dependent Gamma    Holland et al. 1999 
     456                          ! 2 = velocity and stability dependent Gamma    (Holland et al. 1999) 
    459457/ 
    460458!----------------------------------------------------------------------- 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/cfg.txt

    r6006 r6012  
    77AMM12 OPA_SRC 
    88ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
    9 GYRE OPA_SRC 
    109ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC 
    1110ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 
    1211ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC 
     12GYRE OPA_SRC 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/CONFIG/uspcfg.txt

    r5333 r6012  
    11ORCA1_CICE # ORCA2_LIM # OPA_SRC TOP_SRC  # http://gws-access.ceda.ac.uk/public/nemo/uspconfigs/ORCA1_CICE/v3.6.0/ORCA1_CICE_ctl.txt 
    22ISOMIP     # GYRE      # OPA_SRC          # http://gws-access.ceda.ac.uk/public/nemo/uspconfigs/ISOMIP/v3.6.0/ISOMIP_ctl.txt 
     3ISOMIP_r5151_UKMO_ISF     # GYRE      # OPA_SRC          # http://gws-access.ceda.ac.uk/public/nemo/uspconfigs/ISOMIP_r5151_UKMO_ISF/v3.6.0/ISOMIP_ctl.txt 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90

    r5930 r6012  
    489489         CALL bn2    ( pts, rab_n, rn2  ) ! now    Brunt-Vaisala 
    490490 
    491          ! Partial steps: before Horizontal DErivative 
    492          IF( ln_zps  .AND. .NOT. ln_isfcav)                            & 
    493             &            CALL zps_hde    ( kt, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
    494             &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
    495          IF( ln_zps .AND.        ln_isfcav)                            & 
    496             &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv,  &    ! Partial steps for top cell (ISF) 
    497             &                                        rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
    498             &                                 gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level 
     491      ! Partial steps: before Horizontal DErivative 
     492      IF( ln_zps  .AND. .NOT. ln_isfcav)                            & 
     493         &            CALL zps_hde    ( kt, jpts, pts, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient 
     494         &                                        rhd, gru , grv    )  ! of t, s, rd at the last ocean level 
     495      IF( ln_zps .AND.        ln_isfcav)                            & 
     496         &            CALL zps_hde_isf( kt, jpts, pts, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
     497         &                                        rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level 
    499498 
    500499         rn2b(:,:,:) = rn2(:,:,:)         ! need for zdfmxl 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r5997 r6012  
    636636!!gm 
    637637 
    638             IF( ln_zps .AND. .NOT. lk_c1d ) THEN      ! Partial steps: before horizontal gradient 
    639                IF(ln_isfcav) THEN                        ! ocean cavities: top and bottom cells (ISF) 
    640                   CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi,     & 
    641                      &                            rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
    642                      &                     grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) 
    643                ELSE                                      ! no ocean cavities: bottom cells 
    644                   CALL zps_hde    ( kt, jpts, tsb, gtsu, gtsv,        &  !  
    645                      &                        rhd, gru , grv          )  ! of t, s, rd at the last ocean level 
    646                ENDIF 
    647             ENDIF 
    648             ! 
     638            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)      & 
     639               &  CALL zps_hde    ( kt, jpts, tsb, gtsu, gtsv,        &  ! Partial steps: before horizontal gradient 
     640               &                              rhd, gru , grv          )  ! of t, s, rd at the last ocean level 
     641            IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)      & 
     642               &  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, gtui, gtvi,    &    ! Partial steps for top cell (ISF) 
     643               &                                  rhd, gru , grv , grui, grvi     )    ! of t, s, rd at the last ocean level 
     644 
     645#if defined key_zdfkpp 
     646            CALL eos( tsn, rhd, fsdept_n(:,:,:) )                      ! Compute rhd 
     647!!gm fabien            CALL eos( tsn, rhd )                      ! Compute rhd 
     648#endif 
     649 
    649650            DEALLOCATE( t_bkginc ) 
    650651            DEALLOCATE( s_bkginc ) 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90

    r5930 r6012  
    9393      ! ----------------------------------------------------------------------- 
    9494!!gm replace these lines : 
    95       z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:)+fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rau0 
     95      z_cflxemp = SUM ( ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rau0 
    9696      IF( lk_mpp )   CALL mpp_sum( z_cflxemp )     ! sum over the global domain 
    9797!!gm   by : 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90

    r6006 r6012  
    194194                  DO ji = 1,jpi 
    195195                     ! Elevation 
    196                      ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*tmask_i(ji,jj)         
     196                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj)         
    197197                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj) 
    198198                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj) 
    … …  
    358358               X1=ana_amp(ji,jj,jh,1) 
    359359               X2=-ana_amp(ji,jj,jh,2) 
    360                out_v(ji,jj,       jh)=X1 * vmask_i(ji,jj) 
    361                out_v(ji,jj,nb_ana+jh)=X2 * vmask_i(ji,jj) 
     360               out_v(ji,jj,       jh)=X1 * ssvmask(ji,jj) 
     361               out_v(ji,jj,nb_ana+jh)=X2 * ssvmask(ji,jj) 
    362362            END DO 
    363363         END DO 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90

    r6006 r6012  
    101101      IF( ln_rnf_sal)   z_frc_trd_s = z_frc_trd_s + glob_sum( rnf_tsc(:,:,jp_sal) * surf(:,:) ) 
    102102      ! Add ice shelf heat & salt input 
    103       IF( nn_isf .GE. 1 )  THEN 
    104           z_frc_trd_t = z_frc_trd_t + glob_sum( risf_tsc(:,:,jp_tem) * surf(:,:) ) 
    105           z_frc_trd_s = z_frc_trd_s + glob_sum( risf_tsc(:,:,jp_sal) * surf(:,:) ) 
    106       ENDIF 
    107  
     103      IF( ln_isf    ) z_frc_trd_t = z_frc_trd_t + glob_sum( risf_tsc(:,:,jp_tem) * surf(:,:) ) 
    108104      ! Add penetrative solar radiation 
    109105      IF( ln_traqsr )   z_frc_trd_t = z_frc_trd_t + r1_rau0_rcp * glob_sum( qsr     (:,:) * surf(:,:) ) 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r6006 r6012  
    256256   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i            !: interior domain T-point mask 
    257257   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_h            !: internal domain T-point mask (Figure 8.5 NEMO book) 
    258    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask                     !: land/ocean mask of barotropic stream function 
     258   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask              !: land/ocean mask of barotropic stream function 
    259259 
    260260   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   misfdep                 !: top first ocean level                (ISF) 
    … …  
    390390 
    391391      ALLOCATE( mbathy(jpi,jpj) , bathy  (jpi,jpj) ,                                       & 
    392          &     tmask_i(jpi,jpj) , tmask_h(jpi, jpj),                                       &  
     392         &     tmask_i(jpi,jpj) , tmask_h(jpi,jpj) ,                                       &  
    393393         &     ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , & 
    394394         &     bmask(jpi,jpj)   ,                                                          & 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90

    r5836 r6012  
    506506      CALL iom_close( inum ) 
    507507       
    508 !!gm   THIS is TO BE REMOVED !!!!!!! 
    509  
    510 ! need to be define for the extended grid south of -80S 
    511 ! some point are undefined but you need to have e1 and e2 .NE. 0 
    512       WHERE (e1t==0.0_wp) 
    513          e1t=1.0e2 
    514       END WHERE 
    515       WHERE (e1v==0.0_wp) 
    516          e1v=1.0e2 
    517       END WHERE 
    518       WHERE (e1u==0.0_wp) 
    519          e1u=1.0e2 
    520       END WHERE 
    521       WHERE (e1f==0.0_wp) 
    522          e1f=1.0e2 
    523       END WHERE 
    524       WHERE (e2t==0.0_wp) 
    525          e2t=1.0e2 
    526       END WHERE 
    527       WHERE (e2v==0.0_wp) 
    528          e2v=1.0e2 
    529       END WHERE 
    530       WHERE (e2u==0.0_wp) 
    531          e2u=1.0e2 
    532       END WHERE 
    533       WHERE (e2f==0.0_wp) 
    534          e2f=1.0e2 
    535       END WHERE 
    536 !!gm end 
    537         
    538508    END SUBROUTINE hgr_read 
    539509     

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r6006 r6012  
    975975      ! 
    976976      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 
    977       IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 
     977      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 
    978978      ! 
    979979      IF(lwp) THEN                   ! Print the choice 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r6006 r6012  
    12581258      END WHERE 
    12591259 
     1260      ! remove very shallow ice shelf (less than ~ 10m if 75L) 
     1261      WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1) 
     1262         misfdep = 0; risfdep = 0.0_wp; 
     1263         mbathy  = 0; bathy   = 0.0_wp; 
     1264      END WHERE 
     1265      WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp) 
     1266         misfdep = 0; risfdep = 0.0_wp; 
     1267         mbathy  = 0; bathy   = 0.0_wp; 
     1268      END WHERE 
     1269  
    12601270! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation 
    12611271      icompt = 0  
    … …  
    13231333                     zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj)  ) & 
    13241334                         &            - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 
    1325   
    13261335                     IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) <  misfdep(ji,jj)) THEN 
    13271336                        IF (zbathydiff <= zrisfdepdiff) THEN 
    … …  
    13691378         END DO 
    13701379  
    1371  ! point V mbathy(ji,jj) EQ misfdep(ji,jj+1)  
     1380 ! point V mbathy(ji,jj) == misfdep(ji,jj+1)  
    13721381         DO jj = 1, jpjm1 
    13731382            DO ji = 1, jpim1 
    … …  
    14051414            mbathy(:,:)  = INT( zbathy(:,:) ) 
    14061415         ENDIF 
    1407  ! point V misdep(ji,jj) EQ mbathy(ji,jj+1)  
     1416 ! point V misdep(ji,jj) == mbathy(ji,jj+1)  
    14081417         DO jj = 1, jpjm1 
    14091418            DO ji = 1, jpim1 
    … …  
    14431452         ENDIF  
    14441453  
    1445  ! point U mbathy(ji,jj) EQ misfdep(ji,jj+1)  
     1454 ! point U mbathy(ji,jj) == misfdep(ji,jj+1)  
    14461455         DO jj = 1, jpjm1 
    14471456            DO ji = 1, jpim1 
    … …  
    14801489         ENDIF  
    14811490  
    1482  ! point U misfdep(ji,jj) EQ bathy(ji,jj+1)  
     1491 ! point U misfdep(ji,jj) == bathy(ji,jj+1)  
    14831492         DO jj = 1, jpjm1 
    14841493            DO ji = 1, jpim1 
    … …  
    17491758! end check compatibility ice shelf/bathy 
    17501759      ! remove very shallow ice shelf (less than ~ 10m if 75L) 
     1760      WHERE (risfdep(:,:) <= 10._wp) 
     1761         misfdep = 1; risfdep = 0.0_wp; 
     1762      END WHERE 
     1763 
    17511764      IF( icompt == 0 ) THEN  
    17521765         IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry'  
    … …  
    18321845               ENDIF  
    18331846            !       ... on ik / ik-1  
    1834                e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik)  
     1847               e3w_0  (ji,jj,ik  ) = e3t_0  (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
    18351848               e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
    18361849! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90

    r6006 r6012  
    2020   USE oce             ! ocean dynamics and tracers 
    2121   USE dom_oce         ! ocean space and time domain 
    22    USE sbc_oce, ONLY : ln_rnf, nn_isf ! surface boundary condition: ocean 
     22   USE sbc_oce, ONLY : ln_rnf, ln_isf ! surface boundary condition: ocean 
    2323   USE sbcrnf          ! river runoff  
    2424   USE sbcisf          ! ice shelf 
    … …  
    9191      END DO 
    9292      ! 
    93       IF( ln_rnf                     )   CALL sbc_rnf_div( hdivn )      !==  runoffs    ==!   (update hdivn field) 
     93      IF( ln_rnf )   CALL sbc_rnf_div( hdivn )      !==  runoffs    ==!   (update hdivn field) 
    9494      ! 
    95       IF( ln_divisf .AND. nn_isf > 0 )   CALL sbc_isf_div( hdivn )      !==  ice shelf  ==!   (update hdivn field) 
     95      IF( ln_isf )   CALL sbc_isf_div( hdivn )      !==  ice shelf  ==!   (update hdivn field) 
    9696      ! 
    97       IF( ln_iscpl  .AND. ln_hsb     )   CALL iscpl_div( hdivn )   !==  ice sheet  ==!   (update hdivn field) 
     97      IF( ln_iscpl .AND. ln_hsb ) CALL iscpl_div( hdivn )  !==  ice sheet  ==!   (update hdivn field) 
    9898      ! 
    99       CALL lbc_lnk( hdivn, 'T', 1. )                       !==  lateral boundary cond.  ==!   (no sign change) 
     99      CALL lbc_lnk( hdivn, 'T', 1. )                !==  lateral boundary cond.  ==!   (no sign change) 
    100100      ! 
    101101      IF( nn_timing == 1 )  CALL timing_stop('div_hor') 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r5930 r6012  
    4545   USE wrk_nemo        ! Memory Allocation 
    4646   USE timing          ! Timing 
     47   USE iom 
    4748 
    4849   IMPLICIT NONE 
    … …  
    130131      INTEGER ::   ioptio = 0      ! temporary integer 
    131132      INTEGER ::   ios             ! Local integer output status for namelist read 
     133      !! 
     134      INTEGER  ::   ji, jj, jk, ikt    ! dummy loop indices      ISF 
     135      REAL(wp) ::   znad 
     136      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztstop, zrhd ! hypothesys on isf density 
     137      REAL(wp), POINTER, DIMENSION(:,:)     ::  zrhdtop_isf  ! density at bottom of ISF 
     138      REAL(wp), POINTER, DIMENSION(:,:)     ::  ziceload     ! density at bottom of ISF 
    132139      !! 
    133140      NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     & 
    … …  
    190197      IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 
    191198      !  
    192       ! initialisation of ice load 
    193       riceload(:,:)=0.0 
     199      ! initialisation of ice shelf load 
     200      IF ( .NOT. ln_isfcav ) riceload(:,:)=0.0 
     201      IF (       ln_isfcav ) THEN 
     202         CALL wrk_alloc( jpi,jpj, 2,  ztstop)  
     203         CALL wrk_alloc( jpi,jpj,jpk, zrhd  ) 
     204         CALL wrk_alloc( jpi,jpj,     zrhdtop_isf, ziceload)  
     205         ! 
     206         IF(lwp) WRITE(numout,*) 
     207         IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 
     208         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'   
     209 
     210         ! To use density and not density anomaly 
     211         znad=1._wp 
     212          
     213         ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 
     214         ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 
     215 
     216         ! compute density of the water displaced by the ice shelf  
     217         DO jk = 1, jpk 
     218            CALL eos(ztstop(:,:,:),fsdept_n(:,:,jk),zrhd(:,:,jk)) 
     219         END DO 
     220       
     221         ! compute rhd at the ice/oce interface (ice shelf side) 
     222         CALL eos(ztstop,risfdep,zrhdtop_isf) 
     223 
     224         ! Surface value + ice shelf gradient 
     225         ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 
     226         ! divided by 2 later 
     227         ziceload = 0._wp 
     228         DO jj = 1, jpj 
     229            DO ji = 1, jpi 
     230               ikt=mikt(ji,jj) 
     231               ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 
     232               DO jk=2,ikt-1 
     233                  ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 
     234                     &                              * (1._wp - tmask(ji,jj,jk)) 
     235               END DO 
     236               IF (ikt  >=  2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 
     237                                  &                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 
     238            END DO 
     239         END DO 
     240         riceload(:,:)=ziceload(:,:)  ! need to be saved for diaar5 
     241 
     242         CALL wrk_dealloc( jpi,jpj, 2,  ztstop)  
     243         CALL wrk_dealloc( jpi,jpj,jpk, zrhd  ) 
     244         CALL wrk_dealloc( jpi,jpj,     zrhdtop_isf, ziceload)  
     245      END IF 
    194246      ! 
    195247   END SUBROUTINE dyn_hpg_init 
    … …  
    447499   SUBROUTINE hpg_isf( kt ) 
    448500      !!--------------------------------------------------------------------- 
    449       !!                  ***  ROUTINE hpg_sco  *** 
     501      !!                  ***  ROUTINE hpg_isf  *** 
    450502      !! 
    451503      !! ** Method  :   s-coordinate case. Jacobian scheme. 
    … …  
    466518      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    467519      !! 
    468       INTEGER  ::   ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j                 ! dummy loop indices 
    469       REAL(wp) ::   zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1             ! temporary scalars 
    470       REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zhpi, zhpj, zrhd 
     520      INTEGER  ::   ji, jj, jk, ikt, iktp1i, iktp1j   ! dummy loop indices 
     521      REAL(wp) ::   zcoef0, zuap, zvap, znad          ! temporary scalars 
     522      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zhpi, zhpj 
    471523      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztstop 
    472       REAL(wp), POINTER, DIMENSION(:,:)     ::  ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj 
    473       !!---------------------------------------------------------------------- 
    474       ! 
    475       CALL wrk_alloc( jpi,jpj, 2, ztstop)  
    476       CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 
    477       CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj)  
    478       ! 
    479      IF( kt == nit000 ) THEN 
    480          IF(lwp) WRITE(numout,*) 
    481          IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 
    482          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used' 
    483       ENDIF 
    484  
     524      REAL(wp), POINTER, DIMENSION(:,:)     ::  zrhdtop_oce 
     525      !!---------------------------------------------------------------------- 
     526      ! 
     527      CALL wrk_alloc( jpi,jpj,  2, ztstop)  
     528      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj) 
     529      CALL wrk_alloc( jpi,jpj,     zrhdtop_oce ) 
     530      ! 
    485531      ! Local constant initialization 
    486532      zcoef0 = - grav * 0.5_wp 
     533    
    487534      ! To use density and not density anomaly 
    488 !      IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume 
    489 !      ELSE                 ;     znad = 0._wp         ! Fixed volume 
    490 !      ENDIF 
    491535      znad=1._wp 
     536 
    492537      ! iniitialised to 0. zhpi zhpi  
    493538      zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 
    494539 
    495       ! Partial steps: top & bottom before horizontal gradient 
    496 !jc      CALL zps_hde_isf( kt, jpts, tsn, gtsu, gtsv, gtui, gtvi,   &  
    497 !jc               &       rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
    498 !jc               &      grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  ) 
    499  
    500 !==================================================================================      
    501 !=====Compute iceload and contribution of the half first wet layer =================  
    502 !=================================================================================== 
    503  
    504       ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 
    505       ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 
    506  
    507       ! compute density of the water displaced by the ice shelf  
    508       zrhd = rhd ! save rhd 
    509       DO jk = 1, jpk 
    510            zdept(:,:)=gdept_1d(jk) 
    511            CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk)) 
    512       END DO 
    513       WHERE ( tmask(:,:,:) == 1._wp) 
    514         rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd 
    515       END WHERE 
    516        
    517       ! compute rhd at the ice/oce interface (ice shelf side) 
    518       CALL eos(ztstop,risfdep,zrhdtop_isf) 
    519  
    520540      ! compute rhd at the ice/oce interface (ocean side) 
     541      ! usefull to reduce residual current in the test case ISOMIP with no melting 
    521542      DO ji=1,jpi 
    522543        DO jj=1,jpj 
    … …  
    527548      END DO 
    528549      CALL eos(ztstop,risfdep,zrhdtop_oce) 
    529       ! 
    530       ! Surface value + ice shelf gradient 
    531       ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 
    532       ziceload = 0._wp 
    533       DO jj = 1, jpj 
    534          DO ji = 1, jpi   ! vector opt. 
    535             ikt=mikt(ji,jj) 
    536             ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 
    537             DO jk=2,ikt-1 
    538                ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 
    539                   &                              * (1._wp - tmask(ji,jj,jk)) 
    540             END DO 
    541             IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) & 
    542                                &                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 
    543          END DO 
    544       END DO 
    545       riceload(:,:) = 0.0_wp ; riceload(:,:)=ziceload(:,:)  ! need to be saved for diaar5 
    546       ! compute zp from z=0 to first T wet point (correction due to zps not yet applied) 
     550 
     551!==================================================================================      
     552!===== Compute surface value =====================================================  
     553!================================================================================== 
    547554      DO jj = 2, jpjm1 
    548555         DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    550557            ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 
    551558            ! we assume ISF is in isostatic equilibrium 
    552             zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj  ,iktp1i)                                    & 
    553                &                                  * ( 2._wp * znad + rhd(ji+1,jj  ,iktp1i) + zrhdtop_oce(ji+1,jj  ) )   & 
    554                &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    & 
    555                &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   & 
    556                &                                  + ( ziceload(ji+1,jj) - ziceload(ji,jj))                              )  
    557             zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji  ,jj+1,iktp1j)                                    & 
    558                &                                  * ( 2._wp * znad + rhd(ji  ,jj+1,iktp1j) + zrhdtop_oce(ji  ,jj+1) )   & 
    559                &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    &  
    560                &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   & 
    561                &                                  + ( ziceload(ji,jj+1) - ziceload(ji,jj) )                             )  
     559            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj,iktp1i)                                    & 
     560               &                                    * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) )   & 
     561               &                                  - 0.5_wp * fse3w(ji,jj,ikt)                                         & 
     562               &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          & 
     563               &                                  + ( riceload(ji+1,jj) - riceload(ji,jj))                            )  
     564            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji,jj+1,iktp1j)                                    & 
     565               &                                    * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) )   & 
     566               &                                  - 0.5_wp * fse3w(ji,jj,ikt)                                         &  
     567               &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          & 
     568               &                                  + ( riceload(ji,jj+1) - riceload(ji,jj))                            )  
    562569            ! s-coordinate pressure gradient correction (=0 if z coordinate) 
    563570            zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
    … …  
    570577         END DO 
    571578      END DO 
    572 !==================================================================================      
    573 !===== Compute partial cell contribution for the top cell =========================  
    574 !================================================================================== 
    575       DO jj = 2, jpjm1 
    576          DO ji = fs_2, fs_jpim1   ! vector opt. 
    577             iku = miku(ji,jj) ;  
    578             zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp 
    579             ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
    580             ! u direction 
    581             IF ( iku .GT. 1 ) THEN 
    582                ! case iku 
    583                zhpi(ji,jj,iku)   =  zcoef0 / e1u(ji,jj) * ze3wu                                            & 
    584                       &                                 * ( rhd    (ji+1,jj,iku) + rhd   (ji,jj,iku)       & 
    585                       &                                   + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad ) 
    586                ! corrective term ( = 0 if z coordinate ) 
    587                zuap              = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj) 
    588                ! zhpi will be added in interior loop 
    589                ua(ji,jj,iku)     = ua(ji,jj,iku) + zuap 
    590                ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure   
    591                IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj)  = zhpi(ji,jj,iku) 
    592  
    593                ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps) 
    594                zhpiint        =  zcoef0 / e1u(ji,jj)                                                               &     
    595                   &           * (  fse3w(ji+1,jj  ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad)                          & 
    596                   &                                         + (rhd(ji+1,jj,iku  ) + znad) ) * tmask(ji+1,jj,iku)   & 
    597                   &              - fse3w(ji  ,jj  ,iku+1) * ( (rhd(ji  ,jj,iku+1) + znad)                          & 
    598                   &                                         + (rhd(ji  ,jj,iku  ) + znad) ) * tmask(ji  ,jj,iku)   ) 
    599                zhpi(ji,jj,iku+1) =  zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint  
    600             END IF 
    601                 
    602             ! v direction 
    603             ikv = mikv(ji,jj) 
    604             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)) 
    605             IF ( ikv .GT. 1 ) THEN 
    606                ! case ikv 
    607                zhpj(ji,jj,ikv)   =  zcoef0 / e2v(ji,jj) * ze3wv                                            & 
    608                      &                                  * ( rhd(ji,jj+1,ikv) + rhd   (ji,jj,ikv)           & 
    609                      &                                    + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad ) 
    610                ! corrective term ( = 0 if z coordinate ) 
    611                zvap              = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj) 
    612                ! zhpi will be added in interior loop 
    613                va(ji,jj,ikv)      = va(ji,jj,ikv) + zvap 
    614                ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure   
    615                IF (mbkv(ji,jj) == ikv + 1)  zpshpj(ji,jj)  =  zhpj(ji,jj,ikv)  
    616                 
    617                ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps) 
    618                zhpjint        =  zcoef0 / e2v(ji,jj)                                                              & 
    619                   &           * (  fse3w(ji  ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad)                         & 
    620                   &                                       + (rhd(ji,jj+1,ikv  ) + znad) ) * tmask(ji,jj+1,ikv)    & 
    621                   &              - fse3w(ji  ,jj  ,ikv+1) * ( (rhd(ji,jj  ,ikv+1) + znad)                         & 
    622                   &                                       + (rhd(ji,jj  ,ikv  ) + znad) ) * tmask(ji,jj  ,ikv)  ) 
    623                zhpj(ji,jj,ikv+1) =  zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint 
    624             END IF 
    625          END DO 
    626       END DO 
    627579 
    628580!==================================================================================      
    629581!===== Compute interior value =====================================================  
    630582!================================================================================== 
    631  
    632       DO jj = 2, jpjm1 
    633          DO ji = fs_2, fs_jpim1   ! vector opt. 
    634             iku=miku(ji,jj); ikv=mikv(ji,jj) 
    635             DO jk = 2, jpkm1 
     583      ! interior value (2=<jk=<jpkm1) 
     584      DO jk = 2, jpkm1 
     585         DO jj = 2, jpjm1 
     586            DO ji = fs_2, fs_jpim1   ! vector opt. 
    636587               ! hydrostatic pressure gradient along s-surfaces 
    637                ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc) 
    638                zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1)                                                              & 
    639                   &                            + zcoef0 / e1u(ji,jj)                                                           & 
    640                   &                                     * ( fse3w(ji+1,jj  ,jk) * ( (rhd(ji+1,jj,jk  ) + znad)                 & 
    641                   &                                                     + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1)   & 
    642                   &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji  ,jj,jk  ) + znad)                 & 
    643                   &                                                     + (rhd(ji  ,jj,jk-1) + znad) ) * tmask(ji  ,jj,jk-1)   )  
     588               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   & 
     589                  &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   & 
     590                  &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   ) 
     591               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   & 
     592                  &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   & 
     593                  &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   ) 
    644594               ! s-coordinate pressure gradient correction 
    645                ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc)  
    646                zuap = - zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )                    & 
    647                   &            * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1) 
    648                ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 
    649  
    650                ! hydrostatic pressure gradient along s-surfaces 
    651                ! zhpi is masked for the first wet cell  (contribution already done in the upper bloc) 
    652                zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1)                                                              & 
    653                   &                            + zcoef0 / e2v(ji,jj)                                                           & 
    654                   &                                     * ( fse3w(ji  ,jj+1,jk) * ( (rhd(ji,jj+1,jk  ) + znad)                 & 
    655                   &                                                     + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1)   & 
    656                   &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji,jj  ,jk  ) + znad)                 & 
    657                   &                                                     + (rhd(ji,jj  ,jk-1) + znad) ) * tmask(ji,jj  ,jk-1)   ) 
    658                ! s-coordinate pressure gradient correction 
    659                ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 
    660                zvap = - zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )                     & 
    661                   &            * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1) 
     595               zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
     596                  &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 
     597               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
     598                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
    662599               ! add to the general momentum trend 
    663                va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 
     600               ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 
     601               va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 
    664602            END DO 
    665603         END DO 
    666604      END DO 
    667  
    668 !==================================================================================      
    669 !===== Compute bottom cell contribution (partial cell) ============================  
    670 !================================================================================== 
    671  
    672       DO jj = 2, jpjm1 
    673          DO ji = 2, jpim1 
    674             iku = mbku(ji,jj) 
    675             ikv = mbkv(ji,jj) 
    676  
    677             IF (iku .GT. 1) THEN 
    678                ! remove old value (interior case) 
    679                zuap            = -zcoef0 * ( rhd   (ji+1,jj  ,iku) + rhd   (ji,jj,iku) + 2._wp * znad )   & 
    680                      &                   * ( fsde3w(ji+1,jj  ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj) 
    681                ua(ji,jj,iku)   = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 
    682                ! put new value 
    683                ! -zpshpi to avoid double contribution of the partial step in the top layer  
    684                zuap            = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj)  / e1u(ji,jj) 
    685                zhpi(ji,jj,iku) =  zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj)  
    686                ua(ji,jj,iku)   =  ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 
    687             END IF 
    688             ! v direction 
    689             IF (ikv .GT. 1) THEN 
    690                ! remove old value (interior case) 
    691                zvap            = -zcoef0 * ( rhd   (ji  ,jj+1,ikv) + rhd   (ji,jj,ikv) + 2._wp * znad )   & 
    692                      &                   * ( fsde3w(ji  ,jj+1,ikv) - fsde3w(ji,jj,ikv) )   / e2v(ji,jj) 
    693                va(ji,jj,ikv)   = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 
    694                ! put new value 
    695                ! -zpshpj to avoid double contribution of the partial step in the top layer  
    696                zvap            = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj)     / e2v(ji,jj) 
    697                zhpj(ji,jj,ikv) =  zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj)    
    698                va(ji,jj,ikv)   =  va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 
    699             END IF 
    700          END DO 
    701       END DO 
    702       
    703       ! set back to original density value into the ice shelf cell (maybe useless because it is masked) 
    704       rhd = zrhd 
    705       ! 
    706       CALL wrk_dealloc( jpi,jpj,2, ztstop) 
    707       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 
    708       CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 
     605      ! 
     606      CALL wrk_dealloc( jpi,jpj,2  , ztstop) 
     607      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj) 
     608      CALL wrk_dealloc( jpi,jpj    , zrhdtop_oce ) 
    709609      ! 
    710610   END SUBROUTINE hpg_isf 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r6006 r6012  
    9595      ! 
    9696      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     97      INTEGER  ::   ikt          ! local integers 
    9798      REAL(wp) ::   zue3a, zue3n, zue3b, zuf           ! local scalars 
    9899      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      - 
    … …  
    220221               ! Add volume filter correction: compatibility with tracer advection scheme 
    221222               ! => time filter + conservation correction (only at the first level) 
    222                IF ( nn_isf == 0) THEN   ! if no ice shelf melting 
     223               IF ( .NOT. ln_isf ) THEN   ! if no ice shelf melting 
    223224                  fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) & 
    224                                  &                                          -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 
     225                                 &                                         - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 
    225226               ELSE                     ! if ice shelf melting 
    226227                  DO jj = 1,jpj 
    227228                     DO ji = 1,jpi 
    228                         jk = mikt(ji,jj) 
    229                         fse3t_b(ji,jj,jk) = fse3t_b(ji,jj,jk) - atfp * rdt * r1_rau0                       & 
     229                        ikt = mikt(ji,jj) 
     230                        fse3t_b(ji,jj,ikt) = fse3t_b(ji,jj,ikt) - atfp * rdt * r1_rau0                     & 
    230231                                          &                          * ( (emp_b(ji,jj)    - emp(ji,jj)   ) & 
    231232                                          &                            - (rnf_b(ji,jj)    - rnf(ji,jj)   ) & 
    232                                           &                            + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,jk) 
     233                                          &                            + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,ikt) 
    233234                     END DO 
    234235                  END DO 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90

    r5930 r6012  
    231231      IF(  ioptio > 1  .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) )   & 
    232232           &   CALL ctl_stop( ' Choose only one surface pressure gradient scheme' ) 
    233       IF( ln_dynspg_ts .AND. ln_isfcav )   & 
    234            &   CALL ctl_stop( ' dynspg_ts not tested with ice shelf cavity ' ) 
     233      IF( ln_dynspg_exp .AND. ln_isfcav )   & 
     234           &   CALL ctl_stop( ' dynspg_exp not tested with ice shelf cavity ' ) 
    235235      ! 
    236236      IF( ln_dynspg_exp)   nspg =  0 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r5930 r6012  
    145145      INTEGER  ::   ji, jj, jk, jn        ! dummy loop indices 
    146146      INTEGER  ::   ikbu, ikbv, noffset      ! local integers 
     147      INTEGER  ::   iktu, iktv               ! local integers 
    147148      REAL(wp) ::   zraur, z1_2dt_b, z2dt_bf    ! local scalars 
    148149      REAL(wp) ::   zx1, zy1, zx2, zy2          !   -      - 
    … …  
    384385      DO jj = 2, jpjm1                          ! Remove coriolis term (and possibly spg) from barotropic trend 
    385386         DO ji = fs_2, fs_jpim1 
    386              zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1) 
    387              zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) 
     387             zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 
     388             zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 
    388389          END DO 
    389390      END DO  
    … …  
    413414      zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 
    414415      zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 
     416      !        
     417      !                                         ! Add top stress contribution from baroclinic velocities:       
     418      IF (ln_bt_fw) THEN 
     419         DO jj = 2, jpjm1 
     420            DO ji = fs_2, fs_jpim1   ! vector opt. 
     421               iktu = miku(ji,jj) 
     422               iktv = mikv(ji,jj) 
     423               zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities 
     424               zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 
     425            END DO 
     426         END DO 
     427      ELSE 
     428         DO jj = 2, jpjm1 
     429            DO ji = fs_2, fs_jpim1   ! vector opt. 
     430               iktu = miku(ji,jj) 
     431               iktv = mikv(ji,jj) 
     432               zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities 
     433               zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 
     434            END DO 
     435         END DO 
     436      ENDIF 
     437      ! 
     438      ! Note that the "unclipped" top friction parameter is used even with explicit drag 
     439      zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * tfrua(:,:) * zwx(:,:) 
     440      zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * tfrva(:,:) * zwy(:,:) 
    415441      !        
    416442      IF (ln_bt_fw) THEN                        ! Add wind forcing 
    … …  
    544570            DO jj = 2, jpjm1                                    ! Sea Surface Height at u- & v-points 
    545571               DO ji = 2, fs_jpim1   ! Vector opt. 
    546                   zwx(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj)     & 
     572                  zwx(ji,jj) = z1_2 * ssumask(ji,jj)  * r1_e1e2u(ji,jj)     & 
    547573                     &              * ( e1e2t(ji  ,jj) * zsshp2_e(ji  ,jj)  & 
    548574                     &              +   e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 
    549                   zwy(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj)     & 
     575                  zwy(ji,jj) = z1_2 * ssvmask(ji,jj)  * r1_e1e2v(ji,jj)     & 
    550576                     &              * ( e1e2t(ji,jj  ) * zsshp2_e(ji,jj  )  & 
    551577                     &              +   e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 
    … …  
    607633            END DO 
    608634         END DO 
    609          ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1) 
     635         ssha_e(:,:) = (  sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) )  ) * ssmask(:,:) 
    610636         CALL lbc_lnk( ssha_e, 'T',  1._wp ) 
    611637 
    … …  
    622648            DO jj = 2, jpjm1 
    623649               DO ji = 2, jpim1      ! NO Vector Opt. 
    624                   zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1)  * r1_e1e2u(ji,jj)  & 
    625                      &              * ( e1e2t(ji  ,jj  ) * ssha_e(ji  ,jj  ) & 
    626                      &              +   e1e2t(ji+1,jj  ) * ssha_e(ji+1,jj  ) ) 
    627                   zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1)  * r1_e1e2v(ji,jj)  & 
    628                      &              * ( e1e2t(ji  ,jj  ) * ssha_e(ji  ,jj  ) & 
    629                      &              +   e1e2t(ji  ,jj+1) * ssha_e(ji  ,jj+1) ) 
     650                  zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj)    & 
     651                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) & 
     652                     &              +   e1e2t(ji+1,jj  )  * ssha_e(ji+1,jj  ) ) 
     653                  zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj)    & 
     654                     &              * ( e1e2t(ji  ,jj  )  * ssha_e(ji  ,jj  ) & 
     655                     &              +   e1e2t(ji  ,jj+1)  * ssha_e(ji  ,jj+1) ) 
    630656               END DO 
    631657            END DO 
    … …  
    661687            DO jj = 2, jpjm1                             
    662688               DO ji = 2, jpim1 
    663                   zx1 = z1_2 * umask(ji  ,jj,1) *  r1_e1e2u(ji  ,jj)    & 
     689                  zx1 = z1_2 * ssumask(ji  ,jj) *  r1_e1e2u(ji  ,jj)    & 
    664690                     &      * ( e1e2t(ji  ,jj  ) * zsshp2_e(ji  ,jj)    & 
    665691                     &      +   e1e2t(ji+1,jj  ) * zsshp2_e(ji+1,jj  ) ) 
    666                   zy1 = z1_2 * vmask(ji  ,jj,1) *  r1_e1e2v(ji  ,jj  )  & 
     692                  zy1 = z1_2 * ssvmask(ji  ,jj) *  r1_e1e2v(ji  ,jj  )  & 
    667693                     &       * ( e1e2t(ji ,jj  ) * zsshp2_e(ji  ,jj  )  & 
    668694                     &       +   e1e2t(ji ,jj+1) * zsshp2_e(ji  ,jj+1) ) 
    … …  
    736762         zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 
    737763         ! 
     764         ! Add top stresses: 
     765         zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:) 
     766         zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 
     767         ! 
    738768         ! Surface pressure trend: 
    739769         DO jj = 2, jpjm1 
    … …  
    755785                            &                                 + zu_trd(ji,jj)   & 
    756786                            &                                 + zu_frc(ji,jj) ) &  
    757                             &   ) * umask(ji,jj,1) 
     787                            &   ) * ssumask(ji,jj) 
    758788 
    759789                  va_e(ji,jj) = (                                 vn_e(ji,jj)   & 
    … …  
    761791                            &                                 + zv_trd(ji,jj)   & 
    762792                            &                                 + zv_frc(ji,jj) ) & 
    763                             &   ) * vmask(ji,jj,1) 
    764                END DO 
    765             END DO 
    766  
    767          ELSE                 ! Flux form 
     793                            &   ) * ssvmask(ji,jj) 
     794               END DO 
     795            END DO 
     796 
     797         ELSE                                         ! Flux form 
    768798            DO jj = 2, jpjm1 
    769799               DO ji = fs_2, fs_jpim1   ! vector opt. 
    770800 
    771                   zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 
    772                   zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 
     801                  zhura = ssumask(ji,jj)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj)) 
     802                  zhvra = ssvmask(ji,jj)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj)) 
    773803 
    774804                  ua_e(ji,jj) = (                hu_e(ji,jj)  *   un_e(ji,jj)   &  
    … …  
    791821            hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 
    792822            hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 
    793             hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 
    794             hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 
     823            hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 
     824            hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) 
    795825            ! 
    796826         ENDIF 
    … …  
    827857            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:)  
    828858            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:)  
    829          ELSE                                ! Sum transports 
     859         ELSE                                              ! Sum transports 
    830860            ua_b  (:,:) = ua_b  (:,:) + za1 * ua_e  (:,:) * hu_e (:,:) 
    831861            va_b  (:,:) = va_b  (:,:) + za1 * va_e  (:,:) * hv_e (:,:) 
    … …  
    881911         END DO 
    882912         ! Save barotropic velocities not transport: 
    883          ua_b  (:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) ) 
    884          va_b  (:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) ) 
     913         ua_b  (:,:) =  ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 
     914         va_b  (:,:) =  va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 
    885915      ENDIF 
    886916      ! 
    … …  
    897927      ! 
    898928      IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN 
    899          IF ( Agrif_NbStepint().EQ.0 ) THEN 
     929         IF ( Agrif_NbStepint() == 0 ) THEN 
    900930            ub2_i_b(:,:) = 0.e0 
    901931            vb2_i_b(:,:) = 0.e0 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

    r5836 r6012  
    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) ) 
    93                zwvw(ji,jj,jk) = ( zww(ji  ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1)-vn(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    

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

    r5836 r6012  
    112112      !! 
    113113      INTEGER  ::   ji , jj , jk    ! dummy loop indices 
    114       INTEGER  ::   ii0, ii1, iku   ! temporary integer 
    115       INTEGER  ::   ij0, ij1, ikv   ! temporary integer 
     114      INTEGER  ::   ii0, ii1        ! temporary integer 
     115      INTEGER  ::   ij0, ij1        ! temporary integer 
    116116      REAL(wp) ::   zeps, zm1_g, zm1_2g, z1_16, zcofw, z1_slpmax ! local scalars 
    117117      REAL(wp) ::   zci, zfi, zau, zbu, zai, zbi   !   -      - 
    118118      REAL(wp) ::   zcj, zfj, zav, zbv, zaj, zbj   !   -      - 
    119119      REAL(wp) ::   zck, zfk,      zbw             !   -      - 
     120      REAL(wp) ::   zdepu, zdepv                   !   -      - 
     121      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zslpml_hmlpu, zslpml_hmlpv 
    120122      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwz, zww 
    121123      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdzr 
    … …  
    126128      ! 
    127129      CALL wrk_alloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 
     130      CALL wrk_alloc( jpi,jpj, zslpml_hmlpu, zslpml_hmlpv ) 
    128131 
    129132      zeps   =  1.e-20_wp        !==   Local constant initialization   ==! 
    … …  
    149152               zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 
    150153               zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 
     154            END DO 
     155         END DO 
     156      ENDIF 
     157      IF( ln_zps .AND. ln_isfcav ) THEN           ! partial steps correction at the bottom ocean level 
     158         DO jj = 1, jpjm1 
     159            DO ji = 1, jpim1 
     160               IF ( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj)  
     161               IF ( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj) 
    151162            END DO 
    152163         END DO 
    … …  
    171182      ! ===========================      | vslp = d/dj( prd ) / d/dz( prd ) 
    172183      ! 
     184      IF ( ln_isfcav ) THEN 
     185         DO jj = 2, jpjm1 
     186            DO ji = fs_2, fs_jpim1   ! vector opt. 
     187               zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / ( MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp)       & 
     188                  &                                  - 0.5_wp * ( risfdep(ji,jj) + risfdep(ji+1,jj  ) ) ) 
     189               zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / ( MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp)       & 
     190                  &                                  - 0.5_wp * ( risfdep(ji,jj) + risfdep(ji  ,jj+1) ) ) 
     191            END DO 
     192         END DO 
     193      ELSE 
     194         DO jj = 2, jpjm1 
     195            DO ji = fs_2, fs_jpim1   ! vector opt. 
     196               zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji+1,jj  ), 5._wp) 
     197               zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji  ,jj+1), 5._wp) 
     198            END DO 
     199         END DO 
     200      END IF 
     201 
    173202      DO jk = 2, jpkm1                            !* Slopes at u and v points 
    174203         DO jj = 2, jpjm1 
    … …  
    186215               zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 
    187216               zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 
    188                zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps )                                              & 
    189                   &                   + zfi  * uslpml(ji,jj)                                                     & 
    190                   &                          * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) )   & 
    191                   &                          / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk) 
    192                zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps )                                              & 
    193                   &                   + zfj  * vslpml(ji,jj)                                                     & 
    194                   &                          * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) )   & 
    195                   &                          / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk) 
     217               ! thickness of water column between surface and level k at u/v point 
     218               zdepu = 0.5_wp * ( ( fsdept (ji,jj,jk) + fsdept (ji+1,jj  ,jk) )                            & 
     219                                - ( risfdep(ji,jj)    + risfdep(ji+1,jj)    ) - fse3u(ji,jj,miku(ji,jj)) ) 
     220               zdepv = 0.5_wp * ( ( fsdept (ji,jj,jk) + fsdept (ji,jj+1,jk) )                              & 
     221                                - ( risfdep(ji,jj)    + risfdep(ji,jj+1)    ) - fse3v(ji,jj,mikv(ji,jj)) ) 
     222               ! 
     223               zwz(ji,jj,jk) = ( ( 1._wp - zfi) * zau / ( zbu - zeps )                                     & 
     224                  &                      + zfi  * zdepu * zslpml_hmlpu(ji,jj) ) * umask(ji,jj,jk) 
     225               zww(ji,jj,jk) = ( ( 1._wp - zfj) * zav / ( zbv - zeps )                                     & 
     226                  &                      + zfj  * zdepv * zslpml_hmlpv(ji,jj) ) * vmask(ji,jj,jk) 
    196227!!gm  modif to suppress omlmask.... (as in Griffies case) 
    197228!               !                                         ! jk must be >= ML level for zf=1. otherwise  zf=0. 
    … …  
    265296                  &      + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk  ) , zeps  ) * e2t(ji,jj) 
    266297               zai =    (  zgru (ji-1,jj,jk  ) + zgru (ji,jj,jk-1)           & 
    267                   &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci * tmask (ji,jj,jk) 
     298                  &      + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk  )   ) / zci * wmask (ji,jj,jk) 
    268299               zaj =    (  zgrv (ji,jj-1,jk  ) + zgrv (ji,jj,jk-1)           & 
    269                   &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj * tmask (ji,jj,jk) 
     300                  &      + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk  )   ) / zcj * wmask (ji,jj,jk) 
    270301               !                                        ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 
    271302               !                                        ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 
    … …  
    274305               !                                        ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 
    275306               zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) )   ! zfk=1 in the ML otherwise zfk=0 
    276                zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp ) 
    277                zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk  ) * tmask(ji,jj,jk) 
    278                zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk  ) * tmask(ji,jj,jk) 
     307               zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj) - fsdepw(ji,jj,mikt(ji,jj)), 10._wp ) 
     308               zwz(ji,jj,jk) = (  zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk  ) * wmask(ji,jj,jk) 
     309               zww(ji,jj,jk) = (  zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk  ) * wmask(ji,jj,jk) 
    279310 
    280311!!gm  modif to suppress omlmask....  (as in Griffies operator) 
    … …  
    340371      CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. ) 
    341372 
    342  
    343373      IF(ln_ctl) THEN 
    344374         CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp  - u : ', tab3d_2=vslp,  clinfo2=' v : ', kdim=jpk) 
    … …  
    347377      ! 
    348378      CALL wrk_dealloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 
     379      CALL wrk_dealloc( jpi,jpj, zslpml_hmlpu, zslpml_hmlpv ) 
    349380      ! 
    350381      IF( nn_timing == 1 )  CALL timing_stop('ldf_slp') 
    … …  
    486517                  ! 
    487518                  jk = nmln(ji,jj+jp) + 1 
    488                   IF( jk .GT. mbkt(ji,jj+jp) ) THEN  !ML reaches bottom 
     519                  IF( jk > mbkt(ji,jj+jp) ) THEN  !ML reaches bottom 
    489520                     ztj_mlb(ji   ,jj+jp,1-jp,kp) = 0.0_wp 
    490521                  ELSE 
    … …  
    699730            zcj = MAX(   vmask(ji,jj-1,ik  ) + vmask(ji,jj,ik  )           & 
    700731               &       + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps  ) * e2t(ji,jj) 
    701             zai =    (   p_gru(ji-1,jj,ik  ) + p_gru(ji,jj,ik)           & 
     732            zai =    (   p_gru(ji-1,jj,ik  ) + p_gru(ji,jj,ik)             & 
    702733               &       + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1  )  ) / zci  * tmask(ji,jj,ik) 
    703734            zaj =    (   p_grv(ji,jj-1,ik  ) + p_grv(ji,jj,ik  )           & 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r5836 r6012  
    4444   LOGICAL , PUBLIC ::   ln_dm2dc       !: Daily mean to Diurnal Cycle short wave (qsr) 
    4545   LOGICAL , PUBLIC ::   ln_rnf         !: runoffs / runoff mouths 
     46   LOGICAL , PUBLIC ::   ln_isf         !: ice shelf melting 
    4647   LOGICAL , PUBLIC ::   ln_ssr         !: Sea Surface restoring on SST and/or SSS       
    4748   LOGICAL , PUBLIC ::   ln_apr_dyn     !: Atmospheric pressure forcing used on dynamics (ocean & ice) 
    4849   INTEGER , PUBLIC ::   nn_ice         !: flag for ice in the surface boundary condition (=0/1/2/3) 
    49    INTEGER , PUBLIC ::   nn_isf         !: flag for isf in the surface boundary condition (=0/1/2/3/4)  
    5050   INTEGER , PUBLIC ::   nn_ice_embd    !: flag for levitating/embedding sea-ice in the ocean 
    5151   !                                             !: =0 levitating ice (no mass exchange, concentration/dilution effect) 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90

    r6006 r6012  
    3535   ! public in order to be able to output then  
    3636 
    37    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc   
    38    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf              !: net heat flux from ice shelf 
     37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc  !: before and now T & S isf contents [K.m/s & PSU.m/s]  
     38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf                  !: net heat flux from ice shelf      [W/m2] 
    3939   REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m] 
    40    LOGICAL , PUBLIC ::   ln_divisf                   !: flag to correct divergence  
    41    INTEGER , PUBLIC ::   nn_isfblk                   !:  
    42    INTEGER , PUBLIC ::   nn_gammablk                 !: 
    43    LOGICAL , PUBLIC ::   ln_conserve                 !: 
    44    REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient 
    45    REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient  
    46    REAL(wp), PUBLIC ::   rdivisf                     !: flag to test if fwf apply on divergence 
     40   INTEGER , PUBLIC ::   nn_isf                      !: flag to choose between explicit/param/specified   
     41   INTEGER , PUBLIC ::   nn_isfblk                   !: flag to choose the bulk formulation to compute the ice shelf melting 
     42   INTEGER , PUBLIC ::   nn_gammablk                 !: flag to choose how the exchange coefficient is computed 
     43   REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient [] 
     44   REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient [] 
    4745 
    4846   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rzisf_tbl              !:depth of calving front (shallowest point) nn_isf ==2/3 
    49    REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rhisf_tbl, rhisf_tbl_0 !:thickness of tbl 
     47   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rhisf_tbl, rhisf_tbl_0 !:thickness of tbl  [m] 
    5048   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  r1_hisf_tbl            !:1/thickness of tbl 
    5149   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ralpha                 !:proportion of bottom cell influenced by tbl  
    5250   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  risfLeff               !:effective length (Leff) BG03 nn_isf==2 
    5351   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point 
    54    INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base 
    55  
    56    REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ? 
    57    REAL(wp), PUBLIC, SAVE ::   kappa  = 1.54e-6_wp    ! phycst ? 
    58    REAL(wp), PUBLIC, SAVE ::   rhoisf = 920.0_wp      ! phycst ? 
    59    REAL(wp), PUBLIC, SAVE ::   tsurf  = -20.0_wp      ! phycst ? 
    60    REAL(wp), PUBLIC, SAVE ::   lfusisf= 0.334e6_wp    ! phycst ? 
     52   INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)      ::  misfkt, misfkb         !:Level of ice shelf base 
     53 
     54   REAL(wp), PUBLIC, SAVE ::   rcpi     = 2000.0_wp     ! specific heat of ice shelf             [J/kg/K] 
     55   REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    ! heat diffusivity through the ice-shelf [m2/s] 
     56   REAL(wp), PUBLIC, SAVE ::   rhoisf   = 920.0_wp      ! volumic mass of ice shelf              [kg/m3] 
     57   REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      ! air temperature on top of ice shelf    [C] 
     58   REAL(wp), PUBLIC, SAVE ::   rlfusisf = 0.334e6_wp    ! latent heat of fusion of ice shelf     [J/kg] 
    6159 
    6260!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3) 
    63    CHARACTER(len=100), PUBLIC ::   cn_dirisf  = './'    !: Root directory for location of ssr files 
    64    TYPE(FLD_N)       , PUBLIC ::   sn_qisf, sn_fwfisf     !: information about the runoff file to be read 
    65    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qisf, sf_fwfisf 
    66    TYPE(FLD_N)       , PUBLIC ::   sn_rnfisf              !: information about the runoff file to be read 
    67    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnfisf            
    68    TYPE(FLD_N)       , PUBLIC ::   sn_depmax_isf, sn_depmin_isf, sn_Leff_isf     !: information about the runoff file to be read 
     61   CHARACTER(len=100), PUBLIC           :: cn_dirisf  = './' !: Root directory for location of ssr files 
     62   TYPE(FLD_N)       , PUBLIC           :: sn_fwfisf         !: information about the isf melting file to be read 
     63   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf 
     64   TYPE(FLD_N)       , PUBLIC           :: sn_rnfisf         !: information about the isf melting param.   file to be read 
     65   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf            
     66   TYPE(FLD_N)       , PUBLIC           :: sn_depmax_isf     !: information about the grounding line depth file to be read 
     67   TYPE(FLD_N)       , PUBLIC           :: sn_depmin_isf     !: information about the calving   line depth file to be read 
     68   TYPE(FLD_N)       , PUBLIC           :: sn_Leff_isf       !: information about the effective length     file to be read 
    6969    
    7070   !! * Substitutions 
    … …  
    7777CONTAINS 
    7878  
    79    SUBROUTINE sbc_isf(kt) 
     79  SUBROUTINE sbc_isf(kt) 
    8080      !!--------------------------------------------------------------------- 
    81       !!                     ***  ROUTINE sbc_isf  *** 
    82       !!--------------------------------------------------------------------- 
    83       INTEGER, INTENT(in)          ::   kt         ! ocean time step 
    84       ! 
    85       INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror 
    86       INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer 
    87       REAL(wp)                     ::   rmin 
    88       REAL(wp)                     ::   zhk 
    89       REAL(wp)                     ::   zt_frz, zpress 
    90       CHARACTER(len=256)           ::   cfisf , cvarzisf, cvarhisf   ! name for isf file 
    91       CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file 
    92       CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale 
    93       INTEGER           ::   ios           ! Local integer output status for namelist read 
    94       !! 
    95       NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, ln_divisf, ln_conserve, rn_gammat0, rn_gammas0, nn_gammablk, & 
    96          &                 sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 
     81      !!                  ***  ROUTINE sbc_isf  *** 
     82      !! 
     83      !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf  
     84      !!              melting and freezing  
     85      !! 
     86      !! ** Method  :  4 parameterizations are available according to nn_isf  
     87      !!               nn_isf = 1 : Realistic ice_shelf formulation 
     88      !!                        2 : Beckmann & Goose parameterization 
     89      !!                        3 : Specified runoff in deptht (Mathiot & al. ) 
     90      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf 
     91      !!---------------------------------------------------------------------- 
     92      INTEGER, INTENT( in ) :: kt                   ! ocean time step 
     93      ! 
     94      INTEGER               :: ji, jj               ! loop index 
     95      REAL(wp), DIMENSION (:,:), POINTER :: zt_frz, zdep ! freezing temperature (zt_frz) at depth (zdep)  
    9796      !!--------------------------------------------------------------------- 
    9897      ! 
    … …  
    10099      IF( kt == nit000 ) THEN                   !  First call kt=nit000  ! 
    101100         !                                      ! ====================== ! 
    102          REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs  
    103          READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901) 
    104 901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp ) 
    105  
    106          REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs 
    107          READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 ) 
    108 902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp ) 
    109          IF(lwm) WRITE ( numond, namsbc_isf ) 
    110  
    111          IF ( lwp ) WRITE(numout,*) 
    112          IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf' 
    113          IF ( lwp ) WRITE(numout,*) '~~~~~~~~~' 
    114          IF ( lwp ) WRITE(numout,*) 'sbcisf :'  
    115          IF ( lwp ) WRITE(numout,*) '~~~~~~~~' 
    116          IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf 
    117          IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk 
    118          IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl 
    119          IF ( lwp ) WRITE(numout,*) '        ln_divisf   = ', ln_divisf  
    120          IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk  
    121          IF ( lwp ) WRITE(numout,*) '        rn_gammat0  = ', rn_gammat0   
    122          IF ( lwp ) WRITE(numout,*) '        rn_gammas0  = ', rn_gammas0   
    123          IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2  
    124          IF (ln_divisf) THEN       ! keep it in the namelist ??? used true anyway as for runoff ? (PM) 
    125             rdivisf = 1._wp 
    126          ELSE 
    127             rdivisf = 0._wp 
    128          END IF 
    129          ! 
    130          ! Allocate public variable 
    131          IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' ) 
    132          ! 
    133          ! initialisation 
    134          qisf(:,:)        = 0._wp  ; fwfisf(:,:) = 0._wp 
    135          risf_tsc(:,:,:)  = 0._wp 
    136          ! 
    137          ! define isf tbl tickness, top and bottom indice 
    138          IF      (nn_isf == 1) THEN 
    139             rhisf_tbl(:,:) = rn_hisf_tbl 
    140             misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv 
    141          ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN 
    142             ALLOCATE( sf_rnfisf(1), STAT=ierror ) 
    143             ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 
    144             CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 
    145  
    146             !: read effective lenght (BG03) 
    147             IF (nn_isf == 2) THEN 
    148                ! Read Data and save some integral values 
    149                CALL iom_open( sn_Leff_isf%clname, inum ) 
    150                cvarLeff  = 'soLeff'               !: variable name for Efficient Length scale 
    151                CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1) 
    152                CALL iom_close(inum) 
    153                ! 
    154                risfLeff = risfLeff*1000           !: convertion in m 
    155             END IF 
    156  
    157            ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth) 
    158             CALL iom_open( sn_depmax_isf%clname, inum ) 
    159             cvarhisf = TRIM(sn_depmax_isf%clvar) 
    160             CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base 
    161             CALL iom_close(inum) 
    162             ! 
    163             CALL iom_open( sn_depmin_isf%clname, inum ) 
    164             cvarzisf = TRIM(sn_depmin_isf%clvar) 
    165             CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base 
    166             CALL iom_close(inum) 
    167             ! 
    168             rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer 
    169  
    170            !! compute first level of the top boundary layer 
    171            DO ji = 1, jpi 
    172               DO jj = 1, jpj 
    173                   jk = 2 
    174                   DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdepw(ji,jj,jk) < rzisf_tbl(ji,jj) ) ;  jk = jk + 1 ;  END DO 
    175                   misfkt(ji,jj) = jk-1 
    176                END DO 
    177             END DO 
    178  
    179          ELSE IF ( nn_isf == 4 ) THEN 
    180             ! as in nn_isf == 1 
    181             rhisf_tbl(:,:) = rn_hisf_tbl 
    182             misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv 
    183              
    184             ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 
    185             ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror ) 
    186             ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 
    187             ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) ) 
    188             CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 
    189             !CALL fld_fill( sf_qisf  , (/ sn_qisf   /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data'       , 'namsbc_isf' ) 
    190          END IF 
    191           
    192          ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 
    193          rhisf_tbl_0(:,:) = rhisf_tbl(:,:) 
    194          DO jj = 1,jpj 
    195             DO ji = 1,jpi 
    196                ikt = misfkt(ji,jj) 
    197                ikb = misfkt(ji,jj) 
    198                ! thickness of boundary layer at least the top level thickness 
    199                rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt)) 
    200  
    201                ! determine the deepest level influenced by the boundary layer 
    202                ! test on tmask useless ????? 
    203                DO jk = ikt, mbkt(ji,jj) 
    204                   IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
    205                END DO 
    206                rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
    207                misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl 
    208                r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 
    209  
    210                zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1 
    211                ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer 
    212             END DO 
    213          END DO 
     101         CALL sbc_isf_init 
     102      !                                         ! ---------------------------------------- ! 
     103      ELSE                                      !          Swap of forcing fields          ! 
     104         !                                      ! ---------------------------------------- ! 
     105         fwfisf_b  (:,:  ) = fwfisf  (:,:  )    ! Swap the ocean forcing fields except at nit000 
     106         risf_tsc_b(:,:,:) = risf_tsc(:,:,:)    ! where before fields are set at the end of the routine 
    214107         ! 
    215108      END IF 
    216109 
    217       !                                            ! ---------------------------------------- ! 
    218       IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          ! 
    219          !                                         ! ---------------------------------------- ! 
    220          fwfisf_b  (:,:  ) = fwfisf  (:,:  )               ! Swap the ocean forcing fields except at nit000 
    221          risf_tsc_b(:,:,:) = risf_tsc(:,:,:)               ! where before fields are set at the end of the routine 
    222          ! 
    223       ENDIF 
    224  
    225110      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN 
    226  
    227  
    228          ! compute salf and heat flux 
    229          IF (nn_isf == 1) THEN 
    230             ! realistic ice shelf formulation 
     111         ! allocation 
     112         CALL wrk_alloc( jpi,jpj, zt_frz, zdep  ) 
     113 
     114         ! compute salt and heat flux 
     115         SELECT CASE ( nn_isf ) 
     116         CASE ( 1 )    ! realistic ice shelf formulation 
    231117            ! compute T/S/U/V for the top boundary layer 
    232118            CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T') 
    … …  
    242128            CALL sbc_isf_cav (kt) 
    243129 
    244          ELSE IF (nn_isf == 2) THEN 
    245             ! Beckmann and Goosse parametrisation  
     130         CASE ( 2 )    ! Beckmann and Goosse parametrisation  
    246131            stbl(:,:)   = soce 
    247132            CALL sbc_isf_bg03(kt) 
    248133 
    249          ELSE IF (nn_isf == 3) THEN 
    250             ! specified runoff in depth (Mathiot et al., XXXX in preparation) 
     134         CASE ( 3 )    ! specified runoff in depth (Mathiot et al., XXXX in preparation) 
    251135            CALL fld_read ( kt, nn_fsbc, sf_rnfisf   ) 
    252             fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting)  
    253             qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux 
     136            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fwf flux from the isf (fwfisf <0 mean melting)  
     137            qisf(:,:)   = fwfisf(:,:) * rlfusisf             ! heat flux 
    254138            stbl(:,:)   = soce 
    255139 
    256          ELSE IF (nn_isf == 4) THEN 
    257             ! specified fwf and heat flux forcing beneath the ice shelf 
     140         CASE ( 4 )    ! specified fwf and heat flux forcing beneath the ice shelf 
    258141            CALL fld_read ( kt, nn_fsbc, sf_fwfisf   ) 
    259             !CALL fld_read ( kt, nn_fsbc, sf_qisf   ) 
    260             fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf 
    261             qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux 
    262             !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux 
     142            fwfisf(:,:) = - sf_fwfisf(1)%fnow(:,:,1)           ! fwf  flux from the isf (fwfisf <0 mean melting) 
     143            qisf(:,:)   = fwfisf(:,:) * rlfusisf               ! heat flux 
    263144            stbl(:,:)   = soce 
    264145 
    265          END IF 
     146         END SELECT 
     147 
    266148         ! compute tsc due to isf 
    267          ! WARNING water add at temp = 0C, correction term is added, maybe better here but need a 3D variable). 
    268 !         zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 
    269          zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress ) 
    270          risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - rdivisf * fwfisf(:,:) * zt_frz * r1_rau0 ! 
     149         ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU. 
     150         ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rau0). 
     151         ! compute freezing point beneath ice shelf (or top cell if nn_isf = 3) 
     152         DO jj = 1,jpj 
     153            DO ji = 1,jpi 
     154               zdep(ji,jj)=fsdepw_n(ji,jj,misfkt(ji,jj)) 
     155            END DO 
     156         END DO 
     157         CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) ) 
    271158          
    272          ! salt effect already take into account in vertical advection 
    273          risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0 
    274  
    275          ! output 
    276          IF( iom_use('qisf'  ) )   CALL iom_put('qisf'  , qisf) 
    277          IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce ) 
    278  
    279          ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now 
    280          fwfisf(:,:) = rdivisf * fwfisf(:,:)          
    281   
     159         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 ! 
     160         risf_tsc(:,:,jp_sal) = 0.0_wp 
     161 
    282162         ! lbclnk 
    283163         CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.) 
    284164         CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.) 
    285          CALL lbc_lnk(fwfisf(:,:)   ,'T',1.) 
    286          CALL lbc_lnk(qisf(:,:)     ,'T',1.) 
    287  
    288          IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    ! 
     165         CALL lbc_lnk(fwfisf(:,:)         ,'T',1.) 
     166         CALL lbc_lnk(qisf(:,:)           ,'T',1.) 
     167 
     168         IF( kt == nit000 ) THEN                         !   set the forcing field at nit000 - 1    ! 
    289169            IF( ln_rstart .AND.    &                     ! Restart: read in restart file 
    290170                 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN 
    … …  
    297177               risf_tsc_b(:,:,:)= risf_tsc(:,:,:) 
    298178            END IF 
    299          ENDIF 
     179         END IF 
    300180         !  
     181         ! output 
     182         CALL iom_put('qisf'  , qisf) 
     183         CALL iom_put('fwfisf', fwfisf) 
     184 
     185         ! deallocation 
     186         CALL wrk_dealloc( jpi,jpj, zt_frz, zdep  ) 
    301187      END IF 
    302188      !   
    … …  
    317203               &    STAT= sbc_isf_alloc ) 
    318204         ! 
    319          IF( lk_mpp                  )   CALL mpp_sum ( sbc_isf_alloc ) 
     205         IF( lk_mpp             )   CALL mpp_sum ( sbc_isf_alloc ) 
    320206         IF( sbc_isf_alloc /= 0 )   CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.') 
    321207         ! 
    322       ENDIF 
     208      END IF 
    323209  END FUNCTION 
    324210 
    325  
    326    SUBROUTINE sbc_isf_bg03(kt) 
    327       !!========================================================================== 
    328       !!                 *** SUBROUTINE sbcisf_bg03  *** 
    329       !! add net heat and fresh water flux from ice shelf melting 
    330       !! into the adjacent ocean using the parameterisation by 
    331       !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 
    332       !!     interaction for climate models", Ocean Modelling 5(2003) 157-170. 
    333       !!  (hereafter BG) 
    334       !!========================================================================== 
    335       !!---------------------------------------------------------------------- 
    336       !!   sbc_isf_bg03      : routine called from sbcmod 
    337       !!---------------------------------------------------------------------- 
    338       !! 
    339       !! ** Purpose   :   Add heat and fresh water fluxes due to ice shelf melting 
    340       !! ** Reference :   Beckmann et Goosse, 2003, Ocean Modelling 
    341       !! 
     211  SUBROUTINE sbc_isf_init 
     212      !!--------------------------------------------------------------------- 
     213      !!                  ***  ROUTINE sbc_isf_init  *** 
     214      !! 
     215      !! ** Purpose : Initialisation of variables for iceshelf fluxes formulation 
     216      !! 
     217      !! ** Method  :  4 parameterizations are available according to nn_isf  
     218      !!               nn_isf = 1 : Realistic ice_shelf formulation 
     219      !!                        2 : Beckmann & Goose parameterization 
     220      !!                        3 : Specified runoff in deptht (Mathiot & al. ) 
     221      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf 
     222      !!---------------------------------------------------------------------- 
     223      INTEGER               :: ji, jj, jk           ! loop index 
     224      INTEGER               :: ik                   ! current level index 
     225      INTEGER               :: ikt, ikb             ! top and bottom level of the isf boundary layer 
     226      INTEGER               :: inum, ierror 
     227      INTEGER               :: ios                  ! Local integer output status for namelist read 
     228      REAL(wp)              :: zhk 
     229      CHARACTER(len=256)    :: cvarzisf, cvarhisf   ! name for isf file 
     230      CHARACTER(LEN=32 )    :: cvarLeff             ! variable name for efficient Length scale 
     231      !!---------------------------------------------------------------------- 
     232      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf, & 
     233                         & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf 
     234      !!---------------------------------------------------------------------- 
     235 
     236      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs  
     237      READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901) 
     238901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp ) 
     239 
     240      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs 
     241      READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 ) 
     242902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp ) 
     243      IF(lwm) WRITE ( numond, namsbc_isf ) 
     244 
     245      IF ( lwp ) WRITE(numout,*) 
     246      IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf' 
     247      IF ( lwp ) WRITE(numout,*) '~~~~~~~~~' 
     248      IF ( lwp ) WRITE(numout,*) 'sbcisf :'  
     249      IF ( lwp ) WRITE(numout,*) '~~~~~~~~' 
     250      IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf 
     251      IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk 
     252      IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl 
     253      IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk  
     254      IF ( lwp ) WRITE(numout,*) '        rn_gammat0  = ', rn_gammat0   
     255      IF ( lwp ) WRITE(numout,*) '        rn_gammas0  = ', rn_gammas0   
     256      IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2  
     257      ! 
     258      ! Allocate public variable 
     259      IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' ) 
     260      ! 
     261      ! initialisation 
     262      qisf(:,:)        = 0._wp  ; fwfisf  (:,:) = 0._wp 
     263      risf_tsc(:,:,:)  = 0._wp  ; fwfisf_b(:,:) = 0._wp 
     264      ! 
     265      ! define isf tbl tickness, top and bottom indice 
     266      SELECT CASE ( nn_isf ) 
     267      CASE ( 1 )  
     268         rhisf_tbl(:,:) = rn_hisf_tbl 
     269         misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv 
     270 
     271      CASE ( 2 , 3 ) 
     272         ALLOCATE( sf_rnfisf(1), STAT=ierror ) 
     273         ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) ) 
     274         CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 
     275 
     276         !  read effective lenght (BG03) 
     277         IF (nn_isf == 2) THEN 
     278            CALL iom_open( sn_Leff_isf%clname, inum ) 
     279            cvarLeff = TRIM(sn_Leff_isf%clvar) 
     280            CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1) 
     281            CALL iom_close(inum) 
     282            ! 
     283            risfLeff = risfLeff*1000.0_wp           !: convertion in m 
     284         END IF 
     285 
     286         ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth) 
     287         CALL iom_open( sn_depmax_isf%clname, inum ) 
     288         cvarhisf = TRIM(sn_depmax_isf%clvar) 
     289         CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base 
     290         CALL iom_close(inum) 
     291         ! 
     292         CALL iom_open( sn_depmin_isf%clname, inum ) 
     293         cvarzisf = TRIM(sn_depmin_isf%clvar) 
     294         CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base 
     295         CALL iom_close(inum) 
     296         ! 
     297         rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer 
     298 
     299         !! compute first level of the top boundary layer 
     300         DO ji = 1, jpi 
     301            DO jj = 1, jpj 
     302                ik = 2 
     303                DO WHILE ( ik <= mbkt(ji,jj) .AND. fsdepw(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO 
     304                misfkt(ji,jj) = ik-1 
     305            END DO 
     306         END DO 
     307 
     308      CASE ( 4 )  
     309         ! as in nn_isf == 1 
     310         rhisf_tbl(:,:) = rn_hisf_tbl 
     311         misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv 
     312          
     313         ! load variable used in fldread (use for temporal interpolation of isf fwf forcing) 
     314         ALLOCATE( sf_fwfisf(1), STAT=ierror ) 
     315         ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) ) 
     316         CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' ) 
     317 
     318      END SELECT 
     319          
     320      rhisf_tbl_0(:,:) = rhisf_tbl(:,:) 
     321 
     322      ! compute bottom level of isf tbl and thickness of tbl below the ice shelf 
     323      DO jj = 1,jpj 
     324         DO ji = 1,jpi 
     325            ikt = misfkt(ji,jj) 
     326            ikb = misfkt(ji,jj) 
     327            ! thickness of boundary layer at least the top level thickness 
     328            rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt)) 
     329 
     330            ! determine the deepest level influenced by the boundary layer 
     331            DO jk = ikt+1, mbkt(ji,jj) 
     332               IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
     333            END DO 
     334            rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness. 
     335            misfkb(ji,jj) = ikb                                                   ! last wet level of the tbl 
     336            r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 
     337 
     338            zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 
     339            ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer 
     340         END DO 
     341      END DO 
     342 
     343  END SUBROUTINE sbc_isf_init 
     344 
     345  SUBROUTINE sbc_isf_bg03(kt) 
     346      !!--------------------------------------------------------------------- 
     347      !!                  ***  ROUTINE sbc_isf_bg03  *** 
     348      !! 
     349      !! ** Purpose : add net heat and fresh water flux from ice shelf melting 
     350      !!          into the adjacent ocean 
     351      !! 
     352      !! ** Method  :   See reference 
     353      !! 
     354      !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean 
     355      !!         interaction for climate models", Ocean Modelling 5(2003) 157-170. 
     356      !!         (hereafter BG) 
    342357      !! History : 
    343       !!      !  06-02  (C. Wang) Original code 
     358      !!         06-02  (C. Wang) Original code 
    344359      !!---------------------------------------------------------------------- 
    345360      INTEGER, INTENT ( in ) :: kt 
    346361      ! 
    347     INTEGER :: ji, jj, jk, jish  !temporary integer 
    348     INTEGER :: ijkmin 
    349     INTEGER :: ii, ij, ik  
    350     INTEGER :: inum 
    351  
    352     REAL(wp) :: zt_sum      ! sum of the temperature between 200m and 600m 
    353     REAL(wp) :: zt_ave      ! averaged temperature between 200m and 600m 
    354     REAL(wp) :: zt_frz      ! freezing point temperature at depth z 
    355     REAL(wp) :: zpress      ! pressure to compute the freezing point in depth 
    356      
    357     !!---------------------------------------------------------------------- 
    358     IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 
    359      ! 
    360     DO ji = 1, jpi 
    361        DO jj = 1, jpj 
    362           ik = misfkt(ji,jj) 
    363           !! Initialize arrays to 0 (each step) 
    364           zt_sum = 0.e0_wp 
    365           IF ( ik .GT. 1 ) THEN 
    366     ! 3. -----------the average temperature between 200m and 600m --------------------- 
    367              DO jk = misfkt(ji,jj),misfkb(ji,jj) 
    368              ! Calculate freezing temperature 
    369                 CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, fsdept(ji,jj,ik))  
    370                 zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik)  ! sum temp 
    371              ENDDO 
    372              zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 
    373      
    374     ! 4. ------------Net heat flux and fresh water flux due to the ice shelf 
    375           ! For those corresponding to zonal boundary     
    376              qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  & 
    377                          & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,ik)  
     362      INTEGER  :: ji, jj, jk ! dummy loop index 
     363      INTEGER  :: ik         ! current level 
     364      REAL(wp) :: zt_sum     ! sum of the temperature between 200m and 600m 
     365      REAL(wp) :: zt_ave     ! averaged temperature between 200m and 600m 
     366      REAL(wp) :: zt_frz     ! freezing point temperature at depth z 
     367      REAL(wp) :: zpress     ! pressure to compute the freezing point in depth 
     368      !!---------------------------------------------------------------------- 
     369 
     370      IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03') 
     371      ! 
     372      DO ji = 1, jpi 
     373         DO jj = 1, jpj 
     374            ik = misfkt(ji,jj) 
     375            !! Initialize arrays to 0 (each step) 
     376            zt_sum = 0.e0_wp 
     377            IF ( ik > 1 ) THEN 
     378               ! 1. -----------the average temperature between 200m and 600m --------------------- 
     379               DO jk = misfkt(ji,jj),misfkb(ji,jj) 
     380                  ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!) 
     381                  ! after verif with UNESCO, wrong sign in BG eq. 2 
     382                  ! Calculate freezing temperature 
     383                  CALL eos_fzp(stbl(ji,jj), zt_frz, zpress)  
     384                  zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * fse3t(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp 
     385               END DO 
     386               zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value 
     387               ! 2. ------------Net heat flux and fresh water flux due to the ice shelf 
     388               ! For those corresponding to zonal boundary     
     389               qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  & 
     390                           & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk) 
    378391              
    379              fwfisf(ji,jj) = qisf(ji,jj) / lfusisf          !fresh water flux kg/(m2s)                   
    380              fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 
    381              !add to salinity trend 
    382           ELSE 
    383              qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 
    384           END IF 
    385        END DO 
    386     END DO 
    387     ! 
    388     IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03') 
     392               fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf          !fresh water flux kg/(m2s)                   
     393               fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) ) 
     394               !add to salinity trend 
     395            ELSE 
     396               qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp 
     397            END IF 
     398         END DO 
     399      END DO 
     400      ! 
     401      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03') 
    389402      ! 
    390403  END SUBROUTINE sbc_isf_bg03 
    391404 
    392  
    393    SUBROUTINE sbc_isf_cav( kt ) 
     405  SUBROUTINE sbc_isf_cav( kt ) 
    394406      !!--------------------------------------------------------------------- 
    395407      !!                     ***  ROUTINE sbc_isf_cav  *** 
    … …  
    406418      INTEGER, INTENT(in)          ::   kt         ! ocean time step 
    407419      ! 
    408       LOGICAL :: ln_isomip = .true. 
    409       REAL(wp), DIMENSION(:,:), POINTER       ::   zfrz,zpress,zti 
    410       REAL(wp), DIMENSION(:,:), POINTER       ::   zgammat2d, zgammas2d  
    411       !REAL(wp), DIMENSION(:,:), POINTER ::   zqisf, zfwfisf 
     420      INTEGER  ::   ji, jj     ! dummy loop indices 
     421      INTEGER  ::   nit 
    412422      REAL(wp) ::   zlamb1, zlamb2, zlamb3 
    413423      REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7 
    414424      REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac 
    415       REAL(wp) ::   zfwflx, zhtflx, zhtflx_b 
    416       REAL(wp) ::   zgammat, zgammas 
    417       REAL(wp) ::   zeps = 1.e-20_wp        !==   Local constant initialization   ==! 
    418       INTEGER  ::   ji, jj     ! dummy loop indices 
    419       INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers 
    420       INTEGER  ::   ierror     ! return error code 
    421       LOGICAL  ::   lit=.TRUE. 
    422       INTEGER  ::   nit 
     425      REAL(wp) ::   zeps = 1.e-20_wp         
     426      REAL(wp) ::   zerr 
     427      REAL(wp), DIMENSION(:,:), POINTER ::   zfrz 
     428      REAL(wp), DIMENSION(:,:), POINTER ::   zgammat, zgammas  
     429      REAL(wp), DIMENSION(:,:), POINTER ::   zfwflx, zhtflx, zhtflx_b 
     430      LOGICAL  ::   lit 
    423431      !!--------------------------------------------------------------------- 
    424       ! 
    425       ! coeficient for linearisation of tfreez 
    426       zlamb1=-0.0575 
    427       zlamb2=0.0901 
    428       zlamb3=-7.61e-04 
     432      ! coeficient for linearisation of potential tfreez 
     433      ! Crude approximation for pressure (but commonly used) 
     434      zlamb1 =-0.0573_wp 
     435      zlamb2 = 0.0832_wp 
     436      zlamb3 =-7.53e-08_wp * grav * rau0 
    429437      IF( nn_timing == 1 )  CALL timing_start('sbc_isf_cav') 
    430438      ! 
    431       CALL wrk_alloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d ) 
    432  
    433       zcfac=0.0_wp  
    434       IF (ln_conserve)  zcfac=1.0_wp 
    435       zpress(:,:)=0.0_wp 
    436       zgammat2d(:,:)=0.0_wp 
    437       zgammas2d(:,:)=0.0_wp 
    438       ! 
    439       ! 
    440       DO jj = 1, jpj 
    441          DO ji = 1, jpi 
    442             ! Crude approximation for pressure (but commonly used) 
    443             ! 1e-04 to convert from Pa to dBar 
    444             zpress(ji,jj)=grav*rau0*fsdepw(ji,jj,mikt(ji,jj))*1.e-04 
    445             ! 
    446          END DO 
     439      CALL wrk_alloc( jpi,jpj, zfrz  , zgammat, zgammas  ) 
     440      CALL wrk_alloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b ) 
     441 
     442      ! initialisation 
     443      zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0 
     444      zhtflx (:,:) = 0.0_wp     ; zhtflx_b(:,:) = 0.0_wp     
     445      zfwflx (:,:) = 0.0_wp 
     446 
     447      ! compute ice shelf melting 
     448      nit = 1 ; lit = .TRUE. 
     449      DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine 
     450         SELECT CASE ( nn_isfblk ) 
     451         CASE ( 1 )   !  ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006) 
     452            ! Calculate freezing temperature 
     453            CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) ) 
     454 
     455            ! compute gammat every where (2d) 
     456            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 
     457             
     458            ! compute upward heat flux zhtflx and upward water flux zwflx 
     459            DO jj = 1, jpj 
     460               DO ji = 1, jpi 
     461                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj)) 
     462                  zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf 
     463               END DO 
     464            END DO 
     465 
     466            ! Compute heat flux and upward fresh water flux 
     467            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 
     468            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 
     469 
     470         CASE ( 2 )  ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 
     471            ! compute gammat every where (2d) 
     472            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx) 
     473 
     474            ! compute upward heat flux zhtflx and upward water flux zwflx 
     475            ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015) 
     476            DO jj = 1, jpj 
     477               DO ji = 1, jpi 
     478                  ! compute coeficient to solve the 2nd order equation 
     479                  zeps1 = rcp*rau0*zgammat(ji,jj) 
     480                  zeps2 = rlfusisf*rau0*zgammas(ji,jj) 
     481                  zeps3 = rhoisf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps) 
     482                  zeps4 = zlamb2+zlamb3*risfdep(ji,jj) 
     483                  zeps6 = zeps4-ttbl(ji,jj) 
     484                  zeps7 = zeps4-tsurf 
     485                  zaqe  = zlamb1 * (zeps1 + zeps3) 
     486                  zaqer = 0.5_wp/MIN(zaqe,-zeps) 
     487                  zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2 
     488                  zcqe  = zeps2*stbl(ji,jj) 
     489                  zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe                
     490 
     491                  ! Presumably zdis can never be negative because gammas is very small compared to gammat 
     492                  ! compute s freeze 
     493                  zsfrz=(-zbqe-SQRT(zdis))*zaqer 
     494                  IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer 
     495 
     496                  ! compute t freeze (eq. 22) 
     497                  zfrz(ji,jj)=zeps4+zlamb1*zsfrz 
     498   
     499                  ! zfwflx is upward water flux 
     500                  ! zhtflx is upward heat flux (out of ocean) 
     501                  ! compute the upward water and heat flux (eq. 28 and eq. 29) 
     502                  zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) 
     503                  zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) )  
     504               END DO 
     505            END DO 
     506 
     507            ! compute heat and water flux 
     508            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 
     509            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:) 
     510 
     511         END SELECT 
     512 
     513         ! define if we need to iterate (nn_gammablk 0/1 do not need iteration) 
     514         IF ( nn_gammablk <  2 ) THEN ; lit = .FALSE. 
     515         ELSE                            
     516            ! check total number of iteration 
     517            IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 
     518            ELSE                 ; nit = nit + 1 
     519            END IF 
     520 
     521            ! compute error between 2 iterations 
     522            ! if needed save gammat and compute zhtflx_b for next iteration 
     523            zerr = MAXVAL(ABS(zhtflx-zhtflx_b)) 
     524            IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE. 
     525            ELSE                        ; zhtflx_b(:,:) = zhtflx(:,:) 
     526            END IF 
     527         END IF 
    447528      END DO 
    448  
    449 ! Calculate in-situ temperature (ref to surface) 
    450       zti(:,:)=tinsitu( ttbl, stbl, zpress ) 
    451 ! Calculate freezing temperature 
    452       CALL eos_fzp( sss_m(:,:), zfrz(:,:), zpress ) 
    453  
    454        
    455       zhtflx=0._wp ; zfwflx=0._wp 
    456       IF (nn_isfblk == 1) THEN 
    457          DO jj = 1, jpj 
    458             DO ji = 1, jpi 
    459                IF (mikt(ji,jj) > 1 ) THEN 
    460                   nit = 1; lit = .TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp 
    461                   DO WHILE ( lit ) 
    462 ! compute gamma 
    463                      CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit) 
    464 ! zhtflx is upward heat flux (out of ocean) 
    465                      zhtflx = zgammat*rcp*rau0*(zti(ji,jj)-zfrz(ji,jj)) 
    466 ! zwflx is upward water flux 
    467                      zfwflx = - zhtflx/lfusisf 
    468 ! test convergence and compute gammat 
    469                      IF ( (zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE. 
    470  
    471                      nit = nit + 1 
    472                      IF (nit .GE. 100) CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 
    473  
    474 ! save gammat and compute zhtflx_b 
    475                      zgammat2d(ji,jj)=zgammat 
    476                      zhtflx_b = zhtflx 
    477                   END DO 
    478  
    479                   qisf(ji,jj) = - zhtflx 
    480 ! For genuine ISOMIP protocol this should probably be something like 
    481                   fwfisf(ji,jj) = zfwflx  * ( soce / MAX(stbl(ji,jj),zeps)) 
    482                ELSE 
    483                   fwfisf(ji,jj) = 0._wp 
    484                   qisf(ji,jj)   = 0._wp 
    485                END IF 
    486             ! 
    487             END DO 
    488          END DO 
    489  
    490       ELSE IF (nn_isfblk == 2 ) THEN 
    491  
    492 ! More complicated 3 equation thermodynamics as in MITgcm 
    493          DO jj = 2, jpj 
    494             DO ji = 2, jpi 
    495                IF (mikt(ji,jj) > 1 ) THEN 
    496                   nit=1; lit=.TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp; zhtflx=0._wp 
    497                   DO WHILE ( lit ) 
    498                      CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit) 
    499  
    500                      zeps1=rcp*rau0*zgammat 
    501                      zeps2=lfusisf*rau0*zgammas 
    502                      zeps3=rhoisf*rcpi*kappa/MAX(risfdep(ji,jj),zeps) 
    503                      zeps4=zlamb2+zlamb3*risfdep(ji,jj) 
    504                      zeps6=zeps4-zti(ji,jj) 
    505                      zeps7=zeps4-tsurf 
    506                      zaqe=zlamb1 * (zeps1 + zeps3) 
    507                      zaqer=0.5/MIN(zaqe,-zeps) 
    508                      zbqe=zeps1*zeps6+zeps3*zeps7-zeps2 
    509                      zcqe=zeps2*stbl(ji,jj) 
    510                      zdis=zbqe*zbqe-4.0*zaqe*zcqe                
    511 ! Presumably zdis can never be negative because gammas is very small compared to gammat 
    512                      zsfrz=(-zbqe-SQRT(zdis))*zaqer 
    513                      IF (zsfrz .lt. 0.0) zsfrz=(-zbqe+SQRT(zdis))*zaqer 
    514                      zfrz(ji,jj)=zeps4+zlamb1*zsfrz 
    515    
    516 ! zfwflx is upward water flux (MAX usefull if kappa = 0 
    517                      zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps) ) 
    518 ! zhtflx is upward heat flux (out of ocean) 
    519 ! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value 
    520                      zhtflx = ( zgammat*rau0 - zcfac*zfwflx ) * rcp * (zti(ji,jj) - zfrz(ji,jj) )  
    521 ! zwflx is upward water flux 
    522 ! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz 
    523                      zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj) 
    524 ! test convergence and compute gammat 
    525                      IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE. 
    526  
    527                      nit = nit + 1 
    528                      IF (nit .GE. 51) THEN 
    529                         WRITE(numout,*) "sbcisf : too many iteration ... ", & 
    530                             &  zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea 
    531                         CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 
    532                      END IF 
    533 ! save gammat and compute zhtflx_b 
    534                      zgammat2d(ji,jj)=zgammat 
    535                      zgammas2d(ji,jj)=zgammas 
    536                      zhtflx_b = zhtflx 
    537  
    538                   END DO 
    539 ! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value 
    540                   qisf(ji,jj) = - zhtflx  
    541 ! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz 
    542                   fwfisf(ji,jj) = zfwflx  
    543                ELSE 
    544                   fwfisf(ji,jj) = 0._wp 
    545                   qisf(ji,jj)   = 0._wp 
    546                ENDIF 
    547                ! 
    548             END DO 
    549          END DO 
    550       ENDIF 
    551       ! lbclnk 
    552       CALL lbc_lnk(zgammas2d(:,:),'T',1.) 
    553       CALL lbc_lnk(zgammat2d(:,:),'T',1.) 
    554       ! output 
    555       CALL iom_put('isfgammat', zgammat2d) 
    556       CALL iom_put('isfgammas', zgammas2d) 
    557       ! 
    558       CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d ) 
     529      ! 
     530      CALL iom_put('isfgammat', zgammat) 
     531      CALL iom_put('isfgammas', zgammas) 
     532      !  
     533      CALL wrk_dealloc( jpi,jpj, zfrz  , zgammat, zgammas  ) 
     534      CALL wrk_dealloc( jpi,jpj, zfwflx, zhtflx , zhtflx_b ) 
    559535      ! 
    560536      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav') 
    … …  
    562538   END SUBROUTINE sbc_isf_cav 
    563539 
    564  
    565    SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit ) 
     540   SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf ) 
    566541      !!---------------------------------------------------------------------- 
    567542      !! ** Purpose    : compute the coefficient echange for heat flux 
    … …  
    572547      !!                Jenkins et al., 2010, JPO, p2298-2312 
    573548      !!--------------------------------------------------------------------- 
    574       REAL(wp), INTENT(inout) :: gt, gs, zqhisf, zqwisf 
    575       INTEGER , INTENT(in)    :: ji,jj 
    576       LOGICAL , INTENT(inout) :: lit 
    577  
    578       INTEGER  :: ikt                 ! loop index 
    579       REAL(wp) :: zut, zvt, zustar           ! U, V at T point and friction velocity 
     549      REAL(wp), DIMENSION(:,:), INTENT(out) :: pgt, pgs 
     550      REAL(wp), DIMENSION(:,:), INTENT(in ) :: pqhisf, pqwisf 
     551      ! 
     552      INTEGER  :: ikt                         
     553      INTEGER  :: ji, jj                     ! loop index 
     554      REAL(wp), DIMENSION(:,:), POINTER :: zustar           ! U, V at T point and friction velocity 
    580555      REAL(wp) :: zdku, zdkv                 ! U, V shear  
    581556      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number  
    … …  
    587562      REAL(wp) :: zcoef                      ! temporary coef 
    588563      REAL(wp) :: zdep 
    589       REAL(wp), PARAMETER :: zxsiN = 0.052   ! dimensionless constant 
    590       REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number 
    591       REAL(wp), PARAMETER :: znu   = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1) 
    592       REAL(wp) ::   rcs      = 1.0e-3_wp        ! conversion: mm/s ==> m/s 
     564      REAL(wp) :: zeps = 1.0e-20_wp     
     565      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant 
     566      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1) 
    593567      REAL(wp), DIMENSION(2) :: zts, zab 
    594568      !!--------------------------------------------------------------------- 
    595       ! 
    596       IF( nn_gammablk == 0 ) THEN 
    597       !! gamma is constant (specified in namelist) 
    598          gt = rn_gammat0 
    599          gs = rn_gammas0 
    600          lit = .FALSE. 
    601       ELSE IF ( nn_gammablk == 1 ) THEN 
    602       !! gamma is assume to be proportional to u*  
    603       !! WARNING in case of Losh 2008 tbl parametrization,  
    604       !! you have to used the mean value of u in the boundary layer)  
    605       !! not yet coded 
    606       !! Jenkins et al., 2010, JPO, p2298-2312 
    607          ikt = mikt(ji,jj) 
    608       !! Compute U and V at T points 
    609    !      zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) ) 
    610    !      zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) ) 
    611           zut = utbl(ji,jj) 
    612           zvt = vtbl(ji,jj) 
    613  
    614       !! compute ustar 
    615          zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) ) 
    616       !! Compute mean value over the TBL 
    617  
    618       !! Compute gammats 
    619          gt = zustar * rn_gammat0 
    620          gs = zustar * rn_gammas0 
    621          lit = .FALSE. 
    622       ELSE IF ( nn_gammablk == 2 ) THEN 
    623       !! gamma depends of stability of boundary layer 
    624       !! WARNING in case of Losh 2008 tbl parametrization,  
    625       !! you have to used the mean value of u in the boundary layer)  
    626       !! not yet coded 
    627       !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 
    628       !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 
     569      CALL wrk_alloc( jpi,jpj, zustar ) 
     570      ! 
     571      SELECT CASE ( nn_gammablk ) 
     572      CASE ( 0 ) ! gamma is constant (specified in namelist) 
     573         !! ISOMIP formulation (Hunter et al, 2006) 
     574         pgt(:,:) = rn_gammat0 
     575         pgs(:,:) = rn_gammas0 
     576 
     577      CASE ( 1 ) ! gamma is assume to be proportional to u* 
     578         !! Jenkins et al., 2010, JPO, p2298-2312 
     579         !! Adopted by Asay-Davis et al. (2015) 
     580 
     581         !! compute ustar (eq. 24) 
     582         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 
     583 
     584         !! Compute gammats 
     585         pgt(:,:) = zustar(:,:) * rn_gammat0 
     586         pgs(:,:) = zustar(:,:) * rn_gammas0 
     587       
     588      CASE ( 2 ) ! gamma depends of stability of boundary layer 
     589         !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14 
     590         !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO) 
     591         !! compute ustar 
     592         zustar(:,:) = SQRT( rn_tfri2 * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + rn_tfeb2) ) 
     593 
     594         !! compute Pr and Sc number (can be improved) 
     595         zPr =   13.8_wp 
     596         zSc = 2432.0_wp 
     597 
     598         !! compute gamma mole 
     599         zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp 
     600         zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp 
     601 
     602         !! compute gamma 
     603         DO ji=2,jpi 
     604            DO jj=2,jpj 
    629605               ikt = mikt(ji,jj) 
    630606 
    631       !! Compute U and V at T points 
    632                zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) ) 
    633                zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) ) 
    634  
    635       !! compute ustar 
    636                zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) ) 
    637                IF (zustar == 0._wp) THEN           ! only for kt = 1 I think 
    638                  gt = rn_gammat0 
    639                  gs = rn_gammas0 
     607               IF (zustar(ji,jj) == 0._wp) THEN           ! only for kt = 1 I think 
     608                  pgt = rn_gammat0 
     609                  pgs = rn_gammas0 
    640610               ELSE 
    641       !! compute Rc number (as done in zdfric.F90) 
    642                zcoef = 0.5 / fse3w(ji,jj,ikt) 
    643                !                                            ! shear of horizontal velocity 
    644                zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )   & 
    645                   &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  ) 
    646                zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )   & 
    647                   &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  ) 
    648                !                                            ! richardson number (minimum value set to zero) 
    649                zRc = rn2(ji,jj,ikt+1) / ( zdku*zdku + zdkv*zdkv + 1.e-20 ) 
    650  
    651       !! compute Pr and Sc number (can be improved) 
    652                zPr =   13.8 
    653                zSc = 2432.0 
    654  
    655       !! compute gamma mole 
    656                zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0 
    657                zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0 
    658  
    659       !! compute bouyancy  
    660                zts(jp_tem) = ttbl(ji,jj) 
    661                zts(jp_sal) = stbl(ji,jj) 
    662                zdep        = fsdepw(ji,jj,ikt) 
    663                ! 
    664                CALL eos_rab( zts, zdep, zab ) 
     611                  !! compute Rc number (as done in zdfric.F90) 
     612                  zcoef = 0.5_wp / fse3w(ji,jj,ikt) 
     613                  !                                            ! shear of horizontal velocity 
     614                  zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  & 
     615                     &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  ) 
     616                  zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  & 
     617                     &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  ) 
     618                  !                                            ! richardson number (minimum value set to zero) 
     619                  zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps ) 
     620 
     621                  !! compute bouyancy  
     622                  zts(jp_tem) = ttbl(ji,jj) 
     623                  zts(jp_sal) = stbl(ji,jj) 
     624                  zdep        = fsdepw(ji,jj,ikt) 
    665625                  ! 
    666       !! compute length scale  
    667                zbuofdep = grav * ( zab(jp_tem) * zqhisf - zab(jp_sal) * zqwisf )  !!!!!!!!!!!!!!!!!!!!!!!!!!!! 
    668  
    669       !! compute Monin Obukov Length 
    670                ! Maximum boundary layer depth 
    671                zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) -0.001 
    672                ! Compute Monin obukhov length scale at the surface and Ekman depth: 
    673                zmob   = zustar ** 3 / (vkarmn * (zbuofdep + epsln)) 
    674                zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 
    675  
    676       !! compute eta* (stability parameter) 
    677                zetastar = 1 / ( SQRT(1 + MAX(zxsiN * zustar / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0))) 
    678  
    679       !! compute the sublayer thickness 
    680                zhnu = 5 * znu / zustar 
    681       !! compute gamma turb 
    682                zgturb = 1/vkarmn * LOG(zustar * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) & 
    683                &      + 1 / ( 2 * zxsiN * zetastar ) - 1/vkarmn 
    684  
    685       !! compute gammats 
    686                gt = zustar / (zgturb + zgmolet) 
    687                gs = zustar / (zgturb + zgmoles) 
     626                  CALL eos_rab( zts, zdep, zab ) 
     627                  ! 
     628                  !! compute length scale  
     629                  zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     630 
     631                  !! compute Monin Obukov Length 
     632                  ! Maximum boundary layer depth 
     633                  zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) - 0.001_wp 
     634                  ! Compute Monin obukhov length scale at the surface and Ekman depth: 
     635                  zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps)) 
     636                  zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt) 
     637 
     638                  !! compute eta* (stability parameter) 
     639                  zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0_wp))) 
     640 
     641                  !! compute the sublayer thickness 
     642                  zhnu = 5 * znu / zustar(ji,jj) 
     643 
     644                  !! compute gamma turb 
     645                  zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) & 
     646                  &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn 
     647 
     648                  !! compute gammats 
     649                  pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 
     650                  pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 
    688651               END IF 
    689       END IF 
    690       ! 
    691    END SUBROUTINE 
    692  
    693  
    694    SUBROUTINE sbc_isf_tbl( varin, varout, cptin ) 
     652            END DO 
     653         END DO 
     654         CALL lbc_lnk(pgt(:,:),'T',1.) 
     655         CALL lbc_lnk(pgs(:,:),'T',1.) 
     656      END SELECT 
     657      CALL wrk_dealloc( jpi,jpj, zustar ) 
     658      ! 
     659   END SUBROUTINE sbc_isf_gammats 
     660 
     661   SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin ) 
    695662      !!---------------------------------------------------------------------- 
    696663      !!                  ***  SUBROUTINE sbc_isf_tbl  *** 
    697664      !! 
    698       !! ** Purpose : compute mean T/S/U/V in the boundary layer  
    699       !! 
    700       !!---------------------------------------------------------------------- 
    701       REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin 
    702       REAL(wp), DIMENSION(:,:)  , INTENT(out):: varout 
     665      !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point 
     666      !! 
     667      !!---------------------------------------------------------------------- 
     668      REAL(wp), DIMENSION(:,:,:), INTENT( in  ) :: pvarin 
     669      REAL(wp), DIMENSION(:,:)  , INTENT( out ) :: pvarout 
     670      CHARACTER(len=1),           INTENT( in  ) :: cd_ptin ! point of variable in/out 
     671      ! 
     672      REAL(wp) :: ze3, zhk 
     673      REAL(wp), DIMENSION(:,:), POINTER :: zhisf_tbl ! thickness of the tbl 
     674 
     675      INTEGER :: ji, jj, jk                  ! loop index 
     676      INTEGER :: ikt, ikb                    ! top and bottom index of the tbl 
     677      !!---------------------------------------------------------------------- 
     678      ! allocation 
     679      CALL wrk_alloc( jpi,jpj, zhisf_tbl) 
    703680       
    704       CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out 
    705  
    706       REAL(wp) :: ze3, zhk 
    707       REAL(wp), DIMENSION(:,:), POINTER :: zikt 
    708  
    709       INTEGER :: ji,jj,jk 
    710       INTEGER :: ikt,ikb 
    711       INTEGER, DIMENSION(:,:), POINTER :: mkt, mkb 
    712  
    713       CALL wrk_alloc( jpi,jpj, mkt, mkb  ) 
    714       CALL wrk_alloc( jpi,jpj, zikt ) 
    715  
    716       ! get first and last level of tbl 
    717       mkt(:,:) = misfkt(:,:) 
    718       mkb(:,:) = misfkb(:,:) 
    719  
    720       varout(:,:)=0._wp 
    721       DO jj = 2,jpj 
    722          DO ji = 2,jpi 
    723             IF (ssmask(ji,jj) == 1) THEN 
    724                ikt = mkt(ji,jj) 
    725                ikb = mkb(ji,jj) 
     681      ! initialisation 
     682      pvarout(:,:)=0._wp 
     683    
     684      SELECT CASE ( cd_ptin ) 
     685      CASE ( 'U' ) ! compute U in the top boundary layer at T- point  
     686         DO jj = 1,jpj 
     687            DO ji = 1,jpi 
     688               ikt = miku(ji,jj) ; ikb = miku(ji,jj) 
     689               ! thickness of boundary layer at least the top level thickness 
     690               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3u_n(ji,jj,ikt)) 
     691 
     692               ! determine the deepest level influenced by the boundary layer 
     693               DO jk = ikt+1, mbku(ji,jj) 
     694                  IF ( (SUM(fse3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 
     695               END DO 
     696               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
     697 
     698               ! level fully include in the ice shelf boundary layer 
     699               DO jk = ikt, ikb - 1 
     700                  ze3 = fse3u_n(ji,jj,jk) 
     701                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 
     702               END DO 
     703 
     704               ! level partially include in ice shelf boundary layer  
     705               zhk = SUM( fse3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 
     706               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
     707            END DO 
     708         END DO 
     709         DO jj = 2,jpj 
     710            DO ji = 2,jpi 
     711               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj)) 
     712            END DO 
     713         END DO 
     714         CALL lbc_lnk(pvarout,'T',-1.) 
     715       
     716      CASE ( 'V' ) ! compute V in the top boundary layer at T- point  
     717         DO jj = 1,jpj 
     718            DO ji = 1,jpi 
     719               ikt = mikv(ji,jj) ; ikb = mikv(ji,jj) 
     720               ! thickness of boundary layer at least the top level thickness 
     721               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3v_n(ji,jj,ikt)) 
     722 
     723               ! determine the deepest level influenced by the boundary layer 
     724               DO jk = ikt+1, mbkv(ji,jj) 
     725                  IF ( (SUM(fse3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 
     726               END DO 
     727               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(fse3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
     728 
     729               ! level fully include in the ice shelf boundary layer 
     730               DO jk = ikt, ikb - 1 
     731                  ze3 = fse3v_n(ji,jj,jk) 
     732                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 
     733               END DO 
     734 
     735               ! level partially include in ice shelf boundary layer  
     736               zhk = SUM( fse3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj) 
     737               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
     738            END DO 
     739         END DO 
     740         DO jj = 2,jpj 
     741            DO ji = 2,jpi 
     742               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1)) 
     743            END DO 
     744         END DO 
     745         CALL lbc_lnk(pvarout,'T',-1.) 
     746 
     747      CASE ( 'T' ) ! compute T in the top boundary layer at T- point  
     748         DO jj = 1,jpj 
     749            DO ji = 1,jpi 
     750               ikt = misfkt(ji,jj) 
     751               ikb = misfkb(ji,jj) 
    726752 
    727753               ! level fully include in the ice shelf boundary layer 
    728754               DO jk = ikt, ikb - 1 
    729755                  ze3 = fse3t_n(ji,jj,jk) 
    730                   IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 
    731                   IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) & 
    732                      &                                                       * r1_hisf_tbl(ji,jj) * ze3 
    733                   IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) & 
    734                      &                                                       * r1_hisf_tbl(ji,jj) * ze3 
     756                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 
    735757               END DO 
    736758 
    737759               ! level partially include in ice shelf boundary layer  
    738760               zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) 
    739                IF (cptin == 'T') & 
    740                    &  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk) 
    741                IF (cptin == 'U') & 
    742                    &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk) 
    743                IF (cptin == 'V') & 
    744                    &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk) 
    745             END IF 
    746          END DO 
    747       END DO 
    748  
    749       CALL wrk_dealloc( jpi,jpj, mkt, mkb )       
    750       CALL wrk_dealloc( jpi,jpj, zikt )  
    751  
    752       IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.) 
    753       IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.) 
     761               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 
     762            END DO 
     763         END DO 
     764      END SELECT 
     765 
     766      ! mask mean tbl value 
     767      pvarout(:,:) = pvarout(:,:) * ssmask(:,:) 
     768 
     769      ! deallocation 
     770      CALL wrk_dealloc( jpi,jpj, zhisf_tbl )       
    754771      ! 
    755772   END SUBROUTINE sbc_isf_tbl 
    … …  
    768785      !! ** Action  :   phdivn   decreased by the runoff inflow 
    769786      !!---------------------------------------------------------------------- 
    770       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence 
    771       !! 
     787      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdivn   ! horizontal divergence 
     788      !  
    772789      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    773790      INTEGER  ::   ikt, ikb  
    774       INTEGER  ::   nk_isf 
    775       REAL(wp)     ::   zhk, z1_hisf_tbl, zhisf_tbl 
    776       REAL(wp)     ::   zfact     ! local scalar 
     791      REAL(wp) ::   zhk 
     792      REAL(wp) ::   zfact     ! local scalar 
    777793      !!---------------------------------------------------------------------- 
    778794      ! 
    779795      zfact   = 0.5_wp 
    780796      ! 
    781       IF (lk_vvl) THEN     ! need to re compute level distribution of isf fresh water 
     797      IF ( lk_vvl ) THEN     ! need to re compute level distribution of isf fresh water 
    782798         DO jj = 1,jpj 
    783799            DO ji = 1,jpi 
    … …  
    788804 
    789805               ! determine the deepest level influenced by the boundary layer 
    790                ! test on tmask useless ????? 
    791                DO jk = ikt, mbkt(ji,jj) 
    792                   IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
     806               DO jk = ikt+1, mbkt(ji,jj) 
     807                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) <  rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
    793808               END DO 
    794809               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
    … …  
    796811               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 
    797812 
    798                zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1 
     813               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 
    799814               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer 
    800815            END DO 
    801816         END DO 
    802       END IF ! vvl case 
    803       ! 
     817      END IF  
     818      ! 
     819      !==   ice shelf melting distributed over several levels   ==! 
    804820      DO jj = 1,jpj 
    805821         DO ji = 1,jpi 
    … …  
    809825               DO jk = ikt, ikb - 1 
    810826                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) & 
    811                     &               * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 
     827                    &              * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact 
    812828               END DO 
    813829               ! level partially include in ice shelf boundary layer  
    814830               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) & 
    815                   &             + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj)  
    816             !==   ice shelf melting mass distributed over several levels   ==! 
     831                    &            + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj)  
    817832         END DO 
    818833      END DO 
    819834      ! 
    820835   END SUBROUTINE sbc_isf_div 
    821  
    822  
    823    FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti ) 
    824       !!---------------------------------------------------------------------- 
    825       !!                 ***  ROUTINE eos_init  *** 
    826       !! 
    827       !! ** Purpose :   Compute the in-situ temperature [Celcius] 
    828       !! 
    829       !! ** Method  :    
    830       !! 
    831       !! Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408 
    832       !!---------------------------------------------------------------------- 
    833       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptem   ! potential temperature [Celcius] 
    834       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu] 
    835       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ppress ! pressure             [dBar] 
    836       REAL(wp), DIMENSION(:,:), POINTER           ::   pti    ! in-situ temperature [Celcius] 
    837 !      REAL(wp) :: fsatg 
    838 !      REAL(wp) :: pfps, pfpt, pfphp  
    839       REAL(wp) :: zt, zs, zp, zh, zq, zxk 
    840       INTEGER  :: ji, jj            ! dummy loop indices 
    841       ! 
    842       CALL wrk_alloc( jpi,jpj, pti  ) 
    843       !  
    844       DO jj=1,jpj 
    845          DO ji=1,jpi 
    846             zh = ppress(ji,jj) 
    847 ! Theta1 
    848             zt = ptem(ji,jj) 
    849             zs = psal(ji,jj) 
    850             zp = 0.0 
    851             zxk= zh * fsatg( zs, zt, zp ) 
    852             zt = zt + 0.5 * zxk 
    853             zq = zxk 
    854 ! Theta2 
    855             zp = zp + 0.5 * zh 
    856             zxk= zh*fsatg( zs, zt, zp ) 
    857             zt = zt + 0.29289322 * ( zxk - zq ) 
    858             zq = 0.58578644 * zxk + 0.121320344 * zq 
    859 ! Theta3 
    860             zxk= zh * fsatg( zs, zt, zp ) 
    861             zt = zt + 1.707106781 * ( zxk - zq ) 
    862             zq = 3.414213562 * zxk - 4.121320344 * zq 
    863 ! Theta4 
    864             zp = zp + 0.5 * zh 
    865             zxk= zh * fsatg( zs, zt, zp ) 
    866             pti(ji,jj) = zt + ( zxk - 2.0 * zq ) / 6.0 
    867          END DO 
    868       END DO 
    869       ! 
    870       CALL wrk_dealloc( jpi,jpj, pti  ) 
    871       ! 
    872    END FUNCTION tinsitu 
    873  
    874  
    875    FUNCTION fsatg( pfps, pfpt, pfphp ) 
    876       !!---------------------------------------------------------------------- 
    877       !!                 ***  FUNCTION fsatg  *** 
    878       !! 
    879       !! ** Purpose    :   Compute the Adiabatic laspse rate [Celcius].[decibar]^-1 
    880       !! 
    881       !! ** Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408 
    882       !!  
    883       !! ** units      :   pressure        pfphp    decibars 
    884       !!                   temperature     pfpt     deg celsius (ipts-68) 
    885       !!                   salinity        pfps     (ipss-78) 
    886       !!                   adiabatic       fsatg    deg. c/decibar 
    887       !!---------------------------------------------------------------------- 
    888       REAL(wp) :: pfps, pfpt, pfphp  
    889       REAL(wp) :: fsatg 
    890       ! 
    891       fsatg = (((-2.1687e-16*pfpt+1.8676e-14)*pfpt-4.6206e-13)*pfphp         & 
    892         &    +((2.7759e-12*pfpt-1.1351e-10)*(pfps-35.)+((-5.4481e-14*pfpt    & 
    893         &    +8.733e-12)*pfpt-6.7795e-10)*pfpt+1.8741e-8))*pfphp             & 
    894         &    +(-4.2393e-8*pfpt+1.8932e-6)*(pfps-35.)                         & 
    895         &    +((6.6228e-10*pfpt-6.836e-8)*pfpt+8.5258e-6)*pfpt+3.5803e-5 
    896       ! 
    897     END FUNCTION fsatg 
    898     !!====================================================================== 
     836   !!====================================================================== 
    899837END MODULE sbcisf 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r6005 r6012  
    9090      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx, ln_blk_clio, ln_blk_core, ln_mixcpl,   & 
    9191         &             ln_blk_mfs, ln_apr_dyn, nn_ice, nn_ice_embd, ln_dm2dc   , ln_rnf   ,   & 
    92          &             ln_ssr    , nn_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   & 
     92         &             ln_ssr    , ln_isf    , nn_fwb, ln_cdgw    , ln_wave    , ln_sdw   ,   & 
    9393         &             nn_lsm    , nn_limflx , nn_components, ln_cpl 
    9494      INTEGER  ::   ios 
    … …  
    143143         WRITE(numout,*) '              daily mean to diurnal cycle qsr            ln_dm2dc    = ', ln_dm2dc  
    144144         WRITE(numout,*) '              runoff / runoff mouths                     ln_rnf      = ', ln_rnf 
    145          WRITE(numout,*) '              iceshelf formulation                       nn_isf      = ', nn_isf 
     145         WRITE(numout,*) '              iceshelf melting                           ln_isf      = ', ln_isf 
    146146         WRITE(numout,*) '              Sea Surface Restoring on SST and/or SSS    ln_ssr      = ', ln_ssr 
    147147         WRITE(numout,*) '              FreshWater Budget control  (=0/1/2)        nn_fwb      = ', nn_fwb 
    … …  
    181181 
    182182      !                          ! Checks: 
    183       IF( nn_isf .EQ. 0 ) THEN                      ! variable initialisation if no ice shelf  
    184          IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'sbc_init : unable to allocate sbc_isf arrays' ) 
     183      IF( .NOT. ln_isf ) THEN                      ! variable initialisation if no ice shelf  
     184         IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' ) 
    185185         fwfisf  (:,:)   = 0.0_wp ; fwfisf_b  (:,:)   = 0.0_wp 
    186186         risf_tsc(:,:,:) = 0.0_wp ; risf_tsc_b(:,:,:) = 0.0_wp 
    187          rdivisf       = 0.0_wp 
    188187      END IF 
     188 
    189189      IF( nn_ice == 0 .AND. nn_components /= jp_iam_opa )   fr_i(:,:) = 0.e0 ! no ice in the domain, ice fraction is always zero 
    190190 
    … …  
    378378      IF( ln_icebergs    )   CALL icb_stp( kt )                   ! compute icebergs 
    379379 
    380       IF( nn_isf   /= 0  )   CALL sbc_isf( kt )                    ! compute iceshelves 
     380      IF( ln_isf         )   CALL sbc_isf( kt )                   ! compute iceshelves 
    381381 
    382382      IF( ln_rnf         )   CALL sbc_rnf( kt )                   ! add runoffs to fresh water fluxes 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r5836 r6012  
    142142         IF( ln_rnf_tem ) THEN                                       ! use runoffs temperature data 
    143143            rnf_tsc(:,:,jp_tem) = ( sf_t_rnf(1)%fnow(:,:,1) ) * rnf(:,:) * r1_rau0 
     144            CALL eos_fzp( sss_m(:,:), ztfrz(:,:) ) 
    144145            WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp )             ! if missing data value use SST as runoffs temperature 
    145146               rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 
    146147            END WHERE 
    147148            WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )             ! where fwf comes from melting of ice shelves or iceberg 
    148                ztfrz(:,:) = -1.9 !tfreez( sss_m(:,:) ) !PM to be discuss (trouble if sensitivity study) 
    149                rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * lfusisf * r1_rau0_rcp 
     149               rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * rlfusisf * r1_rau0_rcp 
    150150            END WHERE 
    151151         ELSE                                                        ! use SST as runoffs temperature 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90

    r5930 r6012  
    243243      ENDIF 
    244244      IF( ln_isfcav ) THEN                                                       ! ice-shelf cavities 
    245          IF(  ln_traadv_cen .AND. nn_cen_v /= 4    .OR.   &                            ! NO 4th order with ISF 
    246             & ln_traadv_fct .AND. nn_fct_v /= 4   )   CALL ctl_stop( 'tra_adv_init: 4th order COMPACT scheme not allowed with ISF' ) 
     245         IF(  ln_traadv_cen .AND. nn_cen_v == 4    .OR.   &                            ! NO 4th order with ISF 
     246            & ln_traadv_fct .AND. nn_fct_v == 4   )   CALL ctl_stop( 'tra_adv_init: 4th order COMPACT scheme not allowed with ISF' ) 
    247247      ENDIF 
    248248      ! 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90

    r5836 r6012  
    103103      ! 
    104104      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices 
     105      INTEGER  ::  ikt 
    105106      INTEGER  ::  ierr             ! local integer 
    106107      REAL(wp) ::  zmsku, zahu_w, zabe1, zcof1, zcoef3   ! local scalars 
    … …  
    228229            DO jj = 1, jpjm1              ! bottom correction (partial bottom cell) 
    229230               DO ji = 1, fs_jpim1   ! vector opt. 
    230 !!gm  the following anonymous remark is to considered: ! IF useless if zpshde defines pgu everywhere 
    231231                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
    232232                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
    … …  
    255255            ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk) 
    256256            ENDIF 
    257 !!gm I don't understand why we should need this.... since wmask is used instead of tmask 
    258 !         IF ( ln_isfcav ) THEN 
    259 !            DO jj = 1, jpj 
    260 !               DO ji = 1, jpi   ! vector opt. 
    261 !                  ikt = mikt(ji,jj) ! surface level 
    262 !                  zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn  ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1) 
    263 !                  zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt) 
    264 !               END DO 
    265 !            END DO 
    266 !         END IF 
    267 !!gm  
    268257 
    269258            DO jj = 1 , jpjm1            !==  Horizontal fluxes 
    … …  
    272261                  zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk) 
    273262                  ! 
    274                   zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   & 
    275                      &             + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk  ), 1. ) 
    276                   ! 
    277                   zmskv = 1. / MAX(  tmask(ji,jj+1,jk  ) + tmask(ji,jj,jk+1)   & 
    278                      &             + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk  ), 1. ) 
     263                  zmsku = 1. / MAX(  wmask(ji+1,jj,jk  ) + wmask(ji,jj,jk+1)   & 
     264                     &             + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk  ), 1. ) 
     265                  ! 
     266                  zmskv = 1. / MAX(  wmask(ji,jj+1,jk  ) + wmask(ji,jj,jk+1)   & 
     267                     &             + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk  ), 1. ) 
    279268                  ! 
    280269                  zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 
    … …  
    317306               DO ji = fs_2, fs_jpim1   ! vector opt. 
    318307                  ! 
    319                   zmsku = tmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
     308                  zmsku = wmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          & 
    320309                     &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  ) 
    321                   zmskv = tmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
     310                  zmskv = wmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          & 
    322311                     &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  ) 
    323312                     ! 
    … …  
    343332               DO jj = 1, jpjm1 
    344333                  DO ji = fs_2, fs_jpim1 
    345                      ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)   & 
     334                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * wmask(ji,jj,jk)   & 
    346335                        &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             & 
    347336                        &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 
    … …  
    358347                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    & 
    359348                           &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   & 
    360                            &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) / fse3w(ji,jj,jk) 
     349                           &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * wmask(ji,jj,jk) / fse3w(ji,jj,jk) 
    361350                     END DO 
    362351                  END DO 
    … …  
    366355                  DO jj = 1, jpjm1 
    367356                     DO ji = fs_2, fs_jpim1 
    368                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)                      & 
     357                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * wmask(ji,jj,jk)                      & 
    369358                           &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   & 
    370359                           &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   ) 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90

    r5930 r6012  
    2828   USE sbc_oce         ! surface boundary condition: ocean 
    2929   USE sbcrnf          ! river runoffs 
    30    USE sbcisf          ! ice shelf melting/freezing 
     30   USE sbcisf          ! ice shelf melting 
    3131   USE zdf_oce         ! ocean vertical mixing 
    3232   USE domvvl          ! variable volume 
    … …  
    262262         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration 
    263263         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs 
    264          IF (nn_isf .GE. 1) THEN  
    265             ll_isf = .TRUE.            ! active  tracers case  and  ice shelf melting/freezing 
    266          ELSE 
    267             ll_isf = .FALSE. 
    268          END IF 
     264         ll_isf     = ln_isf           ! active  tracers case  and  ice shelf melting 
    269265      ELSE                           
    270266         ll_traqsr  = .FALSE.          ! active  tracers case and NO solar penetration 
    271267         ll_rnf     = .FALSE.          ! passive tracers or NO river runoffs 
    272          ll_isf     = .FALSE.          ! passive tracers or NO ice shelf melting/freezing 
     268         ll_isf     = .FALSE.          ! passive tracers or NO ice shelf melting 
    273269      ENDIF 
    274270      ! 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r6006 r6012  
    119119      INTEGER  ::   ikt, ikb  
    120120      REAL(wp) ::   zfact, z1_e3t, zdep 
    121       REAL(wp) ::   zalpha, zhk 
     121      REAL(wp) ::   zt_frz, zpress 
    122122      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds 
    123123      !!---------------------------------------------------------------------- 
    … …  
    218218      !---------------------------------------- 
    219219      ! 
    220       IF( nn_isf > 0 ) THEN 
     220      IF( ln_isf ) THEN 
    221221         zfact = 0.5_wp 
    222222         DO jj = 2, jpj 
    … …  
    227227    
    228228               ! level fully include in the ice shelf boundary layer 
    229                ! if isfdiv, we have to remove heat flux due to inflow at 0oC (as in rnf when you add rnf at sst) 
    230229               ! sign - because fwf sign of evapo (rnf sign of precip) 
    231230               DO jk = ikt, ikb - 1 
    232231               ! compute trend 
    233                   tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                          & 
    234                      &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) 
    235                   tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                          & 
    236                      &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) 
     232                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                                & 
     233                     &           + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             & 
     234                     &           * r1_hisf_tbl(ji,jj) 
    237235               END DO 
    238236    
    239237               ! level partially include in ice shelf boundary layer  
    240238               ! compute trend 
    241                tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                           & 
    242                   &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
    243                tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                           & 
    244                   &              + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)  
     239               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                                 & 
     240                  &              + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             & 
     241                  &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) 
     242 
    245243            END DO 
    246244         END DO 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r5836 r6012  
    191191      ! 
    192192   END SUBROUTINE zps_hde 
    193  
    194  
    195    SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu , pgtv , pgtui, pgtvi,                                   & 
    196       &                              prd, pgru , pgrv , pmru , pmrv , pgzu , pgzv , pge3ru , pge3rv ,   & 
    197       &                                   pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi ) 
    198       !!---------------------------------------------------------------------- 
    199       !!                     ***  ROUTINE zps_hde  *** 
     193   ! 
     194   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, pgtui, pgtvi,  & 
     195      &                          prd, pgru, pgrv, pgrui, pgrvi ) 
     196      !!---------------------------------------------------------------------- 
     197      !!                     ***  ROUTINE zps_hde_isf  *** 
    200198      !!                     
    201199      !! ** Purpose :   Compute the horizontal derivative of T, S and rho 
    202200      !!      at u- and v-points with a linear interpolation for z-coordinate 
    203       !!      with partial steps. 
     201      !!      with partial steps for top (ice shelf) and bottom. 
    204202      !! 
    205203      !! ** Method  :   In z-coord with partial steps, scale factors on last  
    206204      !!      levels are different for each grid point, so that T, S and rd  
    207205      !!      points are not at the same depth as in z-coord. To have horizontal 
    208       !!      gradients again, we interpolate T and S at the good depth :  
     206      !!      gradients again, we interpolate T and S at the good depth : 
     207      !!      For the bottom case: 
    209208      !!      Linear interpolation of T, S    
    210209      !!         Computation of di(tb) and dj(tb) by vertical interpolation: 
    … …  
    235234      !!          di(rho) = rd~ - rd(i,j,k)   or   rd(i+1,j,k) - rd~ 
    236235      !! 
     236      !!      For the top case (ice shelf): As for the bottom case but upside down 
     237      !! 
    237238      !! ** Action  : compute for top and bottom interfaces 
    238239      !!              - pgtu, pgtv, pgtui, pgtvi: horizontal gradient of tracer at u- & v-points 
    239240      !!              - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points 
    240       !!              - pmru, pmrv, pmrui, pmrvi: horizontal sum of rho at u- & v- point (used in dynhpg with vvl) 
    241       !!              - pgzu, pgzv, pgzui, pgzvi: horizontal gradient of z at u- and v- point (used in dynhpg with vvl) 
    242       !!              - pge3ru, pge3rv, pge3rui, pge3rvi: horizontal gradient of rho weighted by local e3w at u- & v-points  
    243       !!---------------------------------------------------------------------- 
    244       INTEGER                              , INTENT(in   )           ::  kt                ! ocean time-step index 
    245       INTEGER                              , INTENT(in   )           ::  kjpt              ! number of tracers 
    246       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta               ! 4D tracers fields 
    247       !                                                              !!  u-point ! v-point ! 
    248       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu    , pgtv    ! bottom GRADh( ptra )   
    249       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui   , pgtvi   ! top    GRADh( ptra ) 
    250       REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd               ! 3D density anomaly fields 
    251       !                                                              !!  u-point ! v-point ! 
    252       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru    , pgrv    ! bottom GRADh( prd  ) 
    253       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmru    , pmrv    ! bottom SUM  ( prd  ) 
    254       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu    , pgzv    ! bottom GRADh( z    )  
    255       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru  , pge3rv  ! bottom GRADh( prd  ) weighted by e3w 
    256       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui   , pgrvi   ! top    GRADh( prd  )  
    257       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmrui   , pmrvi   ! top    SUM  ( prd  )  
    258       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzui   , pgzvi   ! top    GRADh( z    )  
    259       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3rui , pge3rvi ! top    GRADh( prd  ) weighted by e3w 
     241      !!---------------------------------------------------------------------- 
     242      INTEGER                              , INTENT(in   )           ::  kt           ! ocean time-step index 
     243      INTEGER                              , INTENT(in   )           ::  kjpt         ! number of tracers 
     244      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta          ! 4D tracers fields 
     245      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv   ! hor. grad. of ptra at u- & v-pts  
     246      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF) 
     247      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd          ! 3D density anomaly fields 
     248      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv   ! hor. grad of prd at u- & v-pts (bottom) 
     249      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top) 
    260250      ! 
    261251      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
    262252      INTEGER  ::   iku, ikv, ikum1, ikvm1,ikup1, ikvp1   ! partial step level (ocean bottom level) at u- and v-points 
    263       REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv, zdzwu, zdzwv, zdzwuip1, zdzwvjp1  ! temporary scalars 
     253      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv             ! temporary scalars 
    264254      REAL(wp), DIMENSION(jpi,jpj)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos 
    265255      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::  zti, ztj             !  
    … …  
    277267         DO jj = 1, jpjm1 
    278268            DO ji = 1, jpim1 
    279                iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    280                ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
     269 
     270               iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
     271               ikv = mbkv(ji,jj); ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
     272               ze3wu = fsdept_n(ji+1,jj,iku) - fsdept_n(ji,jj,iku) 
     273               ze3wv = fsdept_n(ji,jj+1,ikv) - fsdept_n(ji,jj,ikv) 
     274               ! 
     275               ! i- direction 
     276               IF( ze3wu >= 0._wp ) THEN      ! case 1 
     277                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku) 
     278                  ! interpolated values of tracers 
     279                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
     280                  ! gradient of  tracers 
     281                  pgtu(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     282               ELSE                           ! case 2 
     283                  zmaxu = -ze3wu / fse3w(ji,jj,iku) 
     284                  ! interpolated values of tracers 
     285                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
     286                  ! gradient of tracers 
     287                  pgtu(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     288               ENDIF 
     289               ! 
     290               ! j- direction 
     291               IF( ze3wv >= 0._wp ) THEN      ! case 1 
     292                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv) 
     293                  ! interpolated values of tracers 
     294                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
     295                  ! gradient of tracers 
     296                  pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     297               ELSE                           ! case 2 
     298                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv) 
     299                  ! interpolated values of tracers 
     300                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
     301                  ! gradient of tracers 
     302                  pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     303               ENDIF 
     304 
     305            END DO 
     306         END DO 
     307         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     308         ! 
     309      END DO 
     310 
     311      ! horizontal derivative of density anomalies (rd) 
     312      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
     313         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ;  
     314         ! 
     315         DO jj = 1, jpjm1 
     316            DO ji = 1, jpim1 
     317 
     318               iku = mbku(ji,jj) 
     319               ikv = mbkv(ji,jj) 
     320               ze3wu = fsdept_n(ji+1,jj,iku) - fsdept_n(ji,jj,iku) 
     321               ze3wv = fsdept_n(ji,jj+1,ikv) - fsdept_n(ji,jj,ikv) 
     322               ! 
     323               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)    ! i-direction: case 1 
     324               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)    ! -     -      case 2 
     325               ENDIF 
     326               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)    ! j-direction: case 1 
     327               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)    ! -     -      case 2 
     328               ENDIF 
     329 
     330            END DO 
     331         END DO 
     332 
     333         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
     334         ! step and store it in  zri, zrj for each  case 
     335         CALL eos( zti, zhi, zri ) 
     336         CALL eos( ztj, zhj, zrj ) 
     337 
     338         DO jj = 1, jpjm1                 ! Gradient of density at the last level  
     339            DO ji = 1, jpim1 
     340 
     341               iku = mbku(ji,jj) 
     342               ikv = mbkv(ji,jj) 
     343               ze3wu = fsdept_n(ji+1,jj,iku) - fsdept_n(ji,jj,iku) 
     344               ze3wv = fsdept_n(ji,jj+1,ikv) - fsdept_n(ji,jj,ikv) 
     345 
     346               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1 
     347               ELSE                        ;   pgru(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2 
     348               ENDIF 
     349               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = ssvmask(ji,jj) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1 
     350               ELSE                        ;   pgrv(ji,jj) = ssvmask(ji,jj) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2 
     351               ENDIF 
     352 
     353            END DO 
     354         END DO 
     355 
     356         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions 
     357         ! 
     358      END IF 
     359      ! 
     360      !     !==  (ISH)  compute grui and gruvi  ==! 
     361      ! 
     362      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            ! 
     363         DO jj = 1, jpjm1 
     364            DO ji = 1, jpim1 
     365               iku = miku(ji,jj); ikup1 = miku(ji,jj) + 1 
     366               ikv = mikv(ji,jj); ikvp1 = mikv(ji,jj) + 1 
     367               ! 
    281368               ! (ISF) case partial step top and bottom in adjacent cell in vertical 
    282369               ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 
    283370               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 
    284371               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 
    285                ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 
    286                ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 
    287                ! 
     372               ze3wu  =  fsdept_n(ji,jj,iku) - fsdept_n(ji+1,jj,iku) 
     373               ze3wv  =  fsdept_n(ji,jj,ikv) - fsdept_n(ji,jj+1,ikv)  
     374 
    288375               ! i- direction 
    289376               IF( ze3wu >= 0._wp ) THEN      ! case 1 
    290                   zmaxu =  ze3wu / fse3w(ji+1,jj,iku) 
    291                   ! interpolated values of tracers 
    292                   zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
     377                  zmaxu = ze3wu / fse3w(ji+1,jj,ikup1) 
     378                  ! interpolated values of tracers 
     379                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikup1,jn) - pta(ji+1,jj,iku,jn) ) 
     380                  ! gradient of tracers 
     381                  pgtui(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     382               ELSE                           ! case 2 
     383                  zmaxu = - ze3wu / fse3w(ji,jj,ikup1) 
     384                  ! interpolated values of tracers 
     385                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikup1,jn) - pta(ji,jj,iku,jn) ) 
    293386                  ! gradient of  tracers 
    294                   pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
    295                ELSE                           ! case 2 
    296                   zmaxu = -ze3wu / fse3w(ji,jj,iku) 
    297                   ! interpolated values of tracers 
    298                   zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
    299                   ! gradient of tracers 
    300                   pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     387                  pgtui(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
    301388               ENDIF 
    302389               ! 
    303390               ! j- direction 
    304391               IF( ze3wv >= 0._wp ) THEN      ! case 1 
    305                   zmaxv =  ze3wv / fse3w(ji,jj+1,ikv) 
    306                   ! interpolated values of tracers 
    307                   ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
    308                   ! gradient of tracers 
    309                   pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
    310                ELSE                           ! case 2 
    311                   zmaxv =  -ze3wv / fse3w(ji,jj,ikv) 
    312                   ! interpolated values of tracers 
    313                   ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
    314                   ! gradient of tracers 
    315                   pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
    316                ENDIF 
    317             END DO 
    318          END DO 
    319          CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     392                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikvp1) 
     393                  ! interpolated values of tracers 
     394                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvp1,jn) - pta(ji,jj+1,ikv,jn) ) 
     395                  ! gradient of tracers 
     396                  pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     397               ELSE                           ! case 2 
     398                  zmaxv =  - ze3wv / fse3w(ji,jj,ikvp1) 
     399                  ! interpolated values of tracers 
     400                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvp1,jn) - pta(ji,jj,ikv,jn) ) 
     401                  ! gradient of tracers 
     402                  pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     403               ENDIF 
     404 
     405            END DO 
     406         END DO 
     407         CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. ); CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
    320408         ! 
    321409      END DO 
    … …  
    323411      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part) 
    324412         ! 
    325          pgru  (:,:)=0._wp   ;   pgrv  (:,:) = 0._wp 
    326          pgzu  (:,:)=0._wp   ;   pgzv  (:,:) = 0._wp  
    327          pmru  (:,:)=0._wp   ;   pmru  (:,:) = 0._wp  
    328          pge3ru(:,:)=0._wp   ;   pge3rv(:,:) = 0._wp  
    329          ! 
    330          DO jj = 1, jpjm1                 ! depth of the partial step level 
    331             DO ji = 1, jpim1 
    332                iku = mbku(ji,jj) 
    333                ikv = mbkv(ji,jj) 
    334                ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 
    335                ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 
    336                ! 
    337                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu     ! i-direction: case 1 
    338                ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) + ze3wu    ! -     -      case 2 
    339                ENDIF 
    340                IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) - ze3wv    ! j-direction: case 1 
    341                ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) + ze3wv    ! -     -      case 2 
    342                ENDIF 
     413         pgrui(:,:)  =0.0_wp; pgrvi(:,:)  =0.0_wp; 
     414         DO jj = 1, jpjm1 
     415            DO ji = 1, jpim1 
     416 
     417               iku = miku(ji,jj) 
     418               ikv = mikv(ji,jj) 
     419               ze3wu  =  fsdept_n(ji,jj,iku) - fsdept_n(ji+1,jj,iku) 
     420               ze3wv  =  fsdept_n(ji,jj,ikv) - fsdept_n(ji,jj+1,ikv)  
     421               ! 
     422               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)    ! i-direction: case 1 
     423               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)    ! -     -      case 2 
     424               ENDIF 
     425 
     426               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)    ! j-direction: case 1 
     427               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)    ! -     -      case 2 
     428               ENDIF 
     429 
    343430            END DO 
    344431         END DO 
    … …  
    346433         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj  
    347434         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj  
    348  
     435         ! 
    349436         DO jj = 1, jpjm1                 ! Gradient of density at the last level  
    350437            DO ji = 1, jpim1 
    351                iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    352                ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    353                ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku)) 
    354                ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv)) 
    355                IF( ze3wu >= 0._wp ) THEN  
    356                   pgzu(ji,jj) = (fsde3w(ji+1,jj,iku) - ze3wu) - fsde3w(ji,jj,iku) 
    357                   pgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1 
    358                   pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) + prd(ji,jj,iku) )   ! i: 1  
    359                   pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  & 
    360                                 * ( (fse3w(ji+1,jj,iku) - ze3wu )* ( zri(ji  ,jj    ) + prd(ji+1,jj,ikum1) + 2._wp) & 
    361                                    - fse3w(ji  ,jj,iku)          * ( prd(ji  ,jj,iku) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2 
    362                ELSE   
    363                   pgzu(ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) + ze3wu) 
    364                   pgru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2 
    365                   pmru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) )   ! i: 2 
    366                   pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  & 
    367                                 * (  fse3w(ji+1,jj,iku)          * ( prd(ji+1,jj,iku) + prd(ji+1,jj,ikum1) + 2._wp) & 
    368                                    -(fse3w(ji  ,jj,iku) + ze3wu) * ( zri(ji  ,jj    ) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2 
    369                ENDIF 
    370                IF( ze3wv >= 0._wp ) THEN 
    371                   pgzv(ji,jj) = (fsde3w(ji,jj+1,ikv) - ze3wv) - fsde3w(ji,jj,ikv)  
    372                   pgrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )   ! j: 1 
    373                   pmrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )   ! j: 1 
    374                   pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  & 
    375                                 * ( (fse3w(ji,jj+1,ikv) - ze3wv )* ( zrj(ji,jj      ) + prd(ji,jj+1,ikvm1) + 2._wp) & 
    376                                    - fse3w(ji,jj  ,ikv)          * ( prd(ji,jj  ,ikv) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2 
    377                ELSE  
    378                   pgzv(ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) + ze3wv) 
    379                   pgrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2 
    380                   pmrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )   ! j: 2 
    381                   pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  & 
    382                                 * (  fse3w(ji,jj+1,ikv)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikvm1) + 2._wp) & 
    383                                    -(fse3w(ji,jj  ,ikv) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2 
    384                ENDIF 
    385             END DO 
    386          END DO 
    387          CALL lbc_lnk( pgru   , 'U', -1. )   ;   CALL lbc_lnk( pgrv   , 'V', -1. )   ! Lateral boundary conditions 
    388          CALL lbc_lnk( pmru   , 'U',  1. )   ;   CALL lbc_lnk( pmrv   , 'V',  1. )   ! Lateral boundary conditions 
    389          CALL lbc_lnk( pgzu   , 'U', -1. )   ;   CALL lbc_lnk( pgzv   , 'V', -1. )   ! Lateral boundary conditions 
    390          CALL lbc_lnk( pge3ru , 'U', -1. )   ;   CALL lbc_lnk( pge3rv , 'V', -1. )   ! Lateral boundary conditions 
    391          ! 
    392       END IF 
    393       ! 
    394       !     !==  (ISH)  compute grui and gruvi  ==! 
    395       ! 
    396       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            ! 
    397          DO jj = 1, jpjm1 
    398             DO ji = 1, jpim1 
    399                iku = miku(ji,jj)   ;  ikup1 = miku(ji,jj) + 1 
    400                ikv = mikv(ji,jj)   ;  ikvp1 = mikv(ji,jj) + 1 
    401                ! 
    402                ! (ISF) case partial step top and bottom in adjacent cell in vertical 
    403                ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 
    404                ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 
    405                ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 
    406                ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))  
    407                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)) 
    408                ! i- direction 
    409                IF( ze3wu >= 0._wp ) THEN      ! case 1 
    410                   zmaxu = ze3wu / fse3w(ji+1,jj,iku+1) 
    411                   ! interpolated values of tracers 
    412                   zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) ) 
    413                   ! gradient of tracers 
    414                   pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
    415                ELSE                           ! case 2 
    416                   zmaxu = - ze3wu / fse3w(ji,jj,iku+1) 
    417                   ! interpolated values of tracers 
    418                   zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) ) 
    419                   ! gradient of  tracers 
    420                   pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
    421                ENDIF 
    422                ! 
    423                ! j- direction 
    424                IF( ze3wv >= 0._wp ) THEN      ! case 1 
    425                   zmaxv =  ze3wv / fse3w(ji,jj+1,ikv+1) 
    426                   ! interpolated values of tracers 
    427                   ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) ) 
    428                   ! gradient of tracers 
    429                   pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
    430                ELSE                           ! case 2 
    431                   zmaxv =  - ze3wv / fse3w(ji,jj,ikv+1) 
    432                   ! interpolated values of tracers 
    433                   ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) ) 
    434                   ! gradient of tracers 
    435                   pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
    436                ENDIF 
    437             END DO!! 
    438          END DO!! 
    439          CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
    440          ! 
    441       END DO 
    442  
    443       IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part) 
    444          ! 
    445          pgrui(:,:)  =0.0_wp ; pgrvi(:,:)  =0.0_wp ; 
    446          pgzui(:,:)  =0.0_wp ; pgzvi(:,:)  =0.0_wp ; 
    447          pmrui(:,:)  =0.0_wp ; pmrui(:,:)  =0.0_wp ; 
    448          pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ; 
    449          ! 
    450          DO jj = 1, jpjm1        ! depth of the partial step level 
    451             DO ji = 1, jpim1 
    452                iku = miku(ji,jj) 
    453                ikv = mikv(ji,jj) 
    454                ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
    455                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)) 
    456                ! 
    457                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1 
    458                ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2 
    459                ENDIF 
    460                IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1 
    461                ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2 
    462                ENDIF 
    463             END DO 
    464          END DO 
    465          ! 
    466          CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj  
    467          CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj  
    468          ! 
    469          DO jj = 1, jpjm1                 ! Gradient of density at the last level  
    470             DO ji = 1, jpim1 
    471                iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1 
    472                ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1 
    473                ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
    474                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)) 
    475                IF( ze3wu >= 0._wp ) THEN 
    476                  pgzui  (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku) 
    477                  pgrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1 
    478                  pmrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1  
    479                  pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                  & 
    480                     &           * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   & 
    481                     &              - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1 
    482                ELSE 
    483                  pgzui  (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu) 
    484                  pgrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) - zri(ji,jj) )      ! i: 2 
    485                  pmrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2 
    486                  pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                   & 
    487                     &           * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  & 
    488                     &              -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2 
    489                ENDIF 
    490                IF( ze3wv >= 0._wp ) THEN 
    491                  pgzvi  (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv)  
    492                  pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )        ! j: 1 
    493                  pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1 
    494                  pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                  &  
    495                      &           * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  & 
    496                      &                - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1 
    497                                   ! + 2 due to the formulation in density and not in anomalie in hpg sco 
    498                ELSE 
    499                  pgzvi  (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv) 
    500                  pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )     ! j: 2 
    501                  pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2 
    502                  pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                   & 
    503                     &           * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) & 
    504                     &              -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2 
    505                ENDIF 
    506             END DO 
    507          END DO 
    508          CALL lbc_lnk( pgrui   , 'U', -1. )   ;   CALL lbc_lnk( pgrvi   , 'V', -1. )   ! Lateral boundary conditions 
    509          CALL lbc_lnk( pmrui   , 'U',  1. )   ;   CALL lbc_lnk( pmrvi   , 'V',  1. )   ! Lateral boundary conditions 
    510          CALL lbc_lnk( pgzui   , 'U', -1. )   ;   CALL lbc_lnk( pgzvi   , 'V', -1. )   ! Lateral boundary conditions 
    511          CALL lbc_lnk( pge3rui , 'U', -1. )   ;   CALL lbc_lnk( pge3rvi , 'V', -1. )   ! Lateral boundary conditions 
     438 
     439               iku = miku(ji,jj)  
     440               ikv = mikv(ji,jj)  
     441               ze3wu  =  fsdept_n(ji,jj,iku) - fsdept_n(ji+1,jj,iku) 
     442               ze3wv  =  fsdept_n(ji,jj,ikv) - fsdept_n(ji,jj+1,ikv)  
     443 
     444               IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj      ) - prd(ji,jj,iku) ) ! i: 1 
     445               ELSE                      ; pgrui(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj  ,iku) - zri(ji,jj    ) ) ! i: 2 
     446               ENDIF 
     447 
     448               IF( ze3wv >= 0._wp ) THEN ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( zrj(ji  ,jj      ) - prd(ji,jj,ikv) ) ! j: 1 
     449               ELSE                      ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( prd(ji  ,jj+1,ikv) - zrj(ji,jj    ) ) ! j: 2 
     450               ENDIF 
     451 
     452            END DO 
     453         END DO 
     454         CALL lbc_lnk( pgrui   , 'U', -1. ); CALL lbc_lnk( pgrvi   , 'V', -1. )   ! Lateral boundary conditions 
    512455         ! 
    513456      END IF   

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

    r4990 r6012  
    3131   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld    !: mixing layer depth (turbocline)      [m] 
    3232   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m] 
    33    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m] 
     33   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: depth of the last T-point inside the mixed layer 
    3434 
    3535   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
    … …  
    111111      END DO 
    112112      ! 
    113       ! w-level of the turbocline 
     113      ! w-level of the turbocline and mixing layer (iom_use) 
    114114      imld(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point 
    115115      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10  
    116116         DO jj = 1, jpj 
    117117            DO ji = 1, jpi 
    118                imkt = mikt(ji,jj) 
    119                IF( avt (ji,jj,jk) < avt_c )   imld(ji,jj) = MAX( imkt, jk )      ! Turbocline  
     118               IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) )   imld(ji,jj) = jk      ! Turbocline  
    120119            END DO 
    121120         END DO 
    … …  
    127126            iikn = nmln(ji,jj) 
    128127            imkt = mikt(ji,jj) 
    129             hmld (ji,jj) = ( fsdepw(ji,jj,iiki  ) - fsdepw(ji,jj,imkt )            )  * ssmask(ji,jj)    ! Turbocline depth  
    130             hmlp (ji,jj) = ( fsdepw(ji,jj,iikn  ) - fsdepw(ji,jj,MAX( imkt,nla10 ) ) ) * ssmask(ji,jj)    ! Mixed layer depth 
    131             hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt )            )  * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer 
     128            hmld (ji,jj) = fsdepw(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth  
     129            hmlp (ji,jj) = fsdepw(ji,jj,iikn ) * ssmask(ji,jj)    ! Mixed layer depth 
     130            hmlpt(ji,jj) = fsdept(ji,jj,iikn-1) * ssmask(ji,jj)    ! depth of the last T-point inside the mixed layer 
    132131         END DO 
    133132      END DO 
    134133      IF( .NOT.lk_offline ) THEN            ! no need to output in offline mode 
    135          CALL iom_put( "mldr10_1", hmlp )   ! mixed layer depth 
    136          CALL iom_put( "mldkz5"  , hmld )   ! turbocline depth 
     134         IF ( iom_use("mldr10_1") ) THEN 
     135            IF( .NOT. ln_isfcav ) CALL iom_put( "mldr10_1", hmlp )            ! mixed layer depth 
     136            IF(       ln_isfcav ) CALL iom_put( "mldr10_1", hmlp - risfdep)   ! mixed layer thickness 
     137         END IF 
     138         IF ( iom_use("mldkz5") ) THEN 
     139            IF( .NOT. ln_isfcav ) CALL iom_put( "mldkz5"  , hmld )             ! turbocline depth 
     140            IF(       ln_isfcav ) CALL iom_put( "mldkz5"  , hmld - risfdep )   ! turbocline thickness 
     141         END IF 
    137142      ENDIF 
    138143       

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/oce.F90

    r5930 r6012  
    6060   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gtsu, gtsv   !: horizontal gradient of T, S bottom u-point  
    6161   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gru , grv    !: horizontal gradient of rd at bottom u-point 
    62    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aru , arv     
    63    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gzu , gzv     
    64    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ge3ru, ge3rv   !: horizontal gradient of T, S and rd at top v-point   
    6562 
    6663   !! (ISF) interpolated gradient (only used for ice shelf case)  
    … …  
    6865   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gtui, gtvi   !: horizontal gradient of T, S and rd at top u-point  
    6966   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   grui, grvi   !: horizontal gradient of T, S and rd at top v-point   
    70    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   arui, arvi   !: horizontal average  of rd          at top v-point   
    71    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gzui, gzvi   !: horizontal gradient of z           at top v-point   
    72    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ge3rui, ge3rvi   !: horizontal gradient of T, S and rd at top v-point   
    7367   !! (ISF) ice load 
    7468   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   riceload 
    … …  
    115109         &     spgu  (jpi,jpj)   , spgv(jpi,jpj)   ,                       & 
    116110         &     gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts),                     & 
    117          &     aru(jpi,jpj)      , arv(jpi,jpj)      ,                     & 
    118          &     gzu(jpi,jpj)      , gzv(jpi,jpj)      ,                     & 
    119111         &     gru(jpi,jpj)      , grv(jpi,jpj)      ,                     & 
    120          &     ge3ru(jpi,jpj)    , ge3rv(jpi,jpj)    ,                     & 
    121112         &     gtui(jpi,jpj,jpts), gtvi(jpi,jpj,jpts),                     & 
    122          &     arui(jpi,jpj)     , arvi(jpi,jpj)     ,                     & 
    123          &     gzui(jpi,jpj)     , gzvi(jpi,jpj)     ,                     & 
    124          &     ge3rui(jpi,jpj)   , ge3rvi(jpi,jpj)   ,                     & 
    125113         &     grui(jpi,jpj)     , grvi(jpi,jpj)     ,                     & 
    126114         &     riceload(jpi,jpj),                             STAT=ierr(2) ) 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/step.F90

    r6010 r6012  
    166166 
    167167         IF( ln_zps .AND.       ln_isfcav)                               & 
    168             &            CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv, gtui, gtvi,   &    ! Partial steps for top cell (ISF) 
    169             &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
    170             &                                   grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level 
    171  
     168            &            CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
     169            &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level 
    172170         IF( ln_traldf_triad ) THEN  
    173171                         CALL ldf_slp_triad( kstp )                       ! before slope for triad operator 
    … …  
    186184      IF( lk_vvl     )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors  
    187185                         CALL wzv           ( kstp )  ! now cross-level velocity  
    188  
    189186!!gm : why also here ???? 
    190187      IF(ln_sto_eos  )   CALL sto_pts( tsn )                             ! Random T/S fluctuations 
    … …  
    196193!!                                         but ensures reproductible results 
    197194!!                                         with previous versions using split-explicit free surface           
    198             IF( ln_zps .AND. .NOT. ln_isfcav)   &                           ! Partial steps: bottom before horizontal gradient 
    199                &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &  ! of t, s, rd at the last ocean level 
    200                &                                          rhd, gru , grv    ) 
    201             IF( ln_zps .AND.       ln_isfcav)   &                           ! Partial steps: top & bottom before horizontal gradient 
    202                &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,   &  
    203                &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   & 
    204                &                                               grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) 
     195            IF( ln_zps .AND. .NOT. ln_isfcav )                               & 
     196               &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,   &  ! Partial steps: before horizontal gradient 
     197               &                                          rhd, gru , grv     )  ! of t, s, rd at the last ocean level 
     198            IF( ln_zps .AND.       ln_isfcav )                                          & 
     199               &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF) 
     200               &                                          rhd, gru , grv , grui, grvi   )  ! of t, s, rd at the first ocean level 
    205201!!jc: fs simplification 
    206202                             

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/TOP_SRC/trcini.F90

    r5836 r6012  
    246246   END SUBROUTINE trc_ini_state 
    247247 
    248  
    249248   SUBROUTINE top_alloc 
    250249      !!---------------------------------------------------------------------- 

  • branches/2015/dev_MetOffice_merge_2015/NEMOGCM/SETTE/input_ORCA2_LIM.cfg

    r5398 r6012  
    1 ORCA2_LIM_nemo_v3.6.tar ORCA2_LIM_nemo_v3.6 
     1ORCA2_LIM_nemo_v3.6st.tar ORCA2_LIM_nemo_v3.6 
Note: See TracChangeset for help on using the changeset viewer.