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 4666 for branches – NEMO

Changeset 4666 for branches


Ignore:
Timestamp:
2014-06-11T14:52:23+02:00 (10 years ago)
Author:
mathiot
Message:

#1331 : add ISOMIP config files + ice shelf code

Location:
branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM
Files:
9 added
51 edited

Legend:

Unmodified
Added
Removed

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/CONFIG/SHARED/field_def.xml

    r4565 r4666  
    156156    <!-- * variable relative to atmospheric pressure forcing : available with ln_apr_dyn --> 
    157157    <field id="ssh_ib"       long_name="Inverse barometer sea surface height"                         unit="m"        /> 
    158     
     158 
     159         <!-- * variable related to ice shelf forcing * --> 
     160         <field id="fwfisf"       long_name="Ice shelf melting"                                            unit="Kg/m2/s"  /> 
     161         <field id="qisf"         long_name="Ice Shelf Heat Flux"                                          unit="W/m2"     /> 
     162         <field id="isfgammat"    long_name="transfert coefficient for isf (temperature) "                 unit="m/s"      /> 
     163         <field id="isfgammas"    long_name="transfert coefficient for isf (salinity)    "                 unit="m/s"      /> 
     164    <field id="stbl"         long_name="salinity in the Losh tbl                    "                 unit="PSU"      /> 
     165    <field id="ttbl"         long_name="temperature in the Losh tbl                 "                 unit="C"        /> 
     166 
    159167         <!-- *_oce variables available with ln_blk_clio or ln_blk_core --> 
    160168         <field id="qns_oce"      long_name="Non solar Downward Heat Flux over open ocean"                 unit="W/m2"     /> 
    … …  
    271279         <field id="suoce"        long_name="ocean surface current along i-axis"          unit="m/s"  /> 
    272280         <field id="uoce"         long_name="ocean current along i-axis"                  unit="m/s"  grid_ref="grid_U_3D" /> 
     281         <field id="ssu"          long_name="ocean surface current along i-axis"          unit="m/s"  grid_ref="none"      /> 
    273282         <field id="uocetr_eff"   long_name="Effective ocean transport along i-axis"      unit="m3/s" grid_ref="grid_U_3D" /> 
    274283         <field id="uocet"        long_name="ocean transport along i-axis times temperature" unit="degC.m/s" grid_ref="grid_U_3D" /> 
    … …  
    281290         <field id="uoce_bbl"     long_name="BBL ocean current along i-axis"              unit="m/s"  grid_ref="grid_U_3D" /> 
    282291    <field id="ahu_bbl"      long_name="BBL diffusive flux along i-axis"             unit="m3/s" /> 
     292         <!-- variable for ice shelves --> 
     293    <field id="utbl"         long_name="zonal current in the Losh tbl"               unit="m/s"  axis_ref="none" /> 
    283294         <!-- variables available with key_diaar5 --> 
    284295         <field id="u_masstr"     long_name="ocean eulerian mass transport along i-axis"  unit="kg/s" grid_ref="grid_U_3D" /> 
    … …  
    294305         <field id="svoce"        long_name="ocean surface current along j-axis"          unit="m/s"  /> 
    295306         <field id="voce"         long_name="ocean current along j-axis"                  unit="m/s"  grid_ref="grid_V_3D" /> 
     307         <field id="ssv"          long_name="ocean surface current along j-axis"          unit="m/s"  grid_ref="none" /> 
    296308         <field id="vocetr_eff"   long_name="Effective ocean transport along j-axis"      unit="m3/s" grid_ref="grid_V_3D" /> 
    297309         <field id="vocet"        long_name="ocean transport along j-axis times temperature" unit="degC.m/s" grid_ref="grid_V_3D" /> 
    … …  
    304316         <field id="voce_bbl"     long_name="BBL ocean current along j-axis"              unit="m/s"  grid_ref="grid_V_3D" /> 
    305317    <field id="ahv_bbl"      long_name="BBL diffusive flux along j-axis"             unit="m3/s"                   /> 
     318         <!-- variable for ice shelves --> 
     319    <field id="vtbl"         long_name="meridional current in the Losh tbl"        unit="m/s"  axis_ref="none" /> 
    306320         <!-- variables available with key_diaar5 --> 
    307321         <field id="v_masstr"     long_name="ocean eulerian mass transport along j-axis"  unit="kg/s" grid_ref="grid_V_3D" /> 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/CONFIG/SHARED/namelist_ref

    r4384 r4666  
    233233   ln_dm2dc    = .false.   !  daily mean to diurnal cycle on short wave 
    234234   ln_rnf      = .true.    !  runoffs                                   (T => fill namsbc_rnf) 
     235   nn_isf      = 1         !  0=no isf                   1 = presence of ISF  
     236                           !  2 = bg03 parametrisation   3 = rnf file for isf    
     237                           !  4 = ISF fwf specified 
    235238   ln_ssr      = .true.    !  Sea Surface Restoring on T and/or S       (T => fill namsbc_ssr) 
    236239   nn_fwb      = 2         !  FreshWater Budget: =0 unchecked 
    … …  
    396399   ln_rnf_tem   = .false.   !  read in temperature information for runoff 
    397400   ln_rnf_sal   = .false.   !  read in salinity information for runoff 
     401/ 
     402!----------------------------------------------------------------------- 
     403&namsbc_isf    !  Top boundary layer (ISF)  
     404!----------------------------------------------------------------------- 
     405!              ! file name ! frequency (hours) ! variable ! time interpol. !  clim   ! 'yearly'/ ! weights  ! rotation ! 
     406!              !           !  (if <0  months)  !   name   !    (logical)   !  (T/F)  ! 'monthly' ! filename ! pairing  ! 
     407! nn_isf == 4 
     408   sn_qisf      = 'rnfisf' ,         -12      ,'sohflisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   '' 
     409   sn_fwfisf    = 'rnfisf' ,         -12      ,'sowflisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   '' 
     410! nn_isf == 3 
     411   sn_rnfisf    = 'runoffs' ,         -12      ,'sofwfisf',    .false.      , .true.  , 'yearly'  ,  ''      ,   '' 
     412! nn_isf == 2 and 3 
     413   sn_depmax_isf = 'runoffs' ,       -12        ,'sozisfmax' ,   .false.  , .true.  , 'yearly'  ,  ''      ,   '' 
     414   sn_depmin_isf = 'runoffs' ,       -12        ,'sozisfmin' ,   .false.  , .true.  , 'yearly'  ,  ''      ,   '' 
     415! nn_isf == 2 
     416   sn_Leff_isf = 'rnfisf' ,       0          ,'Leff'         ,   .false.  , .true.  , 'yearly'  ,  ''      ,   '' 
     417! for all case 
     418   ln_divisf   = .true.  ! apply isf melting as a mass flux or in the salinity trend. (maybe I should remove this option as for runoff?) 
     419! only for nn_isf = 1 or 2 
     420   rn_gammat0  = 1.0e-4   ! gammat coefficient used in blk formula 
     421   rn_gammas0  = 1.0e-4   ! gammas coefficient used in blk formula 
     422! only for nn_isf = 1 
     423   nn_isfblk   =  1       ! 1 ISOMIP ; 2 conservative (3 equation formulation, Jenkins et al. 1991 ??) 
     424   rn_hisf_tbl =  30.      ! thickness of the top boundary layer           (Losh et al. 2008) 
     425                          ! 0 => thickness of the tbl = thickness of the first wet cell 
     426   ln_conserve = .true.   ! conservative case (take into account meltwater advection) 
     427   nn_gammablk = 1        ! 0 = cst Gammat (= gammat/s) 
     428                          ! 1 = velocity dependend Gamma (u* * gammat/s)  (Jenkins et al. 2010) 
     429                          !     if you want to keep the cd as in global config, adjust rn_gammat0 to compensate 
     430                          ! 2 = velocity and stability dependent Gamma    Holland et al. 1999 
    398431/ 
    399432!----------------------------------------------------------------------- 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r4624 r4666  
    759759            ENDIF 
    760760 
    761             tsb(:,:,:,:) = tsn(:,:,:,:)               ! Update before fields 
    762  
    763             CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )                ! Before potential and in situ densities 
     761            tsb(:,:,:,:) = tsn(:,:,:,:)                 ! Update before fields 
     762 
     763            CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )  ! Before potential and in situ densities 
    764764          
    765765            IF( ln_zps .AND. .NOT. lk_c1d ) & 
    766                &  CALL zps_hde( nit000, jpts, tsb, &  ! Partial steps: before horizontal derivative 
    767                &                gtsu, gtsv, rhd,   &  ! of T, S, rd at the bottom ocean level 
    768                &                gru , grv ) 
     766               &  CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
     767               &                rhd, gru , grv, aru, arv, gzu, gzv, ge3ru, ge3rv,  &             ! 
     768               &                gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
    769769 
    770770#if defined key_zdfkpp 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r4624 r4666  
    11521152         END DO 
    11531153 
    1154          tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:) 
     1154         tmask_i (:,:) = lmask(:,:) * tmask_i(:,:) 
    11551155 
    11561156      ENDIF ! ln_mask_file=.TRUE. 
    11571157       
    1158       bdytmask(:,:) = tmask(:,:,1) 
     1158      bdytmask(:,:) = lmask(:,:) 
    11591159      IF( .not. ln_mask_file ) THEN 
    11601160         ! If .not. ln_mask_file then we need to derive mask on U and V grid  

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90

    r4624 r4666  
    2424   USE lib_fortran     ! glob_sum 
    2525   USE restart         ! ocean restart 
    26    USE wrk_nemo         ! work arrays 
    27    USE sbcrnf         ! river runoffd 
     26   USE wrk_nemo        ! work arrays 
     27   USE sbcrnf          ! river runoff 
     28   USE sbcisf          ! ice shelves 
    2829 
    2930   IMPLICIT NONE 
    … …  
    8485      ! 1 - Trends due to forcing ! 
    8586      ! ------------------------- ! 
    86       z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) ) * surf(:,:) ) ! volume fluxes 
     87      z_frc_trd_v = r1_rau0 * glob_sum( - ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:)) * surf(:,:) ) ! volume fluxes 
    8788      z_frc_trd_t =           glob_sum( sbc_tsc(:,:,jp_tem) * surf(:,:) )       ! heat fluxes 
    8889      z_frc_trd_s =           glob_sum( sbc_tsc(:,:,jp_sal) * surf(:,:) )       ! salt fluxes 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90

    r4624 r4666  
    505505            btmsk(:,:,5) = MAX ( btmsk(:,:,3), btmsk(:,:,4) )          ! Indo-Pacific basin 
    506506            WHERE( gphit(:,:) < -30._wp)   ;   btm30(:,:) = 0._wp      ! mask out Southern Ocean 
    507             ELSE WHERE                     ;   btm30(:,:) = tmask(:,:,1) 
     507            ELSE WHERE                     ;   btm30(:,:) = lmask(:,:) 
    508508            END WHERE 
    509509         ENDIF 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r4570 r4666  
    149149         z3d(:,:,:) = tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) 
    150150         CALL iom_put( "toce" , z3d                        )   ! heat content 
    151          CALL iom_put( "sst"  , z3d(:,:,1)                 )   ! sea surface heat content 
    152          z3d(:,:,1) = tsn(:,:,1,jp_tem) * z3d(:,:,1) 
    153          CALL iom_put( "sst2" , z3d(:,:,1)                 )   ! sea surface content of squared temperature 
     151         DO jj = 1, jpj 
     152            DO ji = 1, jpi 
     153               z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * fse3t_n(ji,jj,mikt(ji,jj)) 
     154            END DO 
     155         END DO   
     156         CALL iom_put( "sst"  , z2d(:,:)                 )   ! sea surface heat content       
     157         DO jj = 1, jpj 
     158            DO ji = 1, jpi 
     159               z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem)**2 * fse3t_n(ji,jj,mikt(ji,jj)) 
     160            END DO 
     161         END DO   
     162         CALL iom_put( "sst2" , z2d(:,:)      )   ! sea surface content of squared temperature 
    154163         z3d(:,:,:) = tsn(:,:,:,jp_sal) * fse3t_n(:,:,:)             
    155164         CALL iom_put( "soce" , z3d                        )   ! salinity content 
    156          CALL iom_put( "sss"  , z3d(:,:,1)                 )   ! sea surface salinity content 
    157          z3d(:,:,1) = tsn(:,:,1,jp_sal) * z3d(:,:,1) 
    158          CALL iom_put( "sss2" , z3d(:,:,1)                 )   ! sea surface content of squared salinity 
     165         DO jj = 1, jpj 
     166            DO ji = 1, jpi 
     167               z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * fse3t_n(ji,jj,mikt(ji,jj)) 
     168            END DO 
     169         END DO   
     170         CALL iom_put( "sss"  , z2d(:,:)                 )   ! sea surface salinity content 
     171         DO jj = 1, jpj 
     172            DO ji = 1, jpi 
     173               z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal)**2 * fse3t_n(ji,jj,mikt(ji,jj)) 
     174            END DO 
     175         END DO   
     176         CALL iom_put( "sss2" , z2d(:,:)                 )   ! sea surface content of squared salinity 
    159177      ELSE 
    160          CALL iom_put( "toce" , tsn(:,:,:,jp_tem)          )   ! temperature 
    161          CALL iom_put( "sst"  , tsn(:,:,1,jp_tem)          )   ! sea surface temperature 
    162          CALL iom_put( "sst2" , tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) ) ! square of sea surface temperature 
     178         CALL iom_put( "toce" , tsn(:,:,:,jp_tem)        )   ! temperature 
     179         DO jj = 1, jpj 
     180            DO ji = 1, jpi 
     181               z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) 
     182            END DO 
     183         END DO          
     184         CALL iom_put( "sst"  , z2d(:,:)            ) ! sea surface temperature 
     185         CALL iom_put( "sst2" , z2d(:,:) * z2d(:,:) ) ! square of sea surface temperature 
    163186         CALL iom_put( "soce" , tsn(:,:,:,jp_sal)          )   ! salinity 
    164          CALL iom_put( "sss"  , tsn(:,:,1,jp_sal)          )   ! sea surface salinity 
    165          CALL iom_put( "sss2" , tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) ) ! square of sea surface salinity 
     187         DO jj = 1, jpj 
     188            DO ji = 1, jpi 
     189               z2d(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) 
     190            END DO 
     191         END DO 
     192         CALL iom_put( "sss"  , z2d(:,:)            ) ! sea surface salinity 
     193         CALL iom_put( "sss2" , z2d(:,:) * z2d(:,:) ) ! square of sea surface salinity 
    166194      END IF 
    167195      IF( lk_vvl .AND. (.NOT. ln_dynadv_vec) ) THEN 
    … …  
    171199         CALL iom_put( "uoce" , un                         )    ! i-current 
    172200         CALL iom_put( "voce" , vn                         )    ! j-current 
     201         DO jj = 1, jpj 
     202            DO ji = 1, jpi 
     203               z2d(ji,jj) = un(ji,jj,miku(ji,jj)) 
     204            END DO 
     205         END DO 
     206         CALL iom_put( "ssu"   , z2d                                    )    ! i-current 
     207         DO jj = 1, jpj 
     208            DO ji = 1, jpi 
     209               z2d(ji,jj) = vn(ji,jj,mikv(ji,jj)) 
     210            END DO 
     211         END DO 
     212         CALL iom_put( "ssv"   , z2d                                    )    ! j-current 
    173213      END IF 
    174214      CALL iom_put(    "avt"  , avt                        )    ! T vert. eddy diff. coef. 
    … …  
    623663      ENDIF 
    624664 
    625       ! Write fields on T grid 
    626665      IF( lk_vvl ) THEN 
    627666         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content 
    … …  
    634673         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature 
    635674         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity 
    636  
    637675      ENDIF 
    638676      IF( lk_vvl ) THEN 
    … …  
    685723 
    686724#if ! defined key_coupled 
    687       CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping 
    688       CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping 
    689       IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1) 
    690       CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
     725!      CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping 
     726!      CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping 
     727!      IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1) 
     728!      CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
    691729#endif 
    692730#if ( defined key_coupled && ! defined key_lim3 && ! defined key_lim2 )  
    … …  
    696734      CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping 
    697735#endif 
    698       zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1) 
    699       CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ??? 
     736!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1) 
     737!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ??? 
    700738 
    701739#if defined key_diahth 
    … …  
    717755# endif 
    718756#endif 
    719          ! Write fields on U grid 
     757      ! Write fields on U grid  
    720758      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current 
    721759      IF( ln_traldf_gdia ) THEN 
    … …  
    739777      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress 
    740778 
    741          ! Write fields on V grid 
    742779      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current 
    743780      IF( ln_traldf_gdia ) THEN 
    … …  
    754791      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress 
    755792 
    756          ! Write fields on W grid 
    757793      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current 
    758794      IF( ln_traldf_gdia ) THEN 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90

    r4488 r4666  
    175175   LOGICAL, PUBLIC ::   ln_zps        !: z-coordinate - partial step 
    176176   LOGICAL, PUBLIC ::   ln_sco        !: s-coordinate or hybrid z-s coordinate 
     177   LOGICAL, PUBLIC ::   ln_isf        !: presence of ISF  
    177178 
    178179   !! All coordinates 
    … …  
    250251   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbkt               !: vertical index of the bottom last T- ocean level 
    251252   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbku, mbkv         !: vertical index of the bottom last U- and W- ocean level 
    252    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy              !: ocean depth (meters) 
    253    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i            !: interior domain T-point mask 
    254    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask              !: land/ocean mask of barotropic stream function 
     253   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy                              !: ocean depth (meters) 
     254   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmask_i, umask_i, vmask_i, fmask_i !: interior domain T-point mask 
     255   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bmask                              !: land/ocean mask of barotropic stream function 
     256 
     257   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   micedep                 !: top first ocean level                (ISF) 
     258   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mikt, miku, mikv, mikf  !: first wet T-, U-, V-, F- ocean level (ISF) 
     259   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   icedep, iceh            !: Iceshelf draft, iceshef freeboard    (ISF) 
     260   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   lmask                   !: surface domain T-point mask  
    255261 
    256262   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:), TARGET :: tmask, umask, vmask, fmask   !: land/ocean mask at T-, U-, V- and F-pts 
    … …  
    379385         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(8) ) 
    380386 
    381       ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) ,                     & 
    382          &     tmask_i(jpi,jpj) , bmask(jpi,jpj) ,                     & 
     387      ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) ,                                      & 
     388         &     tmask_i(jpi,jpj) , umask_i(jpi,jpj), vmask_i(jpi,jpj), fmask_i(jpi,jpj), & 
     389         &     bmask(jpi,jpj)   ,                                                       & 
    383390         &     mbkt   (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) ) 
    384391 
     392! (ISF) Allocation of basic array    
     393      ALLOCATE( micedep(jpi,jpj) , icedep(jpi,jpj), iceh(jpi,jpj),      & 
     394         &     mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj) ,           & 
     395         &     mikf(jpi,jpj), lmask(jpi,jpj), STAT=ierr(10) ) 
     396 
    385397      ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk),     &  
    386          &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk), STAT=ierr(10) ) 
     398         &      vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk), STAT=ierr(11) ) 
    387399 
    388400#if defined key_noslip_accurate 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90

    r4624 r4666  
    111111      END DO 
    112112      !                                        ! Inverse of the local depth 
    113       hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1) 
    114       hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1) 
     113      hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 
     114      hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 
    115115 
    116116                             CALL dom_stp      ! time step 
    … …  
    159159      READ  ( numnam_cfg, namrun, IOSTAT = ios, ERR = 902 ) 
    160160902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namrun in configuration namelist', lwp ) 
    161       IF(lwm) WRITE ( numond, namrun ) 
     161      WRITE ( numond, namrun ) 
    162162      ! 
    163163      IF(lwp) THEN                  ! control print 
    … …  
    241241      READ  ( numnam_cfg, namdom, IOSTAT = ios, ERR = 904 ) 
    242242904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp ) 
    243       IF(lwm) WRITE ( numond, namdom ) 
     243      WRITE ( numond, namdom ) 
    244244 
    245245      IF(lwp) THEN 
    … …  
    303303      READ  ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 ) 
    304304906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp ) 
    305       IF(lwm) WRITE( numond, namcla ) 
     305      WRITE( numond, namcla ) 
    306306 
    307307      IF(lwp) THEN 
    … …  
    327327      READ  ( numnam_cfg, namnc4, IOSTAT = ios, ERR = 908 ) 
    328328908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namnc4 in configuration namelist', lwp ) 
    329       IF(lwm) WRITE( numond, namnc4 ) 
     329      WRITE( numond, namnc4 ) 
    330330 
    331331      IF(lwp) THEN                        ! control print 
    … …  
    365365      ! 
    366366      IF(lk_mpp) THEN 
    367          CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 ) 
    368          CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 ) 
    369          CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 ) 
    370          CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 ) 
     367         CALL mpp_minloc( e1t(:,:), tmask_i(:,:), ze1min, iimi1,ijmi1 ) 
     368         CALL mpp_minloc( e2t(:,:), tmask_i(:,:), ze2min, iimi2,ijmi2 ) 
     369         CALL mpp_maxloc( e1t(:,:), tmask_i(:,:), ze1max, iima1,ijma1 ) 
     370         CALL mpp_maxloc( e2t(:,:), tmask_i(:,:), ze2max, iima2,ijma2 ) 
    371371      ELSE 
    372          ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1._wp )     
    373          ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1._wp )     
    374          ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1._wp )     
    375          ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1._wp )     
    376  
    377          iloc  = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp ) 
     372         ze1min = MINVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )     
     373         ze2min = MINVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )     
     374         ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp )     
     375         ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp )     
     376 
     377         iloc  = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp ) 
    378378         iimi1 = iloc(1) + nimpp - 1 
    379379         ijmi1 = iloc(2) + njmpp - 1 
    380          iloc  = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp ) 
     380         iloc  = MINLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp ) 
    381381         iimi2 = iloc(1) + nimpp - 1 
    382382         ijmi2 = iloc(2) + njmpp - 1 
    383          iloc  = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1._wp ) 
     383         iloc  = MAXLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp ) 
    384384         iima1 = iloc(1) + nimpp - 1 
    385385         ijma1 = iloc(2) + njmpp - 1 
    386          iloc  = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1._wp ) 
     386         iloc  = MAXLOC( e2t(:,:), mask = tmask_i(:,:) == 1._wp ) 
    387387         iima2 = iloc(1) + nimpp - 1 
    388388         ijma2 = iloc(2) + njmpp - 1 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90

    r4366 r4666  
    638638      CALL iom_close( inum ) 
    639639       
     640! need to be define for the extended grid south of -80S 
     641! some point are undefined but you need to have e1 and e2 .NE. 0 
     642      WHERE (e1t==0.0_wp) 
     643         e1t=1.0e2 
     644      END WHERE 
     645      WHERE (e1v==0.0_wp) 
     646         e1v=1.0e2 
     647      END WHERE 
     648      WHERE (e1u==0.0_wp) 
     649         e1u=1.0e2 
     650      END WHERE 
     651      WHERE (e1f==0.0_wp) 
     652         e1f=1.0e2 
     653      END WHERE 
     654      WHERE (e2t==0.0_wp) 
     655         e2t=1.0e2 
     656      END WHERE 
     657      WHERE (e2v==0.0_wp) 
     658         e2v=1.0e2 
     659      END WHERE 
     660      WHERE (e2u==0.0_wp) 
     661         e2u=1.0e2 
     662      END WHERE 
     663      WHERE (e2f==0.0_wp) 
     664         e2f=1.0e2 
     665      END WHERE 
     666        
    640667    END SUBROUTINE hgr_read 
    641668     

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90

    r4624 r4666  
    184184         END DO   
    185185      END DO   
     186       
     187      ! (ISF) define barotropic mask and mask the ice shelf point 
     188      lmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked 
     189       
     190      DO jk = 1, jpk 
     191         DO jj = 1, jpj 
     192            DO ji = 1, jpi 
     193               IF( REAL( micedep(ji,jj) - jk, wp ) - 0.1_wp >= 0._wp )   THEN 
     194                  tmask(ji,jj,jk) = 0._wp 
     195               END IF 
     196            END DO   
     197         END DO   
     198      END DO   
    186199 
    187200!!gm  ???? 
    … …  
    207220      ! Interior domain mask (used for global sum) 
    208221      ! -------------------- 
    209       tmask_i(:,:) = tmask(:,:,1) 
     222      tmask_i(:,:) = lmask(:,:)            ! (ISH) tmask_i = 1 even on the ice shelf 
    210223      iif = jpreci                         ! ??? 
    211224      iil = nlci - jpreci + 1 
    … …  
    250263         END DO 
    251264      END DO 
     265      ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at least 1 wet u point 
     266      DO jj = 1, jpjm1 
     267         DO ji = 1, fs_jpim1   ! vector loop 
     268            umask_i(ji,jj)  = lmask(ji,jj) * lmask(ji+1,jj  )  * MIN(1._wp,SUM(umask(ji,jj,:))) 
     269            vmask_i(ji,jj)  = lmask(ji,jj) * lmask(ji  ,jj+1)  * MIN(1._wp,SUM(vmask(ji,jj,:))) 
     270         END DO 
     271         DO ji = 1, jpim1      ! NO vector opt. 
     272            fmask_i(ji,jj) =  lmask(ji,jj  ) * lmask(ji+1,jj  )   & 
     273               &            * lmask(ji,jj+1) * lmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 
     274         END DO 
     275      END DO 
    252276      CALL lbc_lnk( umask, 'U', 1._wp )      ! Lateral boundary conditions 
    253277      CALL lbc_lnk( vmask, 'V', 1._wp ) 
    254278      CALL lbc_lnk( fmask, 'F', 1._wp ) 
     279      CALL lbc_lnk( umask_i, 'U', 1._wp )      ! Lateral boundary conditions 
     280      CALL lbc_lnk( vmask_i, 'V', 1._wp ) 
     281      CALL lbc_lnk( fmask_i, 'F', 1._wp ) 
    255282 
    256283 
    257284      ! 4. ocean/land mask for the elliptic equation 
    258285      ! -------------------------------------------- 
    259       bmask(:,:) = tmask(:,:,1)       ! elliptic equation is written at t-point 
     286      bmask(:,:) = lmask(:,:)       ! elliptic equation is written at t-point 
    260287      ! 
    261288      !                               ! Boundary conditions 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90

    r4624 r4666  
    169169      fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1) 
    170170      fsdepw_b(:,:,1) = 0.0_wp 
    171       DO jk = 2, jpk 
    172          fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 
    173          fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 
    174          fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:) 
    175          fsdept_b(:,:,jk) = fsdept_b(:,:,jk-1) + fse3w_b(:,:,jk) 
    176          fsdepw_b(:,:,jk) = fsdepw_b(:,:,jk-1) + fse3t_b(:,:,jk-1) 
     171      DO jj = 1,jpj 
     172         DO ji = 1,jpi 
     173            DO jk = 1,mikt(ji,jj)-1 
     174               fsdept_n(ji,jj,jk) = gdept_0(ji,jj,jk) 
     175               fsdepw_n(ji,jj,jk) = gdepw_0(ji,jj,jk) 
     176               fsde3w_n(ji,jj,jk) = gdept_0(ji,jj,jk) - sshn(ji,jj) 
     177               fsdept_b(ji,jj,jk) = gdept_0(ji,jj,jk) 
     178               fsdepw_b(ji,jj,jk) = gdepw_0(ji,jj,jk) 
     179            END DO 
     180            DO jk = mikt(ji,jj), jpk 
     181               fsdept_n(ji,jj,jk) = fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk) 
     182               fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 
     183               fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk  ) - sshn   (ji,jj) 
     184               fsdept_b(ji,jj,jk) = fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk) 
     185               fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1) 
     186            END DO 
     187         END DO 
    177188      END DO 
    178189 
    … …  
    185196         hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 
    186197      END DO 
    187       hur_b(:,:) = umask(:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) ) 
    188       hvr_b(:,:) = vmask(:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) ) 
     198      hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) ) 
     199      hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) ) 
    189200 
    190201      ! Restoring frequencies for z_tilde coordinate 
    … …  
    293304      !                                                ! --------------------------------------------- ! 
    294305 
    295       z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) 
     306      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask_i(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask_i(:,:) ) 
    296307      DO jk = 1, jpkm1 
    297308         ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0) 
    … …  
    313324            zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 
    314325         END DO 
    315          zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask(:,:,1) ) 
     326         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 
    316327 
    317328         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only) 
    … …  
    338349         ELSE                         ! layer case 
    339350            DO jk = 1, jpkm1 
    340                tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) 
     351               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 
    341352            END DO 
    342353         END IF 
    … …  
    463474      !                                           ! ---baroclinic part--------- ! 
    464475         DO jk = 1, jpkm1 
    465             fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) 
     476            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 
    466477         END DO 
    467478      ENDIF 
    … …  
    531542      END DO 
    532543      !                                        ! Inverse of the local depth 
    533       hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1) 
    534       hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1) 
     544      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 
     545      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 
    535546 
    536547      CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 
    … …  
    830841               DO jk=1,jpk 
    831842                  fse3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 
    832                       &                            / ( ht_0(:,:) + 1._wp - tmask(:,:,1) ) * tmask(:,:,jk) & 
     843                      &                            / ( ht_0(:,:) + 1._wp - tmask_i(:,:) ) * tmask(:,:,jk) & 
    833844                      &            + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 
    834845               END DO 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90

    r4292 r4666  
    132132       
    133133      CALL dom_uniq( zprw, 'T' ) 
    134       zprt = tmask(:,:,1) * zprw                               !    ! unique point mask 
     134      DO jj = 1, jpj 
     135         DO ji = 1, jpi 
     136            jk=mikt(ji,jj)  
     137            zprt(ji,jj) = tmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask 
     138         END DO 
     139      END DO                             !    ! unique point mask 
    135140      CALL iom_rstput( 0, 0, inum2, 'tmaskutil', zprt, ktype = jp_i1 )   
    136141      CALL dom_uniq( zprw, 'U' ) 
    137       zprt = umask(:,:,1) * zprw 
     142      DO jj = 1, jpj 
     143         DO ji = 1, jpi 
     144            jk=miku(ji,jj)  
     145            zprt(ji,jj) = umask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask 
     146         END DO 
     147      END DO 
    138148      CALL iom_rstput( 0, 0, inum2, 'umaskutil', zprt, ktype = jp_i1 )   
    139149      CALL dom_uniq( zprw, 'V' ) 
    140       zprt = vmask(:,:,1) * zprw 
     150      DO jj = 1, jpj 
     151         DO ji = 1, jpi 
     152            jk=mikv(ji,jj)  
     153            zprt(ji,jj) = vmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask 
     154         END DO 
     155      END DO 
    141156      CALL iom_rstput( 0, 0, inum2, 'vmaskutil', zprt, ktype = jp_i1 )   
    142157      CALL dom_uniq( zprw, 'F' ) 
    143       zprt = fmask(:,:,1) * zprw 
     158      DO jj = 1, jpj 
     159         DO ji = 1, jpi 
     160            jk=mikf(ji,jj)  
     161            zprt(ji,jj) = fmask(ji,jj,jk) * zprw(ji,jj)                        !    ! unique point mask 
     162         END DO 
     163      END DO 
    144164      CALL iom_rstput( 0, 0, inum2, 'fmaskutil', zprt, ktype = jp_i1 )   
    145165 
    … …  
    168188       
    169189      ! note that mbkt is set to 1 over land ==> use surface tmask 
    170       zprt(:,:) = tmask(:,:,1) * REAL( mbkt(:,:) , wp ) 
     190      zprt(:,:) = lmask(:,:) * REAL( mbkt(:,:) , wp ) 
    171191      CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i2 )     !    ! nb of ocean T-points 
     192      zprt(:,:) = lmask(:,:) * REAL( mikt(:,:) , wp ) 
     193      CALL iom_rstput( 0, 0, inum4, 'misf', zprt, ktype = jp_i2 )       !    ! nb of ocean T-points 
     194      zprt(:,:) = lmask(:,:) * REAL( icedep(:,:) , wp ) 
     195      CALL iom_rstput( 0, 0, inum4, 'isfdraft', zprt, ktype = jp_r4 )       !    ! nb of ocean T-points 
    172196             
    173197      IF( ln_sco ) THEN                                         ! s-coordinate 
    … …  
    203227            DO jj = 1,jpj    
    204228               DO ji = 1,jpi 
    205                   e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1) 
    206                   e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * tmask(ji,jj,1) 
     229                  e3tp(ji,jj) = e3t_0(ji,jj,mbkt(ji,jj)) * lmask(ji,jj) 
     230                  e3wp(ji,jj) = e3w_0(ji,jj,mbkt(ji,jj)) * lmask(ji,jj) 
    207231               END DO 
    208232            END DO 
    … …  
    228252            DO jj = 1,jpj    
    229253               DO ji = 1,jpi 
    230                   zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj)  ) * tmask(ji,jj,1) 
    231                   zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * tmask(ji,jj,1) 
     254                  zprt(ji,jj) = gdept_0(ji,jj,mbkt(ji,jj)  ) * lmask(ji,jj) 
     255                  zprw(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) * lmask(ji,jj) 
    232256               END DO 
    233257            END DO 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90

    r4624 r4666  
    3535   USE oce               ! ocean variables 
    3636   USE dom_oce           ! ocean domain 
     37   USE sbc_oce           ! surface variable (isf) 
    3738   USE closea            ! closed seas 
    3839   USE c1d               ! 1D vertical configuration 
    … …  
    102103      ! 
    103104      NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco 
     105      NAMELIST/namsbc/ nn_fsbc   , ln_ana , ln_flx  , ln_blk_clio, ln_blk_core, ln_cpl,   & 
     106         &             ln_blk_mfs, ln_apr_dyn, nn_ice , ln_dm2dc, ln_rnf, ln_ssr, nn_fwb, ln_cdgw, nn_isf 
    104107      !!---------------------------------------------------------------------- 
    105108      ! 
    … …  
    145148      IF( .NOT.lk_c1d )   CALL zgr_bat_ctl      ! check bathymetry (mbathy) and suppress isolated ocean points 
    146149                          CALL zgr_bot_level    ! deepest ocean level for t-, u- and v-points 
     150                          CALL zgr_top_level    ! shallowest ocean level for T-, U-, V- points 
    147151      ! 
    148152      IF( lk_c1d ) THEN                         ! 1D config.: same mbathy value over the 3x3 domain 
    … …  
    294298         gdepw_1d(1) = 0._wp                    ! force first w-level to be exactly at zero 
    295299      ENDIF 
     300 
     301! need to be like this to compute the pressure gradient with ISF 
     302! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 
     303      DO jk = 1, jpkm1 
     304         e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)  
     305      END DO 
     306      e3t_1d(jpk) = e3t_1d(jpk-1)   ! we don't care because this level is masked in NEMO 
     307 
     308      DO jk = 2, jpk 
     309         e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)  
     310      END DO 
     311      e3w_1d(1  ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))  
    296312 
    297313!!gm BUG in s-coordinate this does not work! 
    … …  
    451467         END DO 
    452468         ! 
     469         ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code 
     470         IF( cp_cfg == "isomip" ) THEN  
     471           !  
     472           icedep(:,:)=200.e0  
     473           micedep(:,:)=1  
     474           ij0 = 1 ; ij1 = 40  
     475           DO jj = mj0(ij0), mj1(ij1)  
     476              icedep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp  
     477                END DO  
     478            WHERE( bathy(:,:) <= 0._wp )  icedep(:,:) = 0._wp  
     479           !  
     480         ELSEIF ( cp_cfg == "isomip2" ) THEN 
     481         !  
     482            icedep(:,:)=0.e0 
     483            micedep(:,:)=1 
     484            ij0 = 1 ; ij1 = 40 
     485            DO jj = mj0(ij0), mj1(ij1) 
     486               icedep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp 
     487            END DO 
     488            WHERE( bathy(:,:) <= 0._wp )  icedep(:,:) = 0._wp 
     489         END IF 
     490         ! 
    453491         !                                            ! ================ ! 
    454492      ELSEIF( ntopo == 1 ) THEN                       !   read in file   ! (over the local domain) 
    … …  
    492530            CALL iom_get  ( inum, jpdom_data, 'Bathymetry', bathy ) 
    493531            CALL iom_close( inum ) 
    494             !                                                 
     532            !   
     533            icedep(:,:)=0._wp          
     534            micedep(:,:)=1              
     535            IF ( nn_isf == 1 .OR. nn_isf == 4 ) THEN 
     536               CALL iom_open ( 'isf_draft_meter.nc', inum )  
     537               CALL iom_get  ( inum, jpdom_data, 'isf_draft', icedep ) 
     538               CALL iom_close( inum ) 
     539               WHERE( bathy(:,:) <= 0._wp )  icedep(:,:) = 0._wp 
     540            END IF 
     541            !        
    495542            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration 
    496543               ! 
    … …  
    783830   END SUBROUTINE zgr_bot_level 
    784831 
     832      SUBROUTINE zgr_top_level 
     833      !!---------------------------------------------------------------------- 
     834      !!                    ***  ROUTINE zgr_bot_level  *** 
     835      !! 
     836      !! ** Purpose :   defines the vertical index of ocean top (mik. arrays) 
     837      !! 
     838      !! ** Method  :   computes from micedep with a minimum value of 1 
     839      !! 
     840      !! ** Action  :   mikt, miku, mikv :   vertical indices of the shallowest  
     841      !!                                     ocean level at t-, u- & v-points 
     842      !!                                     (min value = 1) 
     843      !!---------------------------------------------------------------------- 
     844      !! 
     845      INTEGER ::   ji, jj   ! dummy loop indices 
     846      REAL(wp), POINTER, DIMENSION(:,:) ::  zmik 
     847      !!---------------------------------------------------------------------- 
     848      ! 
     849      IF( nn_timing == 1 )  CALL timing_start('zgr_top_level') 
     850      ! 
     851      CALL wrk_alloc( jpi, jpj, zmik ) 
     852      ! 
     853      IF(lwp) WRITE(numout,*) 
     854      IF(lwp) WRITE(numout,*) '    zgr_top_level : ocean top k-index of T-, U-, V- and W-levels ' 
     855      IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~' 
     856      ! 
     857      mikt(:,:) = MAX( micedep(:,:) , 1 )    ! top k-index of T-level (=1) 
     858      !                                      ! top k-index of W-level (=mikt) 
     859      DO jj = 1, jpjm1                       ! top k-index of U- (U-) level 
     860         DO ji = 1, jpim1 
     861            miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  ) 
     862            mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  ) 
     863            mikf(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj), mikt(ji+1,jj  ), mikt(ji+1,jj+1)  ) 
     864         END DO 
     865      END DO 
     866 
     867      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk  
     868      zmik(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk(zmik,'U',1.)   ;   miku  (:,:) = MAX( INT( zmik(:,:) ), 1 ) 
     869      zmik(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk(zmik,'V',1.)   ;   mikv  (:,:) = MAX( INT( zmik(:,:) ), 1 ) 
     870      zmik(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk(zmik,'F',1.)   ;   mikf  (:,:) = MAX( INT( zmik(:,:) ), 1 ) 
     871      ! 
     872      CALL wrk_dealloc( jpi, jpj, zmik ) 
     873      ! 
     874      IF( nn_timing == 1 )  CALL timing_stop('zgr_top_level') 
     875      ! 
     876   END SUBROUTINE zgr_top_level 
    785877 
    786878   SUBROUTINE zgr_zco 
    … …  
    861953      !!---------------------------------------------------------------------- 
    862954      !! 
    863       INTEGER  ::   ji, jj, jk       ! dummy loop indices 
     955      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    864956      INTEGER  ::   ik, it           ! temporary integers 
     957      INTEGER  ::   id, jd, nprocd 
     958      INTEGER  ::   icompt, ibtest   ! (ISF) 
    865959      LOGICAL  ::   ll_print         ! Allow  control print for debugging 
    866960      REAL(wp) ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points 
    867961      REAL(wp) ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t 
    868       REAL(wp) ::   zmax             ! Maximum depth 
     962      REAL(wp) ::   zmax, zmin       ! Maximum and minimum depth 
    869963      REAL(wp) ::   zdiff            ! temporary scalar 
    870964      REAL(wp) ::   zrefdep          ! temporary scalar 
     965      REAL(wp) ::   eps=0.99            ! small offset to avoid large pool in case bathy slightly greater than icedep 
     966      REAL(wp), POINTER, DIMENSION(:,:)   ::   zbathy, zmask   ! 3D workspace (ISH) 
     967      INTEGER , POINTER, DIMENSION(:,:)   ::   zmbathy, zmicedep   ! 3D workspace (ISH) 
    871968      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zprt 
    872969      !!--------------------------------------------------------------------- 
    … …  
    875972      ! 
    876973      CALL wrk_alloc( jpi, jpj, jpk, zprt ) 
     974      CALL wrk_alloc( jpi, jpj, zbathy, zmask) 
     975      CALL wrk_alloc( jpi, jpj, zmbathy, zmicedep) 
    877976      ! 
    878977      IF(lwp) WRITE(numout,*) 
    … …  
    9061005         WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth )   mbathy(:,:) = jk-1 
    9071006      END DO 
     1007      ! (ISF) compute micedep 
     1008      WHERE( icedep(:,:) == 0._wp ) ;   micedep(:,:) = 1   ! no ice shelf : set micedep to 1   
     1009      ELSEWHERE                     ;   micedep(:,:) = 2   ! iceshelf : initialize micedep to second level  
     1010      END WHERE   
     1011 
     1012      ! Compute micedep for ocean points (i.e. first wet level)  
     1013      ! find the first ocean level such that the first level thickness  
     1014      ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where  
     1015      ! e3t_0 is the reference level thickness  
     1016      DO jk = 2, jpkm1  
     1017         zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )  
     1018         WHERE( 0._wp < icedep(:,:) .AND. icedep(:,:) >= zdepth )   micedep(:,:) = jk+1  
     1019      END DO  
     1020      WHERE (icedep <= e3t_1d(1) .AND. icedep .GT. 0._wp) 
     1021         icedep = 0. 
     1022         micedep= 1 
     1023      END WHERE 
     1024     
     1025 
     1026      icompt = 0  
     1027      DO jl = 1, 10  
     1028        IF( lk_mpp ) THEN 
     1029           zbathy(:,:) = FLOAT( micedep(:,:) ) 
     1030           CALL lbc_lnk( zbathy, 'T', 1. ) 
     1031           micedep(:,:) = INT( zbathy(:,:) ) 
     1032           CALL lbc_lnk( icedep, 'T', 1. ) 
     1033           CALL lbc_lnk( bathy, 'T', 1. ) 
     1034           zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1035           CALL lbc_lnk( zbathy, 'T', 1. ) 
     1036           mbathy(:,:) = INT( zbathy(:,:) ) 
     1037        ENDIF 
     1038         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN  
     1039            micedep( 1 ,:) = micedep(jpim1,:)           ! local domain is cyclic east-west  
     1040            micedep(jpi,:) = micedep(  2  ,:)  
     1041         ENDIF 
     1042 
     1043         IF( nperio == 1 .OR. nperio  ==  4 .OR. nperio  ==  6 ) THEN 
     1044            mbathy( 1 ,:) = mbathy(jpim1,:)             ! local domain is cyclic east-west 
     1045            mbathy(jpi,:) = mbathy(  2  ,:) 
     1046         ENDIF 
     1047 
     1048         WHERE (mbathy == 0)  
     1049            icedep = 0._wp 
     1050            micedep= 0 
     1051            bathy  = 0._wp 
     1052         ENDWHERE 
     1053 
     1054         WHERE (bathy(:,:) < icedep(:,:)+eps) 
     1055            micedep(:,:) = 0  
     1056            icedep(:,:)  = 0._wp 
     1057            mbathy(:,:)  = 0 
     1058            bathy(:,:)   = 0._wp 
     1059         END WHERE 
     1060 
     1061         DO jj = 1, jpj 
     1062            DO ji = 1, jpi 
     1063               IF (bathy(ji,jj) .GT. icedep(ji,jj) .AND. mbathy(ji,jj) .LT. micedep(ji,jj)) THEN 
     1064                  bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 
     1065                  mbathy(ji,jj)= mbathy(ji,jj) + 1 
     1066               END IF 
     1067            END DO 
     1068         END DO 
     1069 
     1070 
     1071         ! At least 2 levels for water thickness at T, U, and V point. 
     1072        zmicedep(:,:)=micedep(:,:) 
     1073        zmbathy (:,:)=mbathy (:,:) 
     1074 
     1075         DO jj = 2, jpjm1 
     1076            DO ji = 2, jpim1 
     1077               IF( zmicedep(ji,jj) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN 
     1078                  mbathy(ji,jj) = zmbathy(ji,jj) + 1  
     1079                  bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat ) 
     1080               ENDIF 
     1081               IF( zmicedep(ji,jj+1) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN 
     1082                  mbathy(ji,jj) = zmbathy(ji,jj) + 1  
     1083                  bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat ) 
     1084               ENDIF 
     1085               IF( zmicedep(ji,jj-1) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN 
     1086                  mbathy(ji,jj) = zmbathy(ji,jj) + 1  
     1087                  bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat ) 
     1088               ENDIF 
     1089               IF( zmicedep(ji+1,jj) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN 
     1090                  mbathy(ji,jj) = zmbathy(ji,jj) + 1  
     1091                  bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat ) 
     1092               ENDIF 
     1093               IF( zmicedep(ji-1,jj) == zmbathy(ji,jj) .AND. zmbathy(ji,jj) .GT. 1) THEN 
     1094                  mbathy(ji,jj) = zmbathy(ji,jj) + 1  
     1095                  bathy(ji,jj)=gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj))*e3zps_rat ) 
     1096               ENDIF 
     1097            ENDDO 
     1098         ENDDO 
     1099   IF( lk_mpp ) THEN  
     1100      zbathy(:,:) = FLOAT( micedep(:,:) )  
     1101      CALL lbc_lnk( zbathy, 'T', 1. )  
     1102      micedep(:,:) = INT( zbathy(:,:) )  
     1103      CALL lbc_lnk( icedep, 'T', 1. )  
     1104           CALL lbc_lnk( bathy, 'T', 1. ) 
     1105           zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1106           CALL lbc_lnk( zbathy, 'T', 1. ) 
     1107           mbathy(:,:) = INT( zbathy(:,:) ) 
     1108   ENDIF  
     1109 
     1110         ! if single ocean point put as land 
     1111         zmask=1 
     1112         WHERE (mbathy .EQ. 0) zmask=0 
     1113         DO jj = 2, jpjm1 
     1114            DO ji = 2, jpim1 
     1115               ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
     1116               IF (ibtest .LE. 1) THEN 
     1117                  bathy(ji,jj)=0._wp 
     1118                  mbathy(ji,jj)=0 
     1119                  icedep(ji,jj)=0._wp 
     1120                  micedep(ji,jj)=0 
     1121               END IF 
     1122            END DO 
     1123         END DO 
     1124   IF( lk_mpp ) THEN  
     1125      zbathy(:,:) = FLOAT( micedep(:,:) )  
     1126      CALL lbc_lnk( zbathy, 'T', 1. )  
     1127      micedep(:,:) = INT( zbathy(:,:) )  
     1128      CALL lbc_lnk( icedep, 'T', 1. )  
     1129           CALL lbc_lnk( bathy, 'T', 1. ) 
     1130           zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1131           CALL lbc_lnk( zbathy, 'T', 1. ) 
     1132           mbathy(:,:) = INT( zbathy(:,:) ) 
     1133   ENDIF  
     1134 
     1135         ! if single point on isf coast line 
     1136         DO jk = 1, jpk 
     1137         WHERE (micedep==0) micedep=jpk 
     1138         zmask=0 
     1139         WHERE (micedep .LE. jk) zmask=1 
     1140         DO jj = 2, jpjm1 
     1141            DO ji = 2, jpim1 
     1142               IF (micedep(ji,jj) .EQ. jk) THEN 
     1143               ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 
     1144               IF (ibtest .LE. 1) THEN 
     1145                     icedep(ji,jj)=gdepw_1d(jk+1) ; micedep(ji,jj)=jk+1 
     1146                     IF (micedep(ji,jj) .GT. mbathy(ji,jj)) micedep(ji,jj) = jpk 
     1147                !     bathy(ji,jj)=0.  ; mbathy(ji,jj)=0 
     1148                  !END IF 
     1149               END IF 
     1150               END IF 
     1151            END DO 
     1152         END DO 
     1153         WHERE (micedep==jpk)  
     1154             micedep=0 ; icedep=0._wp ; mbathy=0 ; bathy=0._wp 
     1155         END WHERE 
     1156 
     1157   IF( lk_mpp ) THEN  
     1158      zbathy(:,:) = FLOAT( micedep(:,:) )  
     1159      CALL lbc_lnk( zbathy, 'T', 1. )  
     1160      micedep(:,:) = INT( zbathy(:,:) )  
     1161      CALL lbc_lnk( icedep, 'T', 1. )  
     1162           CALL lbc_lnk( bathy, 'T', 1. ) 
     1163           zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1164           CALL lbc_lnk( zbathy, 'T', 1. ) 
     1165           mbathy(:,:) = INT( zbathy(:,:) ) 
     1166   ENDIF  
     1167         END DO 
     1168 
     1169! fill hole in ice shelf 
     1170         WHERE (micedep==0) micedep=jpk 
     1171        zmicedep(:,:)=micedep(:,:) 
     1172        zmbathy (:,:)=mbathy (:,:) 
     1173         DO jj = 2, jpjm1 
     1174            DO ji = 2, jpim1 
     1175                  ibtest=MIN(zmicedep(ji-1,jj), zmicedep(ji+1,jj), zmicedep(ji,jj-1), zmicedep(ji,jj+1)) 
     1176                  IF( ibtest > zmicedep(ji,jj)) THEN 
     1177                     micedep(ji,jj) = ibtest 
     1178                     icedep(ji,jj)  = gdepw_1d(ibtest) 
     1179                  ENDIF 
     1180            ENDDO 
     1181         ENDDO 
     1182         WHERE (micedep==jpk)  
     1183             micedep=0 ; icedep=0. ; mbathy=0 ; bathy=0 
     1184         END WHERE 
     1185         IF( lk_mpp ) THEN  
     1186            zbathy(:,:) = FLOAT( micedep(:,:) )  
     1187            CALL lbc_lnk( zbathy, 'T', 1. )  
     1188            micedep(:,:) = INT( zbathy(:,:) )  
     1189            CALL lbc_lnk( icedep, 'T', 1. )  
     1190            CALL lbc_lnk( bathy, 'T', 1. ) 
     1191            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1192            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1193            mbathy(:,:) = INT( zbathy(:,:) ) 
     1194        ENDIF  
     1195 
     1196! fill hole in bathymetry 
     1197        zmicedep(:,:)=micedep(:,:) 
     1198        zmbathy (:,:)=mbathy (:,:) 
     1199         DO jj = 2, jpjm1 
     1200            DO ji = 2, jpim1 
     1201               ibtest = MAX(  zmbathy(ji-1,jj), zmbathy(ji+1,jj),   & 
     1202                  &           zmbathy(ji,jj-1), zmbathy(ji,jj+1)  ) 
     1203               IF( ibtest < zmbathy(ji,jj) ) THEN 
     1204                  mbathy(ji,jj) = ibtest 
     1205                  bathy(ji,jj)  = gdepw_1d(ibtest+1)  
     1206               ENDIF 
     1207            END DO 
     1208         END DO 
     1209         IF( lk_mpp ) THEN  
     1210            zbathy(:,:) = FLOAT( micedep(:,:) )  
     1211            CALL lbc_lnk( zbathy, 'T', 1. )  
     1212            micedep(:,:) = INT( zbathy(:,:) )  
     1213            CALL lbc_lnk( icedep, 'T', 1. )  
     1214            CALL lbc_lnk( bathy, 'T', 1. ) 
     1215            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1216            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1217            mbathy(:,:) = INT( zbathy(:,:) ) 
     1218        ENDIF  
     1219        ! remove pool of water stuck between ice shelf and bathymetry 
     1220        DO jk = 1, jpk 
     1221        WHERE (micedep==0) micedep=jpk 
     1222        zmicedep(:,:)=micedep(:,:) 
     1223        zmbathy (:,:)=mbathy (:,:) 
     1224        DO jj = 2, jpjm1 
     1225            DO ji = 2, jpim1 
     1226               IF( jk .GE. zmicedep(ji,jj) .AND. jk .LE. zmbathy(ji,jj) ) THEN 
     1227               IF( (jk > zmbathy(ji,jj+1) .OR. jk < zmicedep(ji,jj+1)) .AND. & 
     1228                 & (jk > zmbathy(ji,jj-1) .OR. jk < zmicedep(ji,jj-1)) .AND. & 
     1229                 & (jk > zmbathy(ji+1,jj) .OR. jk < zmicedep(ji+1,jj)) .AND. & 
     1230                 & (jk > zmbathy(ji-1,jj) .OR. jk < zmicedep(ji-1,jj))       ) THEN 
     1231                  mbathy(ji,jj) = 0 ; micedep(ji,jj) = jpk ; icedep(ji,jj) = 0._wp ; bathy(ji,jj) = 0._wp 
     1232               ENDIF 
     1233               ENDIF 
     1234            ENDDO 
     1235        ENDDO 
     1236        WHERE (micedep==jpk) micedep=0  
     1237        IF( lk_mpp ) THEN  
     1238            zbathy(:,:) = FLOAT( micedep(:,:) )  
     1239            CALL lbc_lnk( zbathy, 'T', 1. )  
     1240            micedep(:,:) = INT( zbathy(:,:) )  
     1241            CALL lbc_lnk( icedep, 'T', 1. )  
     1242            CALL lbc_lnk( bathy, 'T', 1. ) 
     1243            zbathy(:,:) = FLOAT( mbathy(:,:) ) 
     1244            CALL lbc_lnk( zbathy, 'T', 1. ) 
     1245            mbathy(:,:) = INT( zbathy(:,:) ) 
     1246        ENDIF  
     1247        ENDDO 
     1248      END DO  
     1249 
     1250      WHERE (mbathy(:,:) < micedep(:,:)) 
     1251         micedep(:,:) = 0 
     1252         icedep(:,:)  = 0._wp 
     1253         mbathy(:,:)  = 0 
     1254         bathy(:,:)   = 0._wp 
     1255      END WHERE 
     1256 
     1257 
     1258      IF( icompt == 0 ) THEN  
     1259         IF(lwp) WRITE(numout,*)'     no points with ice shelf too close to bathymetry'  
     1260      ELSE  
     1261         IF(lwp) WRITE(numout,*)'    ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'  
     1262      ENDIF  
    9081263 
    9091264      ! Scale factors and depth at T- and W-points 
    … …  
    9381293!gm Bug?  check the gdepw_1d 
    9391294                  !       ... on ik 
    940                   gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0   (ji,jj,ik+1) - gdepw_1d(ik) )   & 
    941                      &                           * ((gdept_1d(      ik  ) - gdepw_1d(ik) )   & 
    942                      &                           / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )) 
     1295                  gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   & 
     1296                     &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   & 
     1297                     &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
    9431298                  e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   &  
    9441299                     &                          / ( gdepw_1d(      ik+1) - gdepw_1d(ik) )  
    … …  
    9731328         END DO 
    9741329      END DO 
     1330      ! 
     1331      ! (ISF) Definition of e3t, u, v, w for ISF case 
     1332      DO jj = 1, jpj  
     1333         DO ji = 1, jpi  
     1334            ik = micedep(ji,jj)  
     1335            IF( ik > 1 ) THEN               ! ice shelf point only  
     1336                IF( icedep(ji,jj) < gdepw_1d(ik) )  icedep(ji,jj)= gdepw_1d(ik)  
     1337                gdepw_0(ji,jj,ik) = icedep(ji,jj)  
     1338!gm Bug?  check the gdepw_0  
     1339                !       ... on ik  
     1340                gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) )   &  
     1341                   &                               * ( gdepw_1d(ik+1) - gdept_1d(ik)      )   &  
     1342                   &                               / ( gdepw_1d(ik+1) - gdepw_1d(ik)      )  
     1343                e3t_0  (ji,jj,ik  ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)  
     1344                e3w_0  (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 
     1345 
     1346                IF( ik + 1 == mbathy(ji,jj) ) THEN               ! ice shelf point only  
     1347                   e3w_0  (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)  
     1348                ENDIF  
     1349                !       ... on ik / ik-1  
     1350                e3w_0  (ji,jj,ik  ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))  
     1351                e3t_0  (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 
     1352! The next line isn't required and doesn't affect results - included for consistency with bathymetry code  
     1353                gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 
     1354             ENDIF  
     1355          END DO  
     1356       END DO  
     1357      !  
     1358      it = 0  
     1359      DO jj = 1, jpj  
     1360         DO ji = 1, jpi  
     1361            ik = micedep(ji,jj)  
     1362            IF( ik > 1 ) THEN               ! ice shelf point only  
     1363               e3tp (ji,jj) = e3t_0(ji,jj,ik  )  
     1364               e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )  
     1365               ! test  
     1366               zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik  )  
     1367               IF( zdiff <= 0. .AND. lwp ) THEN   
     1368                  it = it + 1  
     1369                  WRITE(numout,*) ' it      = ', it, ' ik      = ', ik, ' (i,j) = ', ji, jj  
     1370                  WRITE(numout,*) ' icedep = ', icedep(ji,jj)  
     1371                  WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff  
     1372                  WRITE(numout,*) ' e3tp  = ', e3tp(ji,jj), ' e3wp  = ', e3wp(ji,jj)  
     1373               ENDIF  
     1374            ENDIF  
     1375         END DO  
     1376      END DO  
     1377      ! END (ISF) 
    9751378 
    9761379      ! Scale factors and depth at U-, V-, UW and VW-points 
    … …  
    9911394         END DO 
    9921395      END DO 
     1396      ! (ISF) define e3uw 
     1397      DO jk = 2,jpk                           
     1398         DO jj = 1, jpjm1  
     1399            DO ji = 1, fs_jpim1   ! vector opt.  
     1400! (ISF)** NEEDS CHANGING TO SECOND OPTION FOR ICE SHELF BUT WILL CHANGE RESULTS WITHOUT ICE (DIFFER AT THE 1e-13 LEVEL)  
     1401               e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj  ,jk) ) - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj  ,jk-1) ) 
     1402               e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji  ,jj+1,jk) ) - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji  ,jj+1,jk-1) ) 
     1403            END DO  
     1404         END DO  
     1405      END DO 
     1406      !End (ISF) 
     1407       
    9931408      CALL lbc_lnk( e3u_0 , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw_0, 'U', 1._wp )   ! lateral boundary conditions 
    9941409      CALL lbc_lnk( e3v_0 , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 
    … …  
    10321447      
    10331448      ! Compute gdep3w_0 (vertical sum of e3w) 
    1034       gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 
    1035       DO jk = 2, jpk 
    1036          gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)  
    1037       END DO 
    1038          
     1449      WHERE (micedep == 0) micedep = 1 
     1450      DO jj = 1,jpj 
     1451         DO ji = 1,jpi 
     1452            gdep3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 
     1453            DO jk = 2, micedep(ji,jj) 
     1454               gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1455            END DO 
     1456            IF (micedep(ji,jj) .GE. 2) gdep3w_0(ji,jj,micedep(ji,jj)) = icedep(ji,jj) + 0.5_wp * e3w_0(ji,jj,micedep(ji,jj)) 
     1457            DO jk = micedep(ji,jj) + 1, jpk 
     1458               gdep3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)  
     1459            END DO 
     1460        END DO 
     1461      END DO 
    10391462      !                                               ! ================= ! 
    10401463      IF(lwp .AND. ll_print) THEN                     !   Control print   ! 
    … …  
    10661489      ! 
    10671490      CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 
     1491      CALL wrk_dealloc( jpi, jpj, zmask, zbathy ) 
     1492      CALL wrk_dealloc( jpi, jpj, zmicedep, zmbathy ) 
    10681493      ! 
    10691494      IF( nn_timing == 1 )  CALL timing_stop('zgr_zps') 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90

    r4624 r4666  
    263263                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik-1,jp_sal) 
    264264                  ENDIF 
     265                  ik = mikt(ji,jj) 
     266                  IF( ik > 1 ) THEN 
     267                     zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )  
     268                     ptsd(ji,jj,ik,jp_tem) = (1.-zl) * ptsd(ji,jj,ik,jp_tem) + zl * ptsd(ji,jj,ik+1,jp_tem) 
     269                     ptsd(ji,jj,ik,jp_sal) = (1.-zl) * ptsd(ji,jj,ik,jp_sal) + zl * ptsd(ji,jj,ik+1,jp_sal) 
     270                  END IF 
    265271               END DO 
    266272            END DO 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90

    r4370 r4666  
    7676      ! 
    7777 
    78       IF(lwp) WRITE(numout,*) 
     78      IF(lwp) WRITE(numout,*) ' ' 
    7979      IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers' 
    8080      IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 
    … …  
    110110         ELSEIF( cp_cfg == 'gyre' ) THEN          
    111111            CALL istate_gyre                     ! GYRE  configuration : start from pre-defined T-S fields 
     112        ELSEIF( cp_cfg == 'isomip' .OR. cp_cfg == 'isomip2') THEN 
     113            IF(lwp) WRITE(numout,*) 'Initialization of T+S for ISOMIP domain'  
     114            tsn(:,:,:,jp_tem)=-1.9*tmask(:,:,:)          ! ISOMIP configuration : start from constant T+S fields  
     115            tsn(:,:,:,jp_sal)=34.4*tmask(:,:,:) 
     116            tsb(:,:,:,:)=tsn(:,:,:,:)   
    112117         ELSE                                    ! Initial T-S, U-V fields read in files 
    113118            IF ( ln_tsd_init ) THEN              ! read 3D T and S data at nit000 
    … …  
    129134         CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )        ! before potential and in situ densities 
    130135#if ! defined key_c1d 
    131          IF( ln_zps )   CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv,  & ! zps: before hor. gradient 
    132             &                                       rhd, gru , grv  )   ! of t,s,rd at ocean bottom 
     136         IF( ln_zps )    CALL zps_hde( nit000, jpts, tsb, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
     137            &                                      rhd, gru , grv, aru, arv, gzu, gzv, ge3ru, ge3rv,  &             ! 
     138            &                                      gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
    133139#endif 
    134140         !    

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90

    r3625 r4666  
    5050   REAL(wp), PUBLIC ::   rau0     = 1026._wp         !: volumic mass of reference     [kg/m3] 
    5151#else 
    52    REAL(wp), PUBLIC ::   rau0     = 1035._wp         !: volumic mass of reference     [kg/m3] 
     52   REAL(wp), PUBLIC ::   rau0     = 1028.4_wp       !: volumic mass of reference     [kg/m3] 
    5353#endif 
    5454   REAL(wp), PUBLIC ::   r1_rau0                     !: = 1. / rau0                   [m3/kg] 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90

    r4147 r4666  
    2525   USE oce             ! ocean dynamics and tracers 
    2626   USE dom_oce         ! ocean space and time domain 
    27    USE sbc_oce, ONLY : ln_rnf  ! surface boundary condition: ocean 
     27   USE sbc_oce, ONLY : ln_rnf, nn_isf ! surface boundary condition: ocean 
    2828   USE sbcrnf          ! river runoff  
     29   USE sbcisf          ! ice shelf  
    2930   USE cla             ! cross land advection             (cla_div routine) 
    3031   USE in_out_manager  ! I/O manager 
    … …  
    6667      !!         - compute the now divergence given by : 
    6768      !!         hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 
    68       !!      correct hdiv with runoff inflow (div_rnf) and cross land flow (div_cla)  
     69      !!      correct hdiv with runoff inflow (div_rnf), ice shelf melting (div_isf) 
     70      !!      and cross land flow (div_cla)  
    6971      !!              II. vorticity : 
    7072      !!         - save the curl computed at the previous time-step 
    … …  
    225227      !                                                ! =============== 
    226228 
    227       IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field) 
    228       IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 )   CALL cla_div    ( kt )             ! Cross Land Advection (Update Hor. divergence) 
     229      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs   (update hdivn field) 
     230      IF( ln_divisf .AND. (nn_isf /= 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field) 
     231      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (Update Hor. divergence) 
    229232       
    230233      ! 4. Lateral boundary conditions on hdivn and rotn 
    … …  
    323326      !                                                ! =============== 
    324327 
    325       IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )          ! runoffs (update hdivn field) 
     328      IF( ln_rnf      )   CALL sbc_rnf_div( hdivn )                            ! runoffs (update hdivn field) 
     329      IF( ln_divisf .AND. (nn_isf .GT. 0) )   CALL sbc_isf_div( hdivn )          ! ice shelf (update hdivn field) 
    326330      IF( nn_cla == 1 )   CALL cla_div    ( kt )             ! Cross Land Advection (update hdivn field) 
    327331      ! 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90

    r3294 r4666  
    8282              ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) 
    8383              va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) 
     84               
     85              ! (ISF) stability criteria for top friction 
     86              ikbu = miku(ji,jj)          ! first wet ocean u- & v-levels 
     87              ikbv = mikv(ji,jj) 
     88              ! 
     89              ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
     90              ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  tfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) & 
     91                 &             * (1.-umask(ji,jj,1)) 
     92              va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  tfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) & 
     93                 &             * (1.-vmask(ji,jj,1)) 
     94              ! (ISF) 
     95               
    8496           END DO 
    8597        END DO 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r4624 r4666  
    3737   USE lbclnk          ! lateral boundary condition 
    3838   USE lib_mpp         ! MPP library 
     39   USE eosbn2          ! compute density 
    3940   USE wrk_nemo        ! Memory Allocation 
    4041   USE timing          ! Timing 
    … …  
    135136      READ  ( numnam_cfg, namdyn_hpg, IOSTAT = ios, ERR = 902 ) 
    136137902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist', lwp ) 
    137       IF(lwm) WRITE ( numond, namdyn_hpg ) 
     138      WRITE ( numond, namdyn_hpg ) 
    138139      ! 
    139140      IF(lwp) THEN                   ! Control print 
    … …  
    368369      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
    369370      !! 
    370       INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
    371       REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars 
    372       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj 
    373       !!---------------------------------------------------------------------- 
    374       ! 
    375       CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 
    376       ! 
    377       IF( kt == nit000 ) THEN 
     371      INTEGER  ::   ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j                 ! dummy loop indices 
     372      REAL(wp) ::   zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1             ! temporary scalars 
     373      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  zhpi, zhpj, zrhd 
     374      REAL(wp), POINTER, DIMENSION(:,:,:)   ::  ztstop 
     375      REAL(wp), POINTER, DIMENSION(:,:)     ::  ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj 
     376      !!---------------------------------------------------------------------- 
     377      ! 
     378      CALL wrk_alloc( jpi,jpj, 2, ztstop)  
     379      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 
     380      CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj)  
     381      ! 
     382     IF( kt == nit000 ) THEN 
    378383         IF(lwp) WRITE(numout,*) 
    379384         IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' 
    … …  
    384389      zcoef0 = - grav * 0.5_wp 
    385390      ! To use density and not density anomaly 
    386       IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume 
    387       ELSE                 ;     znad = 0._wp         ! Fixed volume 
    388       ENDIF 
    389  
    390       ! Surface value 
     391!      IF ( lk_vvl ) THEN   ;     znad = 1._wp          ! Variable volume 
     392!      ELSE                 ;     znad = 0._wp         ! Fixed volume 
     393!      ENDIF 
     394      znad=1._wp 
     395      ! iniitialised to 0. zhpi zhpi  
     396      zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 
     397 
     398      ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 
     399      zrhd = rhd 
     400      DO jk = 1, jpk 
     401           zdept(:,:)=gdept_1d(jk) 
     402           CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk)) 
     403      END DO 
     404      WHERE ( tmask(:,:,:) == 1._wp) 
     405        rhd(:,:,:) = zrhd(:,:,:) 
     406      END WHERE 
     407       
     408 
     409      CALL eos(ztstop,icedep,zrhdtop_isf) 
     410 
     411      DO ji=1,jpi 
     412        DO jj=1,jpj 
     413          ikt=mikt(ji,jj) 
     414          ztstop(ji,jj,1)=tsn(ji,jj,ikt,1) 
     415          ztstop(ji,jj,2)=tsn(ji,jj,ikt,2) 
     416        END DO 
     417      END DO 
     418      CALL eos(ztstop,icedep,zrhdtop_oce) 
     419      ! 
     420      ! Surface value + ice shelf gradient 
     421      ! compute pressure (used to compute hpgi/j for all the level from 1 to miku/v) 
     422      ziceload = 0._wp 
     423      DO jj = 1, jpj 
     424         DO ji = 1, jpi   ! vector opt. 
     425            ikt=mikt(ji,jj) 
     426            ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 
     427            DO jk=2,ikt-1 
     428               ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 
     429                  &                              * (1._wp - tmask(ji,jj,jk)) 
     430            END DO 
     431            IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) & 
     432                               &                              * ( fsdepw(ji,jj,ikt) - gdept_1d(ikt-1) ) 
     433         END DO 
     434      END DO 
     435      ! compute zp from first level to first wet cell (blue and purple area) 
    391436      DO jj = 2, jpjm1 
    392437         DO ji = fs_2, fs_jpim1   ! vector opt. 
    393             ! hydrostatic pressure gradient along s-surfaces 
    394             zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   & 
    395                &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) ) 
    396             zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   & 
    397                &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) ) 
     438            ikt=mikt(ji,jj) ; iktp1i=mikt(ji+1,jj); iktp1j=mikt(ji,jj+1) 
     439            ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 
     440            ! we assume ISF is in isostatic equilibrium with a density equal to the reference density 
     441            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj  ,iktp1i)                                    & 
     442               &                                  * ( 2._wp * znad + rhd(ji+1,jj  ,iktp1i) + zrhdtop_oce(ji+1,jj  ) )   & 
     443               &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    & 
     444               &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   & 
     445               &                                  + ( ziceload(ji+1,jj) - ziceload(ji,jj))                              )  
     446            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji  ,jj+1,iktp1j)                                    & 
     447               &                                  * ( 2._wp * znad + rhd(ji  ,jj+1,iktp1j) + zrhdtop_oce(ji  ,jj+1) )   & 
     448               &                                  - 0.5_wp * fse3w(ji  ,jj  ,ikt   )                                    &  
     449               &                                  * ( 2._wp * znad + rhd(ji  ,jj  ,ikt   ) + zrhdtop_oce(ji  ,jj  ) )   & 
     450               &                                  + ( ziceload(ji,jj+1) - ziceload(ji,jj) )                             )  
    398451            ! s-coordinate pressure gradient correction 
    399452            zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   & 
    … …  
    402455               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 
    403456            ! add to the general momentum trend 
    404             ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 
    405             va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 
     457            ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 
     458            va(ji,jj,1) = va(ji,jj,1) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 
     459         END DO 
     460      END DO 
     461      ! 
     462      DO jj = 2, jpjm1 
     463         DO ji = fs_2, fs_jpim1   ! vector opt. 
     464            iku = miku(ji,jj) ;  
     465            zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp 
     466            ze3wu  = fse3w(ji+1,jj,iku+1)-fse3w(ji,jj,iku+1)  
     467            ! u direction 
     468            IF ( iku .GT. 1 ) THEN 
     469               ! case iku 
     470               zhpi(ji,jj,iku)   =  zcoef0 / e1u(ji,jj) * ze3wu                                            & 
     471                      &                                 * ( rhd    (ji+1,jj,iku) + rhd   (ji,jj,iku)       & 
     472                      &                                   + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad ) 
     473  
     474               IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj)  = zhpi(ji,jj,iku) 
     475               zuap              = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj) 
     476               ! zhpi will be added in interior loop and zuapint will be removed in the interior loop 
     477               ua(ji,jj,iku)     = ua(ji,jj,iku) + zuap 
     478 
     479               ! zhpi will be added in interior loop and zuapint will be removed in the interior loop 
     480               ! case iku + 1 (remove the zphi term added in the interior loop) 
     481               zhpiint        =  zcoef0 / e1u(ji,jj)                                                             &     
     482                  &           * (  fse3w(ji+1,jj  ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad)                          & 
     483                  &                                         + (rhd(ji+1,jj,iku  ) + znad) ) * tmask(ji+1,jj,iku)   & 
     484                  &              - fse3w(ji  ,jj  ,iku+1) * ( (rhd(ji  ,jj,iku+1) + znad)                          & 
     485                  &                                         + (rhd(ji  ,jj,iku  ) + znad) ) * tmask(ji  ,jj,iku)   ) 
     486               zhpi(ji,jj,iku+1) =  zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint  
     487            END IF 
     488            ! v direction 
     489            ikv = mikv(ji,jj) 
     490            !ze3wv  = 0.5 * (e3w_0(ikv+1) / e3t_0(ikv)) * ( fse3t(ji,jj+1,ikv)-fse3t(ji,jj,ikv) ) 
     491            ze3wv  = fse3w(ji,jj+1,ikv+1)-fse3w(ji,jj,ikv+1)  
     492            IF ( ikv .GT. 1 ) THEN 
     493               ! case ikv 
     494         !      ze3wv             =  - (fsde3w(ji,jj+1,ikv) - fsde3w(ji,jj,ikv) ) 
     495               zhpj(ji,jj,ikv)   =  zcoef0 / e2v(ji,jj) * ze3wv                                            & 
     496                     &                                  * ( rhd(ji,jj+1,ikv) + rhd   (ji,jj,ikv)           & 
     497                     &                                    + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad ) 
     498 
     499               zvap              = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj) 
     500               ! zhpi will be added in interior loop and zvapint will be removed in the interior loop 
     501               va(ji,jj,ikv)      = va(ji,jj,ikv) + zvap 
     502               IF (mbkv(ji,jj) == ikv + 1)  zpshpj(ji,jj)  =  zhpj(ji,jj,ikv)  
     503               ! zhpi will be added in interior loop and zvapint will be removed in the interior loop 
     504               ! case ikv + 1 (remove the zphj term added in the interior loop) 
     505               zhpjint        =  zcoef0 / e2v(ji,jj)                                                            & 
     506                  &           * (  fse3w(ji  ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad)                         & 
     507                  &                                       + (rhd(ji,jj+1,ikv  ) + znad) ) * tmask(ji,jj+1,ikv)  & 
     508                  &              - fse3w(ji  ,jj  ,ikv+1) * ( (rhd(ji,jj  ,ikv+1) + znad)                         & 
     509                  &                                       + (rhd(ji,jj  ,ikv  ) + znad) ) * tmask(ji,jj  ,ikv)  ) 
     510               zhpj(ji,jj,ikv+1) =  zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint 
     511            END IF 
    406512         END DO 
    407513      END DO 
    408514 
    409515      ! interior value (2=<jk=<jpkm1) 
    410       DO jk = 2, jpkm1 
    411          DO jj = 2, jpjm1 
    412             DO ji = fs_2, fs_jpim1   ! vector opt. 
     516      DO jj = 2, jpjm1 
     517         DO ji = fs_2, fs_jpim1   ! vector opt. 
     518            iku=miku(ji,jj); ikv=mikv(ji,jj) 
     519            DO jk = 2, jpkm1 
    413520               ! hydrostatic pressure gradient along s-surfaces 
    414                zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   & 
    415                   &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   & 
    416                   &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  ) 
    417                zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   & 
    418                   &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   & 
    419                   &              - fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  ) 
    420521               ! s-coordinate pressure gradient correction 
    421                zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    422                   &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 
    423                zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
    424                   &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 
     522               zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1)                                                              & 
     523                  &                            + zcoef0 / e1u(ji,jj)                                                           & 
     524                  &                                     * ( fse3w(ji+1,jj  ,jk) * ( (rhd(ji+1,jj,jk  ) + znad)                 & 
     525                  &                                                     + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1)   & 
     526                  &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji  ,jj,jk  ) + znad)                 & 
     527                  &                                                     + (rhd(ji  ,jj,jk-1) + znad) ) * tmask(ji  ,jj,jk-1)   )  
     528               ! corrective term                
     529               zuap = - zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
     530                  &            * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1) 
     531               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 
     532 
     533               zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1)                                                              & 
     534                  &                            + zcoef0 / e2v(ji,jj)                                                           & 
     535                  &                                     * ( fse3w(ji  ,jj+1,jk) * ( (rhd(ji,jj+1,jk  ) + znad)                 & 
     536                  &                                                     + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1)   & 
     537                  &                                       - fse3w(ji  ,jj  ,jk) * ( (rhd(ji,jj  ,jk  ) + znad)                 & 
     538                  &                                                     + (rhd(ji,jj  ,jk-1) + znad) ) * tmask(ji,jj  ,jk-1)   ) 
     539               ! 
     540               zvap = - zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   & 
     541                  &            * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1) 
    425542               ! add to the general momentum trend 
    426                ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 
    427                va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 
     543               va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 
    428544            END DO 
    429545         END DO 
    430546      END DO 
    431       ! 
    432       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 
     547      
     548      ! partial steps correction at the last level  (use gru & grv computed in zpshde.F90) 
     549# if defined key_vectopt_loop 
     550         jj = 1 
     551         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
     552# else 
     553      DO jj = 2, jpjm1 
     554         DO ji = 2, jpim1 
     555# endif 
     556            iku = mbku(ji,jj) 
     557            ikv = mbkv(ji,jj) 
     558 
     559            IF (iku .GT. 1) THEN 
     560               ! remove old value (interior case) 
     561               zuap            = -zcoef0 * ( rhd   (ji+1,jj  ,iku) + rhd   (ji,jj,iku) + 2._wp * znad )   & 
     562                     &                   * ( fsde3w(ji+1,jj  ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj) 
     563               ua(ji,jj,iku)   = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 
     564               ! put new value 
     565               zuap            = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj)  / e1u(ji,jj) 
     566               zhpi(ji,jj,iku) =  zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj)  
     567               ua(ji,jj,iku)   =  ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 
     568            END IF 
     569            ! v direction 
     570            IF (ikv .GT. 1) THEN 
     571               ! remove old value (interior case) 
     572               zvap            = -zcoef0 * ( rhd   (ji  ,jj+1,ikv) + rhd   (ji,jj,ikv) + 2._wp * znad )   & 
     573                     &                   * ( fsde3w(ji  ,jj+1,ikv) - fsde3w(ji,jj,ikv) )   / e2v(ji,jj) 
     574               va(ji,jj,ikv)   = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 
     575               ! put new value 
     576               zvap            = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj)     / e2v(ji,jj) 
     577               zhpj(ji,jj,ikv) =  zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj)    
     578               va(ji,jj,ikv)   =  va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 
     579            END IF 
     580# if ! defined key_vectopt_loop 
     581         END DO 
     582# endif 
     583      END DO 
     584 
     585 
     586      rhd = zrhd 
     587      ! 
     588      CALL wrk_dealloc( jpi,jpj,2, ztstop) 
     589      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 
     590      CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 
    433591      ! 
    434592   END SUBROUTINE hpg_sco 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90

    r4370 r4666  
    327327            hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 
    328328         END DO 
    329          hur_b(:,:) = umask(:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) ) 
    330          hvr_b(:,:) = vmask(:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) ) 
     329         hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) ) 
     330         hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) ) 
    331331      ENDIF 
    332332      ! 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90

    r3294 r4666  
    9595      DO jj = 2, jpjm1              ! Surface and bottom values set to zero 
    9696         DO ji = fs_2, fs_jpim1           ! vector opt. 
    97             zwuw(ji,jj, 1 ) = 0.e0 
    98             zwvw(ji,jj, 1 ) = 0.e0 
     97            zwuw(ji,jj, 1:miku(ji,jj) ) = 0.e0 
     98            zwvw(ji,jj, 1:mikv(ji,jj) ) = 0.e0 
    9999            zwuw(ji,jj,jpk) = 0.e0 
    100100            zwvw(ji,jj,jpk) = 0.e0 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r4370 r4666  
    112112               avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 
    113113               avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 
     114               ikbu = miku(ji,jj)       ! ocean top level at u- and v-points  
     115               ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points) 
     116               avmu(ji,jj,ikbu-1) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 
     117               avmv(ji,jj,ikbv-1) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 
    114118            END DO 
    115119         END DO 
    … …  
    139143            va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk) 
    140144         ENDDO 
    141          ! Add bottom stress due to barotropic component only: 
     145         ! Add bottom/top stress due to barotropic component only: 
    142146         DO jj = 2, jpjm1         
    143147            DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    148152               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua 
    149153               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va 
     154               ikbu = miku(ji,jj)         ! ocean bottom level at u- and v-points  
     155               ikbv = mikv(ji,jj)         ! (deepest ocean u- and v-points) 
     156               ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl   * fse3u_a(ji,jj,ikbu) 
     157               ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl   * fse3v_a(ji,jj,ikbv) 
     158               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua 
     159               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va 
    150160            END DO 
    151161         END DO 
    … …  
    174184      DO jj = 2, jpjm1        ! Surface boundary conditions 
    175185         DO ji = fs_2, fs_jpim1   ! vector opt. 
    176             zwi(ji,jj,1) = 0._wp 
    177             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
     186            zwi(ji,jj,miku(ji,jj)) = 0._wp 
     187            zwd(ji,jj,miku(ji,jj)) = 1._wp - zws(ji,jj,miku(ji,jj)) 
    178188         END DO 
    179189      END DO 
    … …  
    194204      !----------------------------------------------------------------------- 
    195205      ! 
    196       DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    197          DO jj = 2, jpjm1    
    198             DO ji = fs_2, fs_jpim1   ! vector opt. 
     206      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
     207      DO jj = 2, jpjm1    
     208         DO ji = fs_2, fs_jpim1   ! vector opt. 
     209            DO jk = miku(ji,jj)+1, jpkm1 
    199210               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    200211            END DO 
    … …  
    204215      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    205216         DO ji = fs_2, fs_jpim1   ! vector opt. 
    206             ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl   * fse3u_a(ji,jj,1)  
     217            ze3ua =  ( 1._wp - r_vvl ) * fse3u_n(ji,jj,miku(ji,jj)) + r_vvl   * fse3u_a(ji,jj,miku(ji,jj))  
    207218#if defined key_dynspg_ts 
    208             ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     219            ua(ji,jj,miku(ji,jj)) = ua(ji,jj,miku(ji,jj)) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    209220               &                                      / ( ze3ua * rau0 )  
    210221#else 
    211             ua(ji,jj,1) = ub(ji,jj,1) + p2dt *(ua(ji,jj,1) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    212                &                                                     / ( fse3u(ji,jj,1) * rau0     ) )  
    213 #endif 
    214          END DO 
    215       END DO 
    216       DO jk = 2, jpkm1 
    217          DO jj = 2, jpjm1    
    218             DO ji = fs_2, fs_jpim1   ! vector opt. 
     222            ua(ji,jj,miku(ji,jj)) = ub(ji,jj,miku(ji,jj)) + p2dt *(ua(ji,jj,miku(ji,jj)) +  0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
     223               &                                                     / ( fse3u(ji,jj,miku(ji,jj)) * rau0     ) )  
     224#endif 
     225            DO jk = miku(ji,jj)+1, jpkm1 
    219226#if defined key_dynspg_ts 
    220227               zrhs = ua(ji,jj,jk)   ! zrhs=right hand side 
    … …  
    230237         DO ji = fs_2, fs_jpim1   ! vector opt. 
    231238            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    232          END DO 
    233       END DO 
    234       DO jk = jpk-2, 1, -1 
    235          DO jj = 2, jpjm1    
    236             DO ji = fs_2, fs_jpim1   ! vector opt. 
     239            DO jk = jpk-2, miku(ji,jj), -1 
    237240               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    238241            END DO 
    … …  
    272275      DO jj = 2, jpjm1        ! Surface boundary conditions 
    273276         DO ji = fs_2, fs_jpim1   ! vector opt. 
    274             zwi(ji,jj,1) = 0._wp 
    275             zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
     277            zwi(ji,jj,mikv(ji,jj)) = 0._wp 
     278            zwd(ji,jj,mikv(ji,jj)) = 1._wp - zws(ji,jj,mikv(ji,jj)) 
    276279         END DO 
    277280      END DO 
    … …  
    292295      !----------------------------------------------------------------------- 
    293296      ! 
    294       DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    295          DO jj = 2, jpjm1    
    296             DO ji = fs_2, fs_jpim1   ! vector opt. 
     297      !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
     298      DO jj = 2, jpjm1    
     299         DO ji = fs_2, fs_jpim1   ! vector opt. 
     300            DO jk = mikv(ji,jj)+1, jpkm1         
    297301               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    298302            END DO 
    … …  
    304308            ze3va =  ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl   * fse3v_a(ji,jj,1)  
    305309#if defined key_dynspg_ts             
    306             va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     310            va(ji,jj,mikv(ji,jj)) = va(ji,jj,mikv(ji,jj)) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    307311               &                                      / ( ze3va * rau0 )  
    308312#else 
    309             va(ji,jj,1) = vb(ji,jj,1) + p2dt *(va(ji,jj,1) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    310                &                                                       / ( fse3v(ji,jj,1) * rau0     )  ) 
    311 #endif 
    312          END DO 
    313       END DO 
    314       DO jk = 2, jpkm1 
    315          DO jj = 2, jpjm1 
    316             DO ji = fs_2, fs_jpim1   ! vector opt. 
     313            va(ji,jj,mikv(ji,jj)) = vb(ji,jj,mikv(ji,jj)) + p2dt *(va(ji,jj,mikv(ji,jj)) +  0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
     314               &                                                       / ( fse3v(ji,jj,mikv(ji,jj)) * rau0     )  ) 
     315#endif 
     316            DO jk = mikv(ji,jj)+1, jpkm1 
    317317#if defined key_dynspg_ts 
    318318               zrhs = va(ji,jj,jk)   ! zrhs=right hand side 
    … …  
    328328         DO ji = fs_2, fs_jpim1   ! vector opt. 
    329329            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    330          END DO 
    331       END DO 
    332       DO jk = jpk-2, 1, -1 
    333          DO jj = 2, jpjm1    
    334             DO ji = fs_2, fs_jpim1   ! vector opt. 
     330            DO jk = jpk-2, mikv(ji,jj), -1 
    335331               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    336332            END DO 
    … …  
    363359            avmu(ji,jj,ikbu+1) = 0.e0 
    364360            avmv(ji,jj,ikbv+1) = 0.e0 
     361            ikbu = miku(ji,jj)         ! ocean top level at u- and v-points  
     362            ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points) 
     363            avmu(ji,jj,ikbu-1) = 0.e0 
     364            avmv(ji,jj,ikbv-1) = 0.e0 
    365365         END DO 
    366366      END DO 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90

    r4334 r4666  
    120120                     CALL iom_rstput( kt, nitrst, numrow, 'hdivb'  , hdivb     ) 
    121121                     CALL iom_rstput( kt, nitrst, numrow, 'sshb'   , sshb      ) 
     122      IF( lk_vvl )  THEN     ! need for ISF 
     123                     CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) ) 
     124                     CALL iom_rstput( kt, nitrst, numrow, 'fse3t  ', fse3t  (:,:,:) ) 
     125                     CALL iom_rstput( kt, nitrst, numrow, 'fse3w  ', fse3w  (:,:,:) ) 
     126                     CALL iom_rstput( kt, nitrst, numrow, 'fsdepw ', fsdepw (:,:,:) ) 
     127      END IF 
    122128                     ! 
    123129                     CALL iom_rstput( kt, nitrst, numrow, 'un'     , un        )     ! now fields 
    … …  
    214220      ENDIF 
    215221      ! 
     222      IF( lk_vvl ) THEN     ! need it for (ISF) 
     223                     CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 
     224                     CALL iom_get( numror, jpdom_autoglo, 'fse3t  ', fse3t  (:,:,:) ) 
     225                     CALL iom_get( numror, jpdom_autoglo, 'fse3w  ', fse3w  (:,:,:) ) 
     226                     CALL iom_get( numror, jpdom_autoglo, 'fsdepw ', fsdepw (:,:,:) ) 
     227      END IF 
    216228      CALL iom_get( numror, jpdom_autoglo, 'un'     , un      )   ! now    fields 
    217229      CALL iom_get( numror, jpdom_autoglo, 'vn'     , vn      ) 
    … …  
    245257         hdivb(:,:,:)   = hdivn(:,:,:) 
    246258         sshb (:,:)     = sshn (:,:) 
     259         IF( lk_vvl ) THEN 
     260            DO jk = 1, jpk 
     261               fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 
     262            END DO 
     263         ENDIF 
    247264      ENDIF 
    248265      ! 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90

    r3294 r4666  
    165165            zross(ji,jj) = MAX( MIN( .4 * zn(ji,jj) / zfw, 40.e3 ), 2.e3 ) 
    166166            ! Compute aeiw by multiplying Ro^2 and T^-1 
    167             aeiw(ji,jj) = zross(ji,jj) * zross(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1) 
     167            aeiw(ji,jj) = zross(ji,jj) * zross(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * lmask(ji,jj) 
    168168         END DO 
    169169      END DO 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

    r4488 r4666  
    143143               DO ji = 1, jpim1 
    144144# endif 
    145                   zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 
    146                   zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 
     145! IF should be useless check zpshde (PM) 
     146               IF ( mbku(ji,jj) > 1 ) zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 
     147               IF ( mbkv(ji,jj) > 1 ) zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 
     148               IF ( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj)  
     149               IF ( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj) 
    147150               END DO 
    148151            END DO 
    … …  
    238241               DO ji = fs_2, fs_jpim1   ! vector opt. 
    239242                  uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk  ) ) * 0.5_wp   & 
    240                      &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5_wp 
     243                     &                            * ( umask(ji,jj  ,jk) + umask(ji,jj  ,jk+1) ) * 0.5_wp   & 
     244                     &                            *   tmask(ji,jj,jk-1) 
    241245                  vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk  ) ) * 0.5_wp   & 
    242246                     &                            * ( vmask(ji  ,jj,jk) + vmask(ji  ,jj,jk+1) ) * 0.5_wp 
     247                     &                            *   tmask(ji,jj,jk-1) 
    243248               END DO 
    244249            END DO 
    … …  
    325330                  zck =   ( umask(ji,jj,jk) + umask(ji-1,jj,jk) )   & 
    326331                     &  * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 
    327                   wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 
    328                   wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 
     332                  wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck * tmask(ji,jj,jk-1) 
     333                  wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck * tmask(ji,jj,jk-1) 
    329334               END DO 
    330335            END DO 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90

    r4306 r4666  
    4141   LOGICAL , PUBLIC ::   ln_apr_dyn     !: Atmospheric pressure forcing used on dynamics (ocean & ice) 
    4242   INTEGER , PUBLIC ::   nn_ice         !: flag for ice in the surface boundary condition (=0/1/2/3) 
     43   INTEGER , PUBLIC ::   nn_isf         !: flag for isf in the surface boundary condition (=0/1/2/3/4)  
    4344   INTEGER , PUBLIC ::   nn_ice_embd    !: flag for levitating/embedding sea-ice in the ocean 
    4445   !                                             !: =0 levitating ice (no mass exchange, concentration/dilution effect) 
    … …  
    102103   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sss_m     !: mean (nn_fsbc time-step) surface sea salinity            [psu] 
    103104   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssh_m     !: mean (nn_fsbc time-step) sea surface height                [m] 
    104    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3t_m     !: mean (nn_fsbc time-step) sea surface height                [m] 
     105   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3t_m     !: mean (nn_fsbc time-step) sea surface layer thickness       [m] 
    105106 
    106107   !! * Substitutions 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90

    r4624 r4666  
    371371      ! ... utau, vtau at U- and V_points, resp. 
    372372      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 
     373      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 
    373374      DO jj = 1, jpjm1 
    374375         DO ji = 1, fs_jpim1 
    375             utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) 
    376             vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) 
     376            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) & 
     377               &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
     378            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) & 
     379               &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 
    377380         END DO 
    378381      END DO 
    … …  
    420423         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &    
    421424         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow) 
    422          &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic  
     425         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) 
    423426      ! 
    424427      CALL iom_put( "qlw_oce",   zqlw )                 ! output downward longwave heat over the ocean 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90

    r4347 r4666  
    1919   USE phycst          ! physical constants 
    2020   USE sbcrnf          ! ocean runoffs 
     21   USE sbcisf          ! ice shelf melting contribution 
    2122   USE sbcssr          ! SS damping terms 
    2223   USE in_out_manager  ! I/O manager 
    … …  
    9899         ! 
    99100         IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN 
    100             z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain 
     101            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) -  snwice_fmass(:,:) ) ) / area   ! sum over the global domain 
    101102            zcoef = z_fwf * rcp 
    102103            emp(:,:) = emp(:,:) - z_fwf  

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r4624 r4666  
    4040   USE sbcssr           ! surface boundary condition: sea surface restoring 
    4141   USE sbcrnf           ! surface boundary condition: runoffs 
     42   USE sbcisf           ! surface boundary condition: ice shelf 
    4243   USE sbcfwb           ! surface boundary condition: freshwater budget 
    4344   USE closea           ! closed sea 
    … …  
    8485      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx,  ln_blk_clio, ln_blk_core, ln_cpl,   & 
    8586         &             ln_blk_mfs, ln_apr_dyn, nn_ice,  nn_ice_embd, ln_dm2dc   , ln_rnf,   & 
    86          &             ln_ssr    , nn_fwb    , ln_cdgw , ln_wave , ln_sdw, nn_lsm, cn_iceflx 
     87         &             ln_ssr    , nn_isf    , nn_fwb,  ln_cdgw    , ln_wave    , ln_sdw, nn_lsm, cn_iceflx 
    8788      INTEGER  ::   ios 
    8889      !!---------------------------------------------------------------------- 
    … …  
    131132         WRITE(numout,*) '              daily mean to diurnal cycle qsr            ln_dm2dc    = ', ln_dm2dc  
    132133         WRITE(numout,*) '              runoff / runoff mouths                     ln_rnf      = ', ln_rnf 
     134         WRITE(numout,*) '              iceshelf formulation                       nn_isf      = ', nn_isf 
    133135         WRITE(numout,*) '              Sea Surface Restoring on SST and/or SSS    ln_ssr      = ', ln_ssr 
    134136         WRITE(numout,*) '              FreshWater Budget control  (=0/1/2)        nn_fwb      = ', nn_fwb 
    … …  
    180182         rnfmsk_z(:)   = 0.0_wp 
    181183      ENDIF 
     184      IF( nn_isf .EQ. 0 ) THEN                      ! no specific treatment in vicinity of ice shelf  
     185         IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' ) 
     186         fwfisf  (:,:) = 0.0_wp 
     187      END IF 
    182188      IF( nn_ice == 0  )   fr_i(:,:) = 0.e0        ! no ice in the domain, ice fraction is always zero 
    183189 
    … …  
    347353 
    348354      IF( ln_icebergs    )   CALL icb_stp( kt )                   ! compute icebergs 
     355 
     356      IF( nn_isf   /= 0  )   CALL sbc_isf( kt )                    ! compute iceshelves 
    349357 
    350358      IF( ln_rnf         )   CALL sbc_rnf( kt )                   ! add runoffs to fresh water fluxes 
    … …  
    422430      ! 
    423431      IF(ln_ctl) THEN         ! print mean trends (used for debugging) 
    424          CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i     - : ', mask1=tmask, ovlap=1 ) 
    425          CALL prt_ctl(tab2d_1=(emp-rnf)        , clinfo1=' emp-rnf  - : ', mask1=tmask, ovlap=1 ) 
    426          CALL prt_ctl(tab2d_1=(sfx-rnf)        , clinfo1=' sfx-rnf  - : ', mask1=tmask, ovlap=1 ) 
     432         CALL prt_ctl(tab2d_1=fr_i              , clinfo1=' fr_i     - : ', mask1=tmask, ovlap=1 ) 
     433         CALL prt_ctl(tab2d_1=(emp-rnf + fwfisf), clinfo1=' emp-rnf  - : ', mask1=tmask, ovlap=1 ) 
     434         CALL prt_ctl(tab2d_1=(sfx-rnf + fwfisf), clinfo1=' sfx-rnf  - : ', mask1=tmask, ovlap=1 ) 
    427435         CALL prt_ctl(tab2d_1=qns              , clinfo1=' qns      - : ', mask1=tmask, ovlap=1 ) 
    428436         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr      - : ', mask1=tmask, ovlap=1 ) 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90

    r4624 r4666  
    1919   USE phycst          ! physical constants 
    2020   USE sbc_oce         ! surface boundary condition variables 
     21   USE sbcisf          ! PM we could remove it I think 
    2122   USE closea          ! closed seas 
    2223   USE fldread         ! read input field at current time step 
    … …  
    2425   USE iom             ! I/O module 
    2526   USE lib_mpp         ! MPP library 
     27   USE eosbn2 
     28   USE wrk_nemo        ! Memory allocation 
    2629 
    2730   IMPLICIT NONE 
    … …  
    98101      INTEGER  ::   z_err = 0 ! dummy integer for error handling 
    99102      !!---------------------------------------------------------------------- 
     103      REAL(wp), DIMENSION(:,:), POINTER       ::   ztfrz   ! freezing point used for temperature correction 
     104      ! 
     105      CALL wrk_alloc( jpi,jpj, ztfrz) 
     106 
    100107      ! 
    101108      IF( kt == nit000 )   CALL sbc_rnf_init                           ! Read namelist and allocate structures 
    … …  
    134141               WHERE( sf_t_rnf(1)%fnow(:,:,1) == -999._wp )             ! if missing data value use SST as runoffs temperature 
    135142                   rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 
     143               END WHERE 
     144               WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp )             ! where fwf comes from melting of ice shelves or iceberg 
     145                   ztfrz(:,:) = -1.9 !tfreez( sss_m(:,:) ) !PM to be discuss (trouble if sensitivity study) 
     146                   rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * lfusisf * r1_rau0_rcp 
    136147               END WHERE 
    137148            ELSE                                                        ! use SST as runoffs temperature 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90

    r4292 r4666  
    5454      INTEGER, INTENT(in) ::   kt   ! ocean time step 
    5555      ! 
     56      INTEGER  ::   ji,jj               ! loop index 
    5657      REAL(wp) ::   zcoef, zf_sbc       ! local scalar 
    5758      !!--------------------------------------------------------------------- 
    … …  
    5960      IF( nn_fsbc == 1 ) THEN                             !   Instantaneous surface fields        ! 
    6061         !                                                ! ---------------------------------------- ! 
    61          ssu_m(:,:) = ub(:,:,1) 
    62          ssv_m(:,:) = vb(:,:,1) 
    63          sst_m(:,:) = tsn(:,:,1,jp_tem) 
    64          sss_m(:,:) = tsn(:,:,1,jp_sal) 
     62         DO jj = 1, jpj 
     63            DO ji = 1, jpi 
     64               ssu_m(ji,jj) = ub(ji,jj,miku(ji,jj)) 
     65               ssv_m(ji,jj) = vb(ji,jj,mikv(ji,jj)) 
     66               sst_m(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) 
     67               sss_m(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) 
     68               IF( lk_vvl )   fse3t_m(ji,jj) = fse3t_n(ji,jj,mikt(ji,jj)) 
     69            END DO 
     70         END DO 
    6571         !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics) 
    6672         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) 
    6773         ELSE                    ;   ssh_m(:,:) = sshn(:,:) 
    6874         ENDIF 
    69          ! 
    70          IF( lk_vvl )   fse3t_m(:,:) = fse3t_n(:,:,1) 
    7175         ! 
    7276      ELSE 
    … …  
    7781            IF(lwp) WRITE(numout,*) '~~~~~~~   mean fields initialised to instantaneous values' 
    7882            zcoef = REAL( nn_fsbc - 1, wp ) 
    79             ssu_m(:,:) = zcoef * ub(:,:,1) 
    80             ssv_m(:,:) = zcoef * vb(:,:,1) 
    81             sst_m(:,:) = zcoef * tsn(:,:,1,jp_tem) 
    82             sss_m(:,:) = zcoef * tsn(:,:,1,jp_sal) 
     83            DO jj = 1, jpj 
     84               DO ji = 1, jpi 
     85                  ssu_m(ji,jj) = zcoef * ub(ji,jj,miku(ji,jj)) 
     86                  ssv_m(ji,jj) = zcoef * vb(ji,jj,mikv(ji,jj)) 
     87                  sst_m(ji,jj) = zcoef * tsn(ji,jj,mikt(ji,jj),jp_tem) 
     88                  sss_m(ji,jj) = zcoef * tsn(ji,jj,mikt(ji,jj),jp_sal) 
     89                  IF( lk_vvl )   fse3t_m(ji,jj) = zcoef * fse3t_n(ji,jj,mikt(ji,jj)) 
     90               END DO 
     91            END DO 
    8392            !                          ! removed inverse barometer ssh when Patm forcing is used  
    8493            IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = zcoef * ( sshn(:,:) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) ) 
    8594            ELSE                    ;   ssh_m(:,:) = zcoef *   sshn(:,:) 
    8695            ENDIF 
    87             IF( lk_vvl )   fse3t_m(:,:) = zcoef * fse3t_n(:,:,1) 
    8896            !                                             ! ---------------------------------------- ! 
    8997         ELSEIF( MOD( kt - 2 , nn_fsbc ) == 0 ) THEN      !   Initialisation: New mean computation   ! 
    … …  
    99107         !                                                !        Cumulate at each time step        ! 
    100108         !                                                ! ---------------------------------------- ! 
    101          ssu_m(:,:) = ssu_m(:,:) + ub(:,:,1) 
    102          ssv_m(:,:) = ssv_m(:,:) + vb(:,:,1) 
    103          sst_m(:,:) = sst_m(:,:) + tsn(:,:,1,jp_tem) 
    104          sss_m(:,:) = sss_m(:,:) + tsn(:,:,1,jp_sal) 
    105          !                          ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics) 
     109         DO jj = 1, jpj 
     110            DO ji = 1, jpi 
     111               ssu_m(ji,jj) = ssu_m(ji,jj) + ub(ji,jj,miku(ji,jj)) 
     112               ssv_m(ji,jj) = ssv_m(ji,jj) + vb(ji,jj,mikv(ji,jj)) 
     113               sst_m(ji,jj) = sst_m(ji,jj) + tsn(ji,jj,mikt(ji,jj),jp_tem) 
     114               sss_m(ji,jj) = sss_m(ji,jj) + tsn(ji,jj,mikt(ji,jj),jp_sal) 
     115               IF( lk_vvl )   fse3t_m(ji,jj) = fse3t_m(ji,jj) + fse3t_n(ji,jj,mikt(ji,jj)) 
     116            END DO 
     117         END DO 
     118!                                   ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics) 
    106119         IF( ln_apr_dyn ) THEN   ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) - 0.5 *  ( ssh_ib(:,:) + ssh_ibb(:,:) ) 
    107120         ELSE                    ;   ssh_m(:,:) = ssh_m(:,:) + sshn(:,:) 
    108121         ENDIF 
    109          IF( lk_vvl )   fse3t_m(:,:) = fse3t_m(:,:) + fse3t_n(:,:,1) 
    110122 
    111123         !                                                ! ---------------------------------------- ! 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90

    r4624 r4666  
    5555   PUBLIC   bn2        ! called by step module 
    5656   PUBLIC   eos_alpbet ! called by ldfslp module 
    57    PUBLIC   tfreez     ! called by sbcice_... modules 
     57   PUBLIC   tfreez     ! called by sbcice_... modules and sbcisf module 
     58   PUBLIC   tfreez1D   ! called by trasbc modules 
    5859 
    5960   !                                  !!* Namelist (nameos) * 
    … …  
    396397      ! 
    397398!CDIR NOVERRCHK 
    398          DO jj = 1, jpjm1 
     399         DO jj = 1, jpj 
    399400!CDIR NOVERRCHK 
    400             DO ji = 1, fs_jpim1   ! vector opt. 
     401            DO ji = 1, jpi   ! vector opt. 
    401402               zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) ) 
    402403            END DO 
    403404         END DO 
    404          DO jj = 1, jpjm1 
    405             DO ji = 1, fs_jpim1   ! vector opt. 
    406                zmask = tmask(ji,jj,1)          ! land/sea bottom mask = surf. mask 
     405         DO jj = 1, jpj 
     406            DO ji = 1, jpi   ! vector opt. 
     407               zmask = lmask(ji,jj)          ! land/sea bottom mask = surf. mask 
    407408               zt    = pts  (ji,jj,jp_tem)            ! interpolated T 
    408409               zs    = pts  (ji,jj,jp_sal)            ! interpolated S 
    … …  
    444445         ! 
    445446      CASE( 1 )                !==  Linear formulation = F( temperature )  ==! 
    446          DO jj = 1, jpjm1 
    447             DO ji = 1, fs_jpim1   ! vector opt. 
    448                prd(ji,jj) = ( 0.0285_wp - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 
     447         DO jj = 1, jpj 
     448            DO ji = 1, jpi   ! vector opt. 
     449               prd(ji,jj) = ( 0.0285_wp - rn_alpha * pts(ji,jj,jp_tem) ) * lmask(ji,jj) 
    449450            END DO 
    450451         END DO 
    451452         ! 
    452453      CASE( 2 )                !==  Linear formulation = F( temperature , salinity )  ==! 
    453          DO jj = 1, jpjm1 
    454             DO ji = 1, fs_jpim1   ! vector opt. 
    455                prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 
     454         DO jj = 1, jpj 
     455            DO ji = 1, jpi   ! vector opt. 
     456               prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * lmask(ji,jj) 
    456457            END DO 
    457458         END DO 
    … …  
    522523               DO ji = 1, jpi 
    523524                  zgde3w = grav / fse3w(ji,jj,jk) 
    524                   zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) )         ! potential temperature at w-pt 
    525                   zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35.0  ! salinity anomaly (s-35) at w-pt 
     525                  zt = 0.5_wp * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) )            ! potential temperature at w-pt 
     526                  zs = 0.5_wp * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35.0_wp  ! salinity anomaly (s-35) at w-pt 
    526527                  zh = fsdepw(ji,jj,jk)                                                ! depth in meters  at w-point 
    527528                  ! 
    … …  
    551552                     &                                - 0.121555e-07_wp ) * zh 
    552553                     ! 
    553                   pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2 
     554                  pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) * tmask(ji,jj,jk-1)           &   ! N^2 
    554555                     &          * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )   & 
    555556                     &                     - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) 
    … …  
    566567      CASE( 1 )                !==  Linear formulation = F( temperature )  ==! 
    567568         DO jk = 2, jpkm1 
    568             pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk) 
     569            pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) )   & 
     570               &          / fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1) 
    569571         END DO 
    570572         ! 
    … …  
    573575            pn2(:,:,jk) = grav * (  rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) )      & 
    574576               &                  - rn_beta  * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) )  )   & 
    575                &               / fse3w(:,:,jk) * tmask(:,:,jk) 
     577               &               / fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1) 
    576578         END DO 
    577579#if defined key_zdfddm 
    … …  
    703705   END FUNCTION tfreez 
    704706 
     707   FUNCTION tfreez1D( psal, pdep ) RESULT( ptf ) 
     708      !!---------------------------------------------------------------------- 
     709      !!                 ***  ROUTINE eos_init  *** 
     710      !! 
     711      !! ** Purpose :   Compute the sea surface freezing temperature [Celcius] 
     712      !! 
     713      !! ** Method  :   UNESCO freezing point at the surface (pressure = 0???) 
     714      !!       freezing point [Celcius]=(-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s-7.53e-4*p 
     715      !!       checkvalue: tf= -2.588567 Celsius for s=40.0psu, p=500. decibars 
     716      !! 
     717      !! Reference  :   UNESCO tech. papers in the marine science no. 28. 1978 
     718      !!---------------------------------------------------------------------- 
     719      REAL(wp), INTENT(in   ) ::   psal   ! salinity             [psu] 
     720      REAL(wp), INTENT(in   ), OPTIONAL ::   pdep   ! pressure             [dBar] 
     721      ! Leave result array automatic rather than making explicitly allocated 
     722      REAL(wp)                ::   ptf    ! freezing temperature [Celcius] 
     723      !!---------------------------------------------------------------------- 
     724      ! 
     725      ptf = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal )   & 
     726         &                     - 2.154996e-4_wp *  psal ) * psal 
     727      IF ( PRESENT( pdep ) ) THEN    
     728         ptf = ptf - 7.53e-4_wp * pdep 
     729      ENDIF 
     730      ! 
     731   END FUNCTION tfreez1D 
     732 
     733 
    705734 
    706735   SUBROUTINE eos_init 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90

    r4499 r4666  
    3333   USE wrk_nemo        ! Memory Allocation 
    3434   USE timing          ! Timing 
     35   USE phycst 
    3536 
    3637   IMPLICIT NONE 
    … …  
    121122      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend  
    122123      ! 
    123       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    124       INTEGER  ::   ierr             ! local integer 
     124      INTEGER  ::   ji, jj, jk, jn, ik   ! dummy loop indices 
     125      INTEGER  ::   ierr                 ! local integer 
    125126      REAL(wp) ::   zbtr, ztra                            ! local scalars 
    126127      REAL(wp) ::   zfp_ui, zfp_vj, zfp_w, zcofi          !   -      - 
    … …  
    128129      REAL(wp) ::   zupsut, zcenut, zupst                 !   -      - 
    129130      REAL(wp) ::   zupsvt, zcenvt, zcent, zice           !   -      - 
    130       REAL(wp), POINTER, DIMENSION(:,:  ) :: ztfreez  
     131      REAL(wp), POINTER, DIMENSION(:,:  ) :: ztfreez, zpress  
    131132      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zind 
    132133      !!---------------------------------------------------------------------- 
    … …  
    134135      IF( nn_timing == 1 )  CALL timing_start('tra_adv_cen2') 
    135136      ! 
    136       CALL wrk_alloc( jpi, jpj, ztfreez ) 
     137      CALL wrk_alloc( jpi, jpj, ztfreez, zpress ) 
    137138      CALL wrk_alloc( jpi, jpj, jpk, zwz, zind ) 
    138139      ! 
    … …  
    169170!!gm  not strickly exact : the freezing point should be computed at each ocean levels... 
    170171!!gm  not a big deal since cen2 is no more used in global ice-ocean simulations 
    171       ztfreez(:,:) = tfreez( tsn(:,:,1,jp_sal) ) 
     172!!ch  changes for ice shelf to retain standard behaviour elsewhere, even if not optimal  
     173      DO jj = 1, jpj  
     174         DO ji = 1, jpi  
     175            ik=mikt(ji,jj)  
     176            IF (ik > 1 ) THEN  
     177               zpress(ji,jj) = grav*rau0*fsdept(ji,jj,ik)*1.e-04   
     178            ELSE  
     179               zpress(ji,jj) = 0.0  
     180            ENDIF   
     181         END DO  
     182      END DO  
     183      ztfreez(:,:) = tfreez( tsn(:,:,1, jp_sal), zpress(:,:) ) 
     184       
    172185      DO jk = 1, jpk 
    173186         DO jj = 1, jpj 
    … …  
    224237         !                                                     ! Surface value :  
    225238         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                         ! volume variable 
    226          ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)   ! linear free surface  
     239         ELSE 
     240            DO jj = 1, jpj   ! vector opt. 
     241               DO ji = 1, jpi   ! vector opt. 
     242                  ik=mikt(ji,jj)                 
     243                  zwz(ji,jj,ik ) = pwn(ji,jj,ik) * ptn(ji,jj,ik,jn)   ! linear free surface  
     244                  zwz(ji,jj,1:ik-1) = 0.e0 
     245               END DO 
     246            END DO 
    227247         ENDIF 
    228248         ! 
    … …  
    281301      ENDIF 
    282302      ! 
    283       CALL wrk_dealloc( jpi, jpj, ztfreez ) 
     303      CALL wrk_dealloc( jpi, jpj, ztfreez, zpress ) 
    284304      CALL wrk_dealloc( jpi, jpj, jpk, zwz, zind ) 
    285305      ! 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90

    r4499 r4666  
    7777      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend  
    7878      ! 
    79       INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices   
     79      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices 
     80      INTEGER  ::   ik   
    8081      REAL(wp) ::   z2dtt, zbtr, ztra        ! local scalar 
    8182      REAL(wp) ::   zfp_ui, zfp_vj, zfp_wk   !   -      - 
    … …  
    103104      ENDIF 
    104105      ! 
    105       zwi(:,:,:) = 0.e0 
     106      zwi(:,:,:) = 0.e0 ; zwz(:,:,:) = 0.e0 
    106107      ! 
    107108      !                                                          ! =========== 
    … …  
    132133         ! upstream tracer flux in the k direction 
    133134         ! Surface value 
    134          IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                         ! volume variable 
    135          ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn)   ! linear free surface  
     135         IF( lk_vvl ) THEN    
     136            DO jj = 1, jpj 
     137               DO ji = 1, jpi 
     138                  zwz(ji,jj, mikt(ji,jj) ) = 0.e0                         ! volume variable 
     139               END DO 
     140            END DO 
     141         ELSE                 
     142            DO jj = 1, jpj 
     143               DO ji = 1, jpi 
     144                  zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)   ! linear free surface  
     145               END DO 
     146            END DO    
    136147         ENDIF 
    137148         ! Interior value 
    138          DO jk = 2, jpkm1 
    139             DO jj = 1, jpj 
    140                DO ji = 1, jpi 
     149         DO jj = 1, jpj 
     150            DO ji = 1, jpi 
     151               DO jk = mikt(ji,jj)+1, jpkm1 
    141152                  zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 
    142153                  zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 
    … …  
    157168                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) 
    158169                  ! update and guess with monotonic sheme 
    159                   pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra 
     170                  pta(ji,jj,jk,jn) =   pta(ji,jj,jk,jn)         + ztra   * tmask(ji,jj,jk) 
    160171                  zwi(ji,jj,jk)    = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk) 
    161172               END DO 
    … …  
    191202         zwz(:,:,1) = 0.e0         ! Surface value 
    192203         ! 
    193          DO jk = 2, jpkm1          ! Interior value 
    194             DO jj = 1, jpj 
    195                DO ji = 1, jpi 
     204         DO jj = 1, jpj 
     205            DO ji = 1, jpi 
     206               ik=mikt(ji,jj) 
     207               ! surface value 
     208               zwz(ji,jj,1:ik) = 0.e0 
     209               ! Interior value 
     210               DO jk = mikt(ji,jj)+1, jpkm1                     
    196211                  zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk) 
    197212               END DO 
    … …  
    217232                     &             + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) 
    218233                  ! add them to the general tracer trends 
    219                   pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
     234                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra * tmask(ji,jj,jk) 
    220235               END DO 
    221236            END DO 
    … …  
    281296      zbig  = 1.e+40_wp 
    282297      zrtrn = 1.e-15_wp 
    283       zbetup(:,:,jpk) = 0._wp   ;   zbetdo(:,:,jpk) = 0._wp 
     298      zbetup(:,:,:) = 0._wp   ;   zbetdo(:,:,:) = 0._wp 
    284299 
    285300 
    … …  
    292307         &        paft * tmask + zbig * ( 1.e0 - tmask )  ) 
    293308 
    294       DO jk = 1, jpkm1 
    295          ikm1 = MAX(jk-1,1) 
    296          z2dtt = p2dt(jk) 
    297          DO jj = 2, jpjm1 
    298             DO ji = fs_2, fs_jpim1   ! vector opt. 
    299  
     309      DO jj = 2, jpjm1 
     310         DO ji = fs_2, fs_jpim1   ! vector opt. 
     311            DO jk = mikt(ji,jj), jpkm1 
     312               ikm1 = MAX(jk-1,mikt(ji,jj)) 
     313               z2dtt = p2dt(jk) 
     314                
    300315               ! search maximum in neighbourhood 
    301316               zup = MAX(  zbup(ji  ,jj  ,jk  ),   & 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90

    r4624 r4666  
    420420#endif 
    421421            ik = mbkt(ji,jj)                        ! bottom T-level index 
    422             ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) * tmask(ji,jj,1)      ! bottom before T and S 
    423             zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) * tmask(ji,jj,1) 
     422            ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) * lmask(ji,jj)      ! bottom before T and S 
     423            zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) * lmask(ji,jj) 
    424424            zdep(ji,jj) = gdept_0(ji,jj,ik)         ! bottom T-level reference depth 
    425425            ! 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90

    r4488 r4666  
    7575 
    7676      SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend 
    77       CASE ( 0 )   ;   CALL tra_ldf_lap     ( kt, nit000, 'TRA', gtsu, gtsv, tsb, tsa, jpts        )  ! iso-level laplacian 
     77      CASE ( 0 )   ;   CALL tra_ldf_lap     ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        & 
     78                               &                                   tsb, tsa, jpts        )  ! iso-level laplacian 
    7879      CASE ( 1 )                                                                              ! rotated laplacian 
    7980         IF( ln_traldf_grif ) THEN                                                           
    8081                       CALL tra_ldf_iso_grif( kt, nit000,'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )      ! Griffies operator 
    8182         ELSE                                                                                 
    82                        CALL tra_ldf_iso     ( kt, nit000, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )      ! Madec operator 
    83          ENDIF 
    84       CASE ( 2 )   ;   CALL tra_ldf_bilap   ( kt, nit000, 'TRA', gtsu, gtsv, tsb, tsa, jpts        )  ! iso-level bilaplacian 
     83                       CALL tra_ldf_iso     ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        & 
     84                               &                                  tsb, tsa, jpts, ahtb0 )      ! Madec operator 
     85         ENDIF 
     86      CASE ( 2 )   ;   CALL tra_ldf_bilap   ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        & 
     87                               &                                   tsb, tsa, jpts        )  ! iso-level bilaplacian 
    8588      CASE ( 3 )   ;   CALL tra_ldf_bilapg  ( kt, nit000, 'TRA',             tsb, tsa, jpts        )  ! s-coord. geopot. bilap. 
    8689         ! 
    8790      CASE ( -1 )                                ! esopa: test all possibility with control print 
    88          CALL tra_ldf_lap   ( kt, nit000, 'TRA', gtsu, gtsv, tsb, tsa, jpts        )  
     91         CALL tra_ldf_lap   ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        & 
     92         &                                       tsb, tsa, jpts        )  
    8993         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf0 - Ta: ', mask1=tmask,               & 
    9094         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    … …  
    9296            CALL tra_ldf_iso_grif( kt, nit000, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 ) 
    9397         ELSE 
    94             CALL tra_ldf_iso     ( kt, nit000, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )   
     98            CALL tra_ldf_iso     ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        & 
     99            &                                               tsb, tsa, jpts, ahtb0 )   
    95100         ENDIF 
    96101         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf1 - Ta: ', mask1=tmask,               & 
    97102         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    98          CALL tra_ldf_bilap ( kt, nit000, 'TRA', gtsu, gtsv, tsb, tsa, jpts        )  
     103         CALL tra_ldf_bilap ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        & 
     104         &                                       tsb, tsa, jpts        )  
    99105         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf2 - Ta: ', mask1=tmask,               & 
    100106         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' ) 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90

    r4292 r4666  
    4949CONTAINS 
    5050  
    51    SUBROUTINE tra_ldf_bilap( kt, kit000, cdtype, pgu, pgv,      & 
     51   SUBROUTINE tra_ldf_bilap( kt, kit000, cdtype, pgu, pgv,            & 
     52      &                                          pgui, pgvi,          & 
    5253      &                                  ptb, pta, kjpt )   
    5354      !!---------------------------------------------------------------------- 
    … …  
    8283      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator) 
    8384      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
    84       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
     85      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgu , pgv  ! tracer gradient at pstep levels 
     86      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at pstep levels 
    8587      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields 
    8688      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
    … …  
    128130                     IF( mbku(ji,jj) == jk )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgu(ji,jj,jn) 
    129131                     IF( mbkv(ji,jj) == jk )  ztv(ji,jj,jk) = zeev(ji,jj) * pgv(ji,jj,jn) 
     132                     ! (ISH) 
     133                     IF( miku(ji,jj) == jk )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgui(ji,jj,jn) 
     134                     IF( mikv(ji,jj) == jk )  ztu(ji,jj,jk) = zeev(ji,jj) * pgvi(ji,jj,jn) 
    130135                  END DO 
    131136               END DO 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90

    r4292 r4666  
    5252 
    5353   SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pgu, pgv,              & 
     54      &                                pgui, pgvi,                    & 
    5455      &                                ptb, pta, kjpt, pahtb0 ) 
    5556      !!---------------------------------------------------------------------- 
    … …  
    110111      REAL(wp)                         ::   zztmp               ! local scalar 
    111112#endif 
    112       REAL(wp), POINTER, DIMENSION(:,:  ) ::  zdkt, zdk1t, z2d 
    113       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdit, zdjt, ztfw  
     113      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d 
     114      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdkt, zdk1t, zdit, zdjt, ztfw  
    114115      !!---------------------------------------------------------------------- 
    115116      ! 
    116117      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_iso') 
    117118      ! 
    118       CALL wrk_alloc( jpi, jpj,      zdkt, zdk1t, z2d )  
    119       CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt, ztfw  )  
     119      CALL wrk_alloc( jpi, jpj,      z2d )  
     120      CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t )  
    120121      ! 
    121122 
    … …  
    150151            DO jj = 1, jpjm1 
    151152               DO ji = 1, fs_jpim1   ! vector opt. 
    152                   zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
    153                   zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)       
     153! IF useless if zpshde defines pgu everywhere 
     154                  IF (mbku(ji,jj) > 1) zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)           
     155                  IF (mbkv(ji,jj) > 1) zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 
     156                  ! (ISF) 
     157                  IF (miku(ji,jj) > 1) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)           
     158                  IF (mikv(ji,jj) > 1) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)      
    154159               END DO 
    155160            END DO 
    … …  
    161166!CDIR PARALLEL DO PRIVATE( zdk1t )  
    162167         !                                                ! =============== 
    163          DO jk = 1, jpkm1                                 ! Horizontal slab 
     168         DO jj = 1, jpj                                 ! Horizontal slab 
    164169            !                                             ! =============== 
    165             ! 1. Vertical tracer gradient at level jk and jk+1 
    166             ! ------------------------------------------------ 
    167             ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 
    168             zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1) 
    169             ! 
    170             IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:) 
    171             ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk) 
    172             ENDIF 
     170            DO ji = 1, jpi   ! vector opt. 
     171               DO jk = mikt(ji,jj), jpkm1 
     172               ! 1. Vertical tracer gradient at level jk and jk+1 
     173               ! ------------------------------------------------ 
     174               ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2) 
     175                  zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn) - ptb(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1) 
     176               ! 
     177                  IF( jk == mikt(ji,jj) ) THEN  ;   zdkt(ji,jj,jk) = zdk1t(ji,jj,jk) 
     178                  ELSE                          ;   zdkt(ji,jj,jk) = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
     179                  ENDIF 
     180               END DO 
     181            END DO 
     182         END DO 
    173183 
    174184            ! 2. Horizontal fluxes 
    175185            ! --------------------    
    176             DO jj = 1 , jpjm1 
    177                DO ji = 1, fs_jpim1   ! vector opt. 
     186         DO jj = 1 , jpjm1 
     187            DO ji = 1, fs_jpim1   ! vector opt. 
     188               DO jk = mikt(ji,jj), jpkm1 
    178189                  zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,jk) 
    179190                  zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,jk) 
    … …  
    189200                  ! 
    190201                  zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   & 
    191                      &              + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      & 
    192                      &                         + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk) 
     202                     &              + zcof1 * (  zdkt (ji+1,jj,jk) + zdk1t(ji,jj,jk)      & 
     203                     &                         + zdk1t(ji+1,jj,jk) + zdkt (ji,jj,jk)  )  ) * umask(ji,jj,jk) 
    193204                  zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   & 
    194                      &              + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      & 
    195                      &                         + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)                   
    196                END DO 
    197             END DO 
     205                     &              + zcof2 * (  zdkt (ji,jj+1,jk) + zdk1t(ji,jj,jk)      & 
     206                     &                         + zdk1t(ji,jj+1,jk) + zdkt (ji,jj,jk)  )  ) * vmask(ji,jj,jk)                   
     207               END DO 
     208            END DO 
     209         END DO 
    198210 
    199211            ! II.4 Second derivative (divergence) and add to the general trend 
    200212            ! ---------------------------------------------------------------- 
    201             DO jj = 2 , jpjm1 
    202                DO ji = fs_2, fs_jpim1   ! vector opt. 
    203                   zbtr = 1.0 / ( e12t(ji,jj) * fse3t_n(ji,jj,jk) ) 
     213         DO jj = 2 , jpjm1 
     214            DO ji = fs_2, fs_jpim1   ! vector opt. 
     215               DO jk = mikt(ji,jj), jpkm1 
     216                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 
    204217                  ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  ) 
    205218                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 
    … …  
    264277            DO jj = 2, jpjm1 
    265278               DO ji = fs_2, fs_jpim1   ! vector opt. 
    266                   zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk) 
     279                  zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
    267280                  ! 
    268281                  zmsku = 1./MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)      & 
    … …  
    297310      END DO 
    298311      ! 
    299       CALL wrk_dealloc( jpi, jpj,      zdkt, zdk1t, z2d )  
    300       CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt, ztfw  )  
     312      CALL wrk_dealloc( jpi, jpj, z2d )  
     313      CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t )  
    301314      ! 
    302315      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_iso') 
    … …  
    309322   !!---------------------------------------------------------------------- 
    310323CONTAINS 
    311    SUBROUTINE tra_ldf_iso( kt, kit000,cdtype, pgu, pgv, ptb, pta, kjpt, pahtb0 )      ! Empty routine 
     324   SUBROUTINE tra_ldf_iso( kt, kit000,cdtype, pgu, pgv, pgui, pgvi, ptb, pta, kjpt, pahtb0 )      ! Empty routine 
    312325      INTEGER:: kt, kit000 
    313326      CHARACTER(len=3) ::   cdtype 
    314       REAL, DIMENSION(:,:,:) ::   pgu, pgv   ! tracer gradient at pstep levels 
     327      REAL, DIMENSION(:,:,:) ::   pgu, pgv, pgui, pgvi    ! tracer gradient at pstep levels 
    315328      REAL, DIMENSION(:,:,:,:) ::   ptb, pta 
    316329      WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', kt, kit000, cdtype,   & 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90

    r4364 r4666  
    4343CONTAINS 
    4444 
    45    SUBROUTINE tra_ldf_lap( kt, kit000, cdtype, pgu, pgv,      & 
     45   SUBROUTINE tra_ldf_lap( kt, kit000, cdtype, pgu , pgv ,    & 
     46      &                                        pgui, pgvi,    & 
    4647      &                                ptb, pta, kjpt )  
    4748      !!---------------------------------------------------------------------- 
    … …  
    6970      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers 
    7071      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels 
     72      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top levels 
    7173      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! before and now tracer fields 
    7274      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend  
    … …  
    114116                        ztv(ji,jj,jk) = zabe2 * pgv(ji,jj,jn) 
    115117                     ENDIF 
     118                      
     119                     ! (ISH) 
     120                     ! ice shelf level level  
     121                     iku = miku(ji,jj)  
     122                     ikv = mikv(ji,jj)  
     123                     IF( iku == jk ) THEN  
     124                        zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * re2u_e1u(ji,jj) * fse3u_n(ji,jj,iku)  
     125                        ztu(ji,jj,jk) = zabe1 * pgui(ji,jj,jn)  
     126                     ENDIF  
     127                     IF( ikv == jk ) THEN  
     128                        zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * re1v_e2v(ji,jj) * fse3v_n(ji,jj,ikv)  
     129                        ztv(ji,jj,jk) = zabe2 * pgvi(ji,jj,jn)  
     130                     END IF  
    116131                  END DO 
    117132               END DO 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r4624 r4666  
    294294                  DO jj = 2, jpjm1 
    295295                     DO ji = fs_2, fs_jpim1   ! vector opt. 
    296                         qsr_hc(ji,jj,jk) =  etot3(ji,jj,jk) * qsr(ji,jj) 
     296                        ! (ISF) no light penetration below the ice shelves          
     297                        qsr_hc(ji,jj,jk) =  etot3(ji,jj,jk) * qsr(ji,jj) * tmask(ji,jj,1) 
    297298                     END DO 
    298299                  END DO 
    … …  
    520521                  ! 
    521522                  DO jk = 1, nksr 
    522                      etot3(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) )  
     523                     ! (ISF) no light penetration below the ice shelves 
     524                     etot3(:,:,jk) = r1_rau0_rcp * ( zea(:,:,jk) - zea(:,:,jk+1) ) * tmask(:,:,1) 
    523525                  END DO 
    524526                  etot3(:,:,nksr+1:jpk) = 0.e0                ! below 400m set to zero 
    … …  
    548550                        zc0 = zz0 * EXP( -fsdepw(ji,jj,jk  )*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk  )*xsi1r ) 
    549551                        zc1 = zz0 * EXP( -fsdepw(ji,jj,jk+1)*xsi0r ) + zz1 * EXP( -fsdepw(ji,jj,jk+1)*xsi1r ) 
    550                         etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  )  
     552                        etot3(ji,jj,jk) = (  zc0 * tmask(ji,jj,jk) - zc1 * tmask(ji,jj,jk+1)  ) * tmask(ji,jj,1)  
    551553                     END DO 
    552554                  END DO 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/trasbc.F90

    r3764 r4666  
    2424   USE prtctl          ! Print control 
    2525   USE sbcrnf          ! River runoff   
     26   USE sbcisf          ! Ice shelf    
    2627   USE sbcmod          ! ln_rnf   
    2728   USE iom 
    … …  
    2930   USE wrk_nemo        ! Memory Allocation 
    3031   USE timing          ! Timing 
     32   USE eosbn2 
    3133 
    3234   IMPLICIT NONE 
    … …  
    112114      !! 
    113115      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices   
     116      INTEGER  ::   ikt, ikb  
     117      INTEGER  ::   nk_isf 
    114118      REAL(wp) ::   zfact, z1_e3t, zdep 
     119      REAL(wp) ::   zalpha, zhk 
     120      REAL(wp) ::  zt_frz, zpress 
    115121      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds 
    116122      !!---------------------------------------------------------------------- 
    … …  
    205211      ENDIF 
    206212      ! 
     213      ! 
     214      !---------------------------------------- 
     215      !       Ice Shelf effects (ISF) 
     216      !     tbl treated as in Losh (2008) JGR 
     217      !---------------------------------------- 
     218      ! 
     219      IF (nn_isf .GT. 0) THEN 
     220         zfact = 0.5e0 
     221         DO jj = 2, jpj 
     222            DO ji = fs_2, fs_jpim1 
     223          
     224               ikt = misfkt(ji,jj) 
     225               ikb = misfkb(ji,jj) 
     226    
     227               ! level fully include in the ice shelf boundary layer 
     228               ! if isfdiv, we have to remove heat flux due to inflow at 0oC (as in rnf when you add rnf at sst) 
     229               ! sign - because fwf sign of evapo (rnf sign of precip) 
     230               DO jk = ikt, ikb - 1 
     231               ! compute tfreez for the temperature correction (we add water at freezing temperature) 
     232                  zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04 
     233                  zt_frz = -1.9 !tfreez1D( tsn(ji,jj,jk,jp_sal), zpress ) 
     234               ! compute trend 
     235                  tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem)                                        & 
     236                     &           + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)        & 
     237                     &               - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) & 
     238                     &           * r1_hisf_tbl(ji,jj) 
     239                  tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal)                                        & 
     240                     &           + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) 
     241               END DO 
     242    
     243               ! level partially include in ice shelf boundary layer  
     244               zhk   = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1 
     245               zalpha = rhisf_tbl(ji,jj) * ( 1 - zhk ) / fse3t(ji,jj,ikb)     ! proportion of bottom cell influenced by boundary layer 
     246               ! compute tfreez for the temperature correction (we add water at freezing temperature) 
     247               zpress = grav*rau0*fsdept(ji,jj,ikb)*1.e-04 
     248               zt_frz = -1.9 !tfreez1D( tsn(ji,jj,ikb,jp_sal), zpress ) 
     249               ! compute trend 
     250               tsa(ji,jj,ikb,jp_tem) = tsa(ji,jj,ikb,jp_tem)                                         & 
     251                  &              + zfact * (risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem)        & 
     252                  &                  - rdivisf * (fwfisf(ji,jj) + fwfisf_b(ji,jj)) * zt_frz * r1_rau0) &  
     253                  &              * r1_hisf_tbl(ji,jj) * zalpha 
     254               tsa(ji,jj,ikb,jp_sal) = tsa(ji,jj,ikb,jp_sal)                                         & 
     255                  &              + zfact * (risf_tsc_b(ji,jj,jp_sal) + risf_tsc(ji,jj,jp_sal)) * r1_hisf_tbl(ji,jj) * zalpha  
     256            END DO 
     257         END DO 
     258         IF( lrst_oce ) THEN 
     259            IF(lwp) WRITE(numout,*) 
     260            IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ',   & 
     261               &                    'at it= ', kt,' date= ', ndastp 
     262            IF(lwp) WRITE(numout,*) '~~~~' 
     263            CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:) ) 
     264            CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b', risf_tsc(:,:,jp_tem) ) 
     265            CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b', risf_tsc(:,:,jp_sal) ) 
     266         ENDIF 
     267      END IF 
     268      ! 
    207269      !---------------------------------------- 
    208270      !        River Runoff effects 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90

    r3294 r4666  
    120120            ELSE                                            ;   zwt(:,:,2:jpk) = fsavs(:,:,2:jpk) 
    121121            ENDIF 
    122             zwt(:,:,1) = 0._wp 
    123             ! 
     122            DO jj=1, jpj 
     123               DO ji=1, jpi 
     124                  zwt(ji,jj,1:mikt(ji,jj)) = 0._wp 
     125               END DO 
     126            END DO 
     127! 
    124128#if defined key_ldfslp 
    125129            ! isoneutral diffusion: add the contribution  
    … …  
    180184            DO jj = 2, jpjm1 
    181185               DO ji = fs_2, fs_jpim1 
    182                   zwt(ji,jj,1) = zwd(ji,jj,1) 
    183                END DO 
    184             END DO 
    185             DO jk = 2, jpkm1 
    186                DO jj = 2, jpjm1 
    187                   DO ji = fs_2, fs_jpim1 
    188                     zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 
     186                  zwt(ji,jj,1:mikt(ji,jj)) = zwd(ji,jj,1:mikt(ji,jj)) 
     187                  DO jk = mikt(ji,jj)+1, jpkm1 
     188                     zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 
    189189                  END DO 
    190190               END DO 
    … …  
    196196         DO jj = 2, jpjm1 
    197197            DO ji = fs_2, fs_jpim1 
    198                ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1) 
    199                ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1) 
    200                pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn) 
    201             END DO 
    202          END DO 
    203          DO jk = 2, jpkm1 
    204             DO jj = 2, jpjm1 
    205                DO ji = fs_2, fs_jpim1 
     198               ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,mikt(ji,jj)) 
     199               ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,mikt(ji,jj)) 
     200               pta(ji,jj,mikt(ji,jj),jn) = ze3tb * ptb(ji,jj,mikt(ji,jj),jn)                     & 
     201                  &                      + p2dt(mikt(ji,jj)) * ze3tn * pta(ji,jj,mikt(ji,jj),jn) 
     202               DO jk = mikt(ji,jj)+1, jpkm1 
    206203                  ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,jk) 
    207204                  ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t  (ji,jj,jk) 
    … …  
    216213            DO ji = fs_2, fs_jpim1 
    217214               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 
    218             END DO 
    219          END DO 
    220          DO jk = jpk-2, 1, -1 
    221             DO jj = 2, jpjm1 
    222                DO ji = fs_2, fs_jpim1 
     215               DO jk = jpk-2, mikt(ji,jj), -1 
    223216                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   & 
    224217                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk) 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r3294 r4666  
    4040 
    4141   SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv,   & 
    42                                  prd, pgru, pgrv    ) 
     42      &                          prd, pgru, pgrv, pmru, pmrv, pgzu, pgzv, pge3ru, pge3rv,  & 
     43      &                   sgtu, sgtv, sgru, sgrv, smru, smrv, sgzu, sgzv, sge3ru, sge3rv ) 
    4344      !!---------------------------------------------------------------------- 
    4445      !!                     ***  ROUTINE zps_hde  *** 
    … …  
    8889      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields 
    8990      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts  
     91      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  sgtu, sgtv  ! hor. grad. of stra at u- & v-pts (ISF) 
    9092      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
    91       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad. of prd at u- & v-pts  
     93      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv      ! hor. grad of prd at u- & v-pts (bottom) 
     94      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmru, pmrv      ! hor. sum  of prd at u- & v-pts (bottom) 
     95      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu, pgzv      ! hor. grad of z   at u- & v-pts (bottom) 
     96      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru, pge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (bottom) 
     97      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  sgru, sgrv      ! hor. grad of prd at u- & v-pts (top) 
     98      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  smru, smrv      ! hor. sum  of prd at u- & v-pts (top) 
     99      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  sgzu, sgzv      ! hor. grad of z   at u- & v-pts (top) 
     100      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  sge3ru, sge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (top) 
    92101      ! 
    93102      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
    94       INTEGER  ::   iku, ikv, ikum1, ikvm1   ! partial step level (ocean bottom level) at u- and v-points 
    95       REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars 
     103      INTEGER  ::   iku, ikv, ikum1, ikvm1,ikup1, ikvp1   ! partial step level (ocean bottom level) at u- and v-points 
     104      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv, zdzwu, zdzwv, zdzwuip1, zdzwvjp1  ! temporary scalars 
    96105      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zri, zrj, zhi, zhj 
    97106      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zti, ztj    ! interpolated value of tracer 
    … …  
    103112      CALL wrk_alloc( jpi, jpj, kjpt, zti, ztj           )  
    104113      ! 
     114      pgru(:,:)=0.0_wp ; pgrv(:,:)=0.0_wp ; pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ; 
     115      ! 
    105116      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
    106117         ! 
    … …  
    114125               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    115126               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
    116                ze3wu = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
    117                ze3wv = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
     127               ! (ISF) case partial step top and bottom in adjacent cell in vertical 
     128               ze3wu  = (fsdept(ji+1,jj,iku) - fsdepw(ji+1,jj,iku)) - (fsdept(ji,jj,iku) - fsdepw(ji,jj,iku)) 
     129               ze3wv  = (fsdept(ji,jj+1,ikv) - fsdepw(ji,jj+1,ikv)) - (fsdept(ji,jj,ikv) - fsdepw(ji,jj,ikv)) 
    118130               ! 
    119131               ! i- direction 
     132               IF (iku .GT. 1) THEN 
    120133               IF( ze3wu >= 0._wp ) THEN      ! case 1 
    121134                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku) 
    … …  
    123136                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
    124137                  ! gradient of  tracers 
    125                   pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     138                  pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     139                  pgzu(ji,jj)    = (fsde3w(ji+1,jj,iku) - ze3wu) - fsde3w(ji,jj,iku) 
    126140               ELSE                           ! case 2 
    127141                  zmaxu = -ze3wu / fse3w(ji,jj,iku) 
    … …  
    129143                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
    130144                  ! gradient of tracers 
    131                   pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     145                  pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     146                  pgzu(ji,jj)    = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) + ze3wu) 
     147               ENDIF 
    132148               ENDIF 
    133149               ! 
    134150               ! j- direction 
     151               IF (ikv .GT. 1) THEN 
    135152               IF( ze3wv >= 0._wp ) THEN      ! case 1 
    136153                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv) 
    … …  
    138155                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
    139156                  ! gradient of tracers 
    140                   pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     157                  pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     158                  pgzv(ji,jj)    = (fsde3w(ji,jj+1,ikv) - ze3wv) - fsde3w(ji,jj,ikv)  
    141159               ELSE                           ! case 2 
    142160                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv) 
    … …  
    144162                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
    145163                  ! gradient of tracers 
    146                   pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     164                  pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     165                  pgzv(ji,jj)    = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) + ze3wv) 
     166               ENDIF 
    147167               ENDIF 
    148168# if ! defined key_vectopt_loop 
    … …  
    165185               iku = mbku(ji,jj) 
    166186               ikv = mbkv(ji,jj) 
    167                ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
    168                ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
     187               ze3wu  = (fsdept(ji+1,jj,iku) - fsdepw(ji+1,jj,iku)) - (fsdept(ji,jj,iku) - fsdepw(ji,jj,iku)) 
     188               ze3wv  = (fsdept(ji,jj+1,ikv) - fsdepw(ji,jj+1,ikv)) - (fsdept(ji,jj,ikv) - fsdepw(ji,jj,ikv)) 
     189 
    169190               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1 
    170191               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2 
    … …  
    177198# endif 
    178199         END DO 
    179  
     200          
    180201         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
    181202         ! step and store it in  zri, zrj for each  case 
    … …  
    191212            DO ji = 1, jpim1 
    192213# endif 
    193                iku = mbku(ji,jj) 
    194                ikv = mbkv(ji,jj) 
    195                ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
    196                ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
    197                IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1 
    198                ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2 
    199                ENDIF 
    200                IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )   ! j: 1 
    201                ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2 
    202                ENDIF 
    203 # if ! defined key_vectopt_loop 
    204             END DO 
    205 # endif 
    206          END DO 
    207          CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions 
     214               iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
     215               ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 )    ! last and before last ocean level at u- & v-points 
     216               ze3wu  = (fsdept(ji+1,jj,iku) - fsdepw(ji+1,jj,iku)) - (fsdept(ji,jj,iku) - fsdepw(ji,jj,iku)) 
     217               ze3wv  = (fsdept(ji,jj+1,ikv) - fsdepw(ji,jj+1,ikv)) - (fsdept(ji,jj,ikv) - fsdepw(ji,jj,ikv)) 
     218               IF( ze3wu >= 0._wp ) THEN  
     219                  pgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1 
     220                  pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) + prd(ji,jj,iku) )   ! i:  
     221                  pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  & 
     222                                * ( (fse3w(ji+1,jj,iku) - ze3wu )* ( zri(ji  ,jj    ) + prd(ji+1,jj,ikum1) + 2._wp) & 
     223                                   - fse3w(ji  ,jj,iku)          * ( prd(ji  ,jj,iku) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2 
     224               ELSE   
     225                  pgru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2 
     226                  pmru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) )   ! i: 2 
     227                  pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  & 
     228                                * (  fse3w(ji+1,jj,iku)          * ( prd(ji+1,jj,iku) + prd(ji+1,jj,ikum1) + 2._wp) & 
     229                                   -(fse3w(ji  ,jj,iku) + ze3wu) * ( zri(ji  ,jj    ) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2 
     230               ENDIF 
     231               IF( ze3wv >= 0._wp ) THEN 
     232                  pgrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )   ! j: 1 
     233                  pmrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )   ! j: 1 
     234                  pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  & 
     235                                * ( (fse3w(ji,jj+1,ikv) - ze3wv )* ( zrj(ji,jj      ) + prd(ji,jj+1,ikvm1) + 2._wp) & 
     236                                   - fse3w(ji,jj  ,ikv)          * ( prd(ji,jj  ,ikv) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2 
     237               ELSE  
     238                  pgrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2 
     239                  pmrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )   ! j: 2 
     240                  pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  & 
     241                                * (  fse3w(ji,jj+1,ikv)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikvm1) + 2._wp) & 
     242                                   -(fse3w(ji,jj  ,ikv) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2 
     243               ENDIF 
     244# if ! defined key_vectopt_loop 
     245            END DO 
     246# endif 
     247         END DO 
     248         CALL lbc_lnk( pgru   , 'U', -1. )   ;   CALL lbc_lnk( pgrv   , 'V', -1. )   ! Lateral boundary conditions 
     249         CALL lbc_lnk( pmru   , 'U',  1. )   ;   CALL lbc_lnk( pmrv   , 'V',  1. )   ! Lateral boundary conditions 
     250         CALL lbc_lnk( pgzu   , 'U', -1. )   ;   CALL lbc_lnk( pgzv   , 'V', -1. )   ! Lateral boundary conditions 
     251         CALL lbc_lnk( pge3ru , 'U', -1. )   ;   CALL lbc_lnk( pge3rv , 'V', -1. )   ! Lateral boundary conditions 
    208252         ! 
    209253      END IF 
    210       ! 
    211       CALL wrk_dealloc( jpi, jpj,       zri, zrj, zhi, zhj )  
     254         ! (ISH)  compute grui and gruvi 
     255      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            ! 
     256# if defined key_vectopt_loop 
     257         jj = 1 
     258         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled) 
     259# else 
     260         DO jj = 1, jpjm1 
     261            DO ji = 1, jpim1 
     262# endif 
     263               iku = miku(ji,jj)   ;  ikup1 = miku(ji,jj) + 1 
     264               ikv = mikv(ji,jj)   ;  ikvp1 = mikv(ji,jj) + 1 
     265               ze3wu  = fse3w(ji+1,jj  ,iku+1) - fse3w(ji,jj,iku+1) 
     266               ze3wv  = fse3w(ji  ,jj+1,ikv+1) - fse3w(ji,jj,ikv+1) 
     267               ! 
     268               ! i- direction 
     269               IF( ze3wu >= 0._wp ) THEN      ! case 1 
     270                  zmaxu = ze3wu / fse3w(ji+1,jj,iku+1) 
     271                  ! interpolated values of tracers 
     272                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) ) 
     273                  ! gradient of tracers 
     274                  sgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     275                  sgzu(ji,jj)    = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku) 
     276               ELSE                           ! case 2 
     277                  zmaxu = - ze3wu / fse3w(ji,jj,iku+1) 
     278                  ! interpolated values of tracers 
     279                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) ) 
     280                  sgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     281                  ! gradient of  tracers 
     282                  sgzu(ji,jj)    = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu) 
     283               ENDIF 
     284               ! 
     285               ! j- direction 
     286               IF( ze3wv >= 0._wp ) THEN      ! case 1 
     287                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv+1) 
     288                  ! interpolated values of tracers 
     289                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) ) 
     290                  ! gradient of tracers 
     291                  sgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     292                  sgzv(ji,jj)    = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv)  
     293               ELSE                           ! case 2 
     294                  zmaxv =  - ze3wv / fse3w(ji,jj,ikv+1) 
     295                  ! interpolated values of tracers 
     296                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) ) 
     297                  ! gradient of tracers 
     298                  sgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     299                  sgzv(ji,jj)    = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv) 
     300               ENDIF 
     301# if ! defined key_vectopt_loop 
     302            END DO 
     303# endif 
     304         END DO!! 
     305         CALL lbc_lnk( sgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( sgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     306         ! 
     307      END DO 
     308 
     309      ! horizontal derivative of density anomalies (rd) 
     310      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
     311# if defined key_vectopt_loop 
     312         jj = 1 
     313         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled) 
     314# else 
     315         DO jj = 1, jpjm1 
     316            DO ji = 1, jpim1 
     317# endif 
     318               iku = miku(ji,jj) 
     319               ikv = mikv(ji,jj) 
     320               ze3wu  = fse3w(ji+1,jj  ,iku+1) - fse3w(ji,jj,iku+1) 
     321               ze3wv  = fse3w(ji  ,jj+1,ikv+1) - fse3w(ji,jj,ikv+1) 
     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# if ! defined key_vectopt_loop 
     330            END DO 
     331# endif 
     332         END DO 
     333 
     334         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
     335         ! step and store it in  zri, zrj for each  case 
     336         CALL eos( zti, zhi, zri )   
     337         CALL eos( ztj, zhj, zrj ) 
     338 
     339         ! Gradient of density at the last level  
     340# if defined key_vectopt_loop 
     341         jj = 1 
     342         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled) 
     343# else 
     344         DO jj = 1, jpjm1 
     345            DO ji = 1, jpim1 
     346# endif 
     347               iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1 
     348               ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1 
     349               ze3wu  = fse3w(ji+1,jj  ,iku+1) - fse3w(ji,jj,iku+1) 
     350               ze3wv  = fse3w(ji  ,jj+1,ikv+1) - fse3w(ji,jj,ikv+1) 
     351               IF( ze3wu >= 0._wp ) THEN 
     352                 sgru  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1 
     353                 smru  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1  
     354                 sge3ru(ji,jj) = umask(ji,jj,iku+1)                                                                  & 
     355                                * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   & 
     356                                   - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1 
     357               ELSE 
     358                 sgru  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) - zri(ji,jj) )      ! i: 2 
     359                 smru  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2 
     360                 sge3ru(ji,jj) = umask(ji,jj,iku+1)                                                                   & 
     361                                * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  & 
     362                                   -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2 
     363               ENDIF 
     364               IF( ze3wv >= 0._wp ) THEN 
     365                 sgrv  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )        ! j: 1 
     366                 smrv  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1 
     367                 sge3rv(ji,jj) = vmask(ji,jj,ikv+1)                                                                  &  
     368                                * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  & 
     369                                   - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1 
     370                                  ! + 2 due to the formulation in density and not in anomalie in hpg sco 
     371               ELSE 
     372                 sgrv  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )     ! j: 2 
     373                 smrv  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2 
     374                 sge3rv(ji,jj) = vmask(ji,jj,ikv+1)                                                                   & 
     375                                * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) & 
     376                                   -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2 
     377               ENDIF 
     378# if ! defined key_vectopt_loop 
     379            END DO 
     380# endif 
     381         END DO 
     382         CALL lbc_lnk( sgru   , 'U', -1. )   ;   CALL lbc_lnk( sgrv   , 'V', -1. )   ! Lateral boundary conditions 
     383         CALL lbc_lnk( smru   , 'U',  1. )   ;   CALL lbc_lnk( smrv   , 'V',  1. )   ! Lateral boundary conditions 
     384         CALL lbc_lnk( sgzu   , 'U', -1. )   ;   CALL lbc_lnk( sgzv   , 'V', -1. )   ! Lateral boundary conditions 
     385         CALL lbc_lnk( sge3ru , 'U', -1. )   ;   CALL lbc_lnk( sge3rv , 'V', -1. )   ! Lateral boundary conditions 
     386         ! 
     387      END IF   
     388      ! 
     389      CALL wrk_dealloc( jpi, jpj,       zri, zrj, zhi, zhj)  
    212390      CALL wrk_dealloc( jpi, jpj, kjpt, zti, ztj           )  
    213391      ! 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90

    r4147 r4666  
    4040   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   avtb_2d        !: horizontal shape of background Kz profile 
    4141   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   bfrua, bfrva   !: Bottom friction coefficients set in zdfbfr 
     42   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   tfrua, tfrva   !: top friction coefficients set in zdfbfr 
    4243   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avmu , avmv    !: vertical viscosity coef at uw- & vw-pts       [m2/s] 
    4344   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avm  , avt     !: vertical viscosity & diffusivity coef at w-pt [m2/s] 
    … …  
    5758      ALLOCATE(avmb(jpk) , bfrua(jpi,jpj) ,                         & 
    5859         &     avtb(jpk) , bfrva(jpi,jpj) , avtb_2d(jpi,jpj) ,      & 
     60         &     tfrua(jpi, jpj), tfrva(jpi, jpj)              ,      & 
    5961         &     avmu(jpi,jpj,jpk), avm(jpi,jpj,jpk)           ,      & 
    6062         &     avmv(jpi,jpj,jpk), avt(jpi,jpj,jpk)           , STAT = zdf_oce_alloc ) 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r4624 r4666  
    4141   REAL(wp), PUBLIC ::   rn_bfrien    ! local factor to enhance coefficient bfri (PUBLIC for TAM) 
    4242   LOGICAL , PUBLIC ::   ln_bfr2d     ! logical switch for 2D enhancement (PUBLIC for TAM) 
     43   REAL(wp), PUBLIC ::   rn_tfri1     ! top drag coefficient (linear case)  (PUBLIC for TAM) 
     44   REAL(wp), PUBLIC ::   rn_tfri2     ! top drag coefficient (non linear case) (PUBLIC for TAM) 
     45   REAL(wp), PUBLIC ::   rn_tfri2_max ! Maximum top drag coefficient (non linear case and ln_loglayer=T) (PUBLIC for TAM) 
     46   REAL(wp), PUBLIC ::   rn_tfeb2     ! background top turbulent kinetic energy  [m2/s2] (PUBLIC for TAM) 
     47   REAL(wp), PUBLIC ::   rn_tfrien    ! local factor to enhance coefficient tfri (PUBLIC for TAM) 
     48   LOGICAL , PUBLIC ::   ln_tfr2d     ! logical switch for 2D enhancement (PUBLIC for TAM) 
     49 
    4350   LOGICAL , PUBLIC ::   ln_loglayer  ! switch for log layer bfr coeff. (PUBLIC for TAM) 
    4451   REAL(wp), PUBLIC ::   rn_bfrz0     ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM) 
     52   REAL(wp), PUBLIC ::   rn_tfrz0     ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM) 
    4553   LOGICAL , PUBLIC ::   ln_bfrimp    ! logical switch for implicit bottom friction 
    46    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d            ! 2D bottom drag coefficient (PUBLIC for TAM) 
     54   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d, tfrcoef2d   ! 2D bottom/top drag coefficient (PUBLIC for TAM) 
    4755 
    4856   !! * Substitutions 
    … …  
    6068      !!                ***  FUNCTION zdf_bfr_alloc  *** 
    6169      !!---------------------------------------------------------------------- 
    62       ALLOCATE( bfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc ) 
     70      ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc ) 
    6371      ! 
    6472      IF( lk_mpp             )   CALL mpp_sum ( zdf_bfr_alloc ) 
    … …  
    8896      INTEGER  ::   ikbt, ikbu, ikbv             ! local integers 
    8997      REAL(wp) ::   zvu, zuv, zecu, zecv, ztmp   ! temporary scalars 
    90       REAL(wp), POINTER, DIMENSION(:,:) ::  zbfrt 
     98      REAL(wp), POINTER, DIMENSION(:,:) ::  zbfrt, ztfrt 
    9199      !!---------------------------------------------------------------------- 
    92100      ! 
    … …  
    101109      IF( nn_bfr == 2 ) THEN                 ! quadratic bottom friction only 
    102110         ! 
    103          CALL wrk_alloc( jpi, jpj, zbfrt ) 
     111         CALL wrk_alloc( jpi, jpj, zbfrt, ztfrt ) 
    104112 
    105113         IF ( ln_loglayer.AND.lk_vvl ) THEN ! "log layer" bottom friction coefficient 
    … …  
    120128                  zbfrt(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp) 
    121129                  zbfrt(ji,jj) = MIN(zbfrt(ji,jj), rn_bfri2_max) 
     130! (ISF) 
     131                  ikbt = mikt(ji,jj) 
     132! JC: possible WAD implementation should modify line below if layers vanish 
     133                  ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp 
     134                  ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 
     135                  ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max) 
     136 
    122137               END DO 
    123138            END DO 
    … …  
    125140         ELSE 
    126141            zbfrt(:,:) = bfrcoef2d(:,:) 
     142            ztfrt(:,:) = tfrcoef2d(:,:) 
    127143         ENDIF 
    128144 
    … …  
    150166               bfrua(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) * zecu 
    151167               bfrva(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) * zecv 
     168               ! 
     169               ! (ISF) ======================================================================== 
     170               ikbu = miku(ji,jj)         ! ocean bottom level at u- and v-points  
     171               ikbv = mikv(ji,jj)         ! (deepest ocean u- and v-points) 
     172               ! 
     173               zvu  = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     & 
     174                  &           + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  ) 
     175               zuv  = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     & 
     176                  &           + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  ) 
     177               ! 
     178               zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 ) 
     179               zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 ) 
     180               ! 
     181               tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj  ) ) * zecu  
     182               tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji  ,jj+1) ) * zecv 
     183               ! (ISF) END ==================================================================== 
     184 
    152185            END DO 
    153186         END DO 
    … …  
    158191         IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        & 
    159192            &                       tab2d_2=bfrva, clinfo2=       ' v: ', mask2=vmask,ovlap=1 ) 
    160          CALL wrk_dealloc( jpi,jpj,zbfrt ) 
     193         CALL wrk_dealloc( jpi,jpj,zbfrt, ztfrt ) 
    161194      ENDIF 
    162195      ! 
    … …  
    183216      INTEGER   ::   ios 
    184217      REAL(wp)  ::   zminbfr, zmaxbfr   ! temporary scalars 
     218      REAL(wp)  ::   zmintfr, zmaxtfr   ! temporary scalars 
    185219      REAL(wp)  ::   ztmp, zfru, zfrv   !    -         - 
    186220      !! 
    187221      NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfri2_max, rn_bfeb2, rn_bfrz0, ln_bfr2d, & 
    188                     &  rn_bfrien, ln_bfrimp, ln_loglayer 
     222                    &          rn_tfri1, rn_tfri2, rn_tfri2_max, rn_tfeb2, rn_tfrz0, ln_tfr2d, & 
     223                    &  rn_bfrien, rn_tfrien, ln_bfrimp, ln_loglayer 
    189224      !!---------------------------------------------------------------------- 
    190225      ! 
    … …  
    215250         bfrua(:,:) = 0.e0 
    216251         bfrva(:,:) = 0.e0 
     252         tfrua(:,:) = 0.e0 
     253         tfrva(:,:) = 0.e0 
    217254         ! 
    218255      CASE( 1 ) 
    219256         IF(lwp) WRITE(numout,*) '      linear botton friction' 
    220          IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri1  = ', rn_bfri1 
     257         IF(lwp) WRITE(numout,*) '      bottom friction coef.   rn_bfri1  = ', rn_bfri1 
    221258         IF( ln_bfr2d ) THEN 
    222259            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d 
    223260            IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien 
     261         ENDIF 
     262         IF(lwp) WRITE(numout,*) '      top    friction coef.   rn_bfri1  = ', rn_bfri1 
     263         IF( ln_tfr2d ) THEN 
     264            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d 
     265            IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien 
    224266         ENDIF 
    225267         ! 
    … …  
    236278         bfrua(:,:) = - bfrcoef2d(:,:) 
    237279         bfrva(:,:) = - bfrcoef2d(:,:) 
     280         ! 
     281         IF(ln_tfr2d) THEN 
     282            ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 
     283            CALL iom_open('tfr_coef.nc',inum) 
     284            CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array 
     285            CALL iom_close(inum) 
     286            tfrcoef2d(:,:) = rn_tfri1 * ( 1 + rn_tfrien * tfrcoef2d(:,:) ) 
     287         ELSE 
     288            tfrcoef2d(:,:) = rn_tfri1  ! initialize tfrcoef2d to the namelist variable 
     289         ENDIF 
     290         ! 
     291         tfrua(:,:) = - tfrcoef2d(:,:) 
     292         tfrva(:,:) = - tfrcoef2d(:,:) 
    238293         ! 
    239294      CASE( 2 ) 
    … …  
    252307            IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien 
    253308         ENDIF 
     309         IF(lwp) WRITE(numout,*) '      quadratic top    friction' 
     310         IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri2  = ', rn_tfri2 
     311         IF(lwp) WRITE(numout,*) '      Max. coef. (log case)   rn_tfri2_max  = ', rn_tfri2_max 
     312         IF(lwp) WRITE(numout,*) '      background tke   rn_tfeb2  = ', rn_tfeb2 
     313         IF(lwp) WRITE(numout,*) '      log formulation   ln_tfr2d = ', ln_loglayer 
     314         IF(lwp) WRITE(numout,*) '      bottom roughness  rn_tfrz0 [m] = ', rn_tfrz0 
     315         IF( rn_tfrz0<=0.e0 ) THEN 
     316            WRITE(ctmp1,*) '      bottom roughness must be strictly positive' 
     317            CALL ctl_stop( ctmp1 ) 
     318         ENDIF 
     319         IF( ln_tfr2d ) THEN 
     320            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_tfr2d  = ', ln_tfr2d 
     321            IF(lwp) WRITE(numout,*) '      coef rn_tfri2 enhancement factor                rn_tfrien  = ',rn_tfrien 
     322         ENDIF 
    254323         ! 
    255324         IF(ln_bfr2d) THEN 
    … …  
    263332            bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable 
    264333         ENDIF 
     334 
     335         IF(ln_tfr2d) THEN 
     336            ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 
     337            CALL iom_open('tfr_coef.nc',inum) 
     338            CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! bfrcoef2d is used as tmp array 
     339            CALL iom_close(inum) 
     340            ! 
     341            tfrcoef2d(:,:) = rn_tfri2 * ( 1 + rn_tfrien * tfrcoef2d(:,:) ) 
     342         ELSE 
     343            tfrcoef2d(:,:) = rn_tfri2  ! initialize tfrcoef2d to the namelist variable 
     344         ENDIF 
    265345         ! 
    266346         IF ( ln_loglayer.AND.(.NOT.lk_vvl) ) THEN ! set "log layer" bottom friction once for all 
    … …  
    279359                  bfrcoef2d(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp) 
    280360                  bfrcoef2d(ji,jj) = MIN(bfrcoef2d(ji,jj), rn_bfri2_max) 
     361                  ikbt = mikt(ji,jj) 
     362                  ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_tfrz0 ))**2._wp 
     363                  tfrcoef2d(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 
     364                  tfrcoef2d(ji,jj) = MIN(tfrcoef2d(ji,jj), rn_tfri2_max) 
    281365               END DO 
    282366            END DO 
    … …  
    308392      zminbfr =  1.e10_wp    ! initialise tracker for minimum of bottom friction coefficient 
    309393      zmaxbfr = -1.e10_wp    ! initialise tracker for maximum of bottom friction coefficient 
     394      zmintfr =  1.e10_wp    ! initialise tracker for minimum of bottom friction coefficient 
     395      zmaxtfr = -1.e10_wp    ! initialise tracker for maximum of bottom friction coefficient 
    310396      ! 
    311397#  if defined key_vectopt_loop 
    … …  
    339425             zminbfr = MIN(  zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) )  ) 
    340426             zmaxbfr = MAX(  zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) )  ) 
     427! (ISF) 
     428             ikbu = miku(ji,jj)       ! deepest ocean level at u- and v-points 
     429             ikbv = mikv(ji,jj) 
     430             zfru = 0.5 * fse3u(ji,jj,ikbu) / rdt 
     431             zfrv = 0.5 * fse3v(ji,jj,ikbv) / rdt 
     432             IF( ABS( tfrcoef2d(ji,jj) ) > zfru ) THEN 
     433                IF( ln_ctl ) THEN 
     434                   WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbu 
     435                   WRITE(numout,*) 'BFR ', ABS( tfrcoef2d(ji,jj) ), zfru 
     436                ENDIF 
     437                ictu = ictu + 1 
     438             ENDIF 
     439             IF( ABS( tfrcoef2d(ji,jj) ) > zfrv ) THEN 
     440                 IF( ln_ctl ) THEN 
     441                     WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbv 
     442                     WRITE(numout,*) 'BFR ', tfrcoef2d(ji,jj), zfrv 
     443                 ENDIF 
     444                 ictv = ictv + 1 
     445             ENDIF 
     446             zmintfr = MIN(  zmintfr, MIN( zfru, ABS( tfrcoef2d(ji,jj) ) )  ) 
     447             zmaxtfr = MAX(  zmaxtfr, MIN( zfrv, ABS( tfrcoef2d(ji,jj) ) )  ) 
     448 
    341449         END DO 
    342450      END DO 
    … …  
    346454         CALL mpp_min( zminbfr ) 
    347455         CALL mpp_max( zmaxbfr ) 
     456         CALL mpp_min( zmintfr ) 
     457         CALL mpp_max( zmaxtfr ) 
    348458      ENDIF 
    349459      IF( .NOT.ln_bfrimp) THEN 
    … …  
    352462         WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points ' 
    353463         WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr 
     464         WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr 
    354465         WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary' 
    355466      ENDIF 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

    r4624 r4666  
    107107 
    108108      !                                                ! =============== 
    109       DO jk = 2, jpkm1                                 ! Horizontal slab 
     109      !                                                ! Horizontal slab 
    110110         !                                             ! =============== 
     111      DO jj = 1, jpj                                     ! indicators: 
     112         DO ji = 1, jpi 
     113            DO jk = mikt(ji,jj)+1, jpkm1                                 ! Horizontal slab 
    111114         ! Define the mask  
    112115         ! --------------- 
    113          rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) )         ! only retains positive value of rrau 
    114  
    115          DO jj = 1, jpj                                     ! indicators: 
    116             DO ji = 1, jpi 
     116               rrau(ji,jj,jk) = MAX( 1.e-20, rrau(ji,jj,jk) )         ! only retains positive value of rrau 
     117 
    117118               ! stability indicator: msks=1 if rn2>0; 0 elsewhere 
    118119               IF( rn2(ji,jj,jk) + 1.e-12  <= 0. ) THEN   ;   zmsks(ji,jj) = 0._wp 
    … …  
    136137               ELSE                                                         ;   zmskd3(ji,jj) = 1._wp 
    137138               ENDIF 
    138             END DO 
    139          END DO 
    140139         ! mask zmsk in order to have avt and avs masked 
    141          zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk) 
     140               zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk) 
    142141 
    143142 
    … …  
    145144         ! ------------------ 
    146145         ! Constant eddy coefficient: reset to the background value 
    147 !CDIR NOVERRCHK 
    148          DO jj = 1, jpj 
    149 !CDIR NOVERRCHK 
    150             DO ji = 1, jpi 
    151146               zinr = 1./rrau(ji,jj,jk) 
    152147               ! salt fingering 
    … …  
    165160            END DO 
    166161         END DO 
    167  
    168  
     162      END DO 
    169163         ! Increase avmu, avmv if necessary 
    170164         ! -------------------------------- 
    171165!!gm to be changed following the definition of avm. 
    172          DO jj = 1, jpjm1 
    173             DO ji = 1, fs_jpim1   ! vector opt. 
     166      DO jj = 1, jpjm1 
     167         DO ji = 1, fs_jpim1   ! vector opt. 
     168            DO jk = miku(ji,jj)+1, jpkm1                                 ! Horizontal slab 
    174169               avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk),    & 
    175170                  &                  avt(ji,jj,jk), avt(ji+1,jj,jk),   & 
    176171                  &                  avs(ji,jj,jk), avs(ji+1,jj,jk) )  * umask(ji,jj,jk) 
     172            END DO 
     173            DO jk = mikv(ji,jj)+1, jpkm1 
    177174               avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk),    & 
    178175                  &                  avt(ji,jj,jk), avt(ji,jj+1,jk),   & 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

    r4624 r4666  
    241241      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0) 
    242242         DO ji = fs_2, fs_jpim1   ! vector opt. 
    243             en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
     243            IF (mikt(ji,jj) .GT. 1) THEN 
     244               en(ji,jj,mikt(ji,jj))=rn_emin 
     245            ELSE 
     246               en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 
     247            END IF 
    244248         END DO 
    245249      END DO 
    … …  
    342346      END DO 
    343347      ! 
    344       DO jk = 2, jpkm1           !* Matrix and right hand side in en 
    345          DO jj = 2, jpjm1 
    346             DO ji = fs_2, fs_jpim1   ! vector opt. 
     348      DO jj = 2, jpjm1 
     349         DO ji = fs_2, fs_jpim1   ! vector opt. 
     350            DO jk = mikt(ji,jj)+1, jpkm1           !* Matrix and right hand side in en 
    347351               zcof   = zfact1 * tmask(ji,jj,jk) 
    348352               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal 
    … …  
    360364               !                                   ! right hand side in en 
    361365               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    & 
    362                   &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) * tmask(ji,jj,jk) 
    363             END DO 
    364          END DO 
    365       END DO 
    366       !                          !* Matrix inversion from level 2 (tke prescribed at level 1) 
    367       DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    368          DO jj = 2, jpjm1 
    369             DO ji = fs_2, fs_jpim1    ! vector opt. 
     366                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) & 
     367                  &                                 * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 
     368            END DO 
     369            !                          !* Matrix inversion from level 2 (tke prescribed at level 1) 
     370            DO jk = mikt(ji,jj)+2, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    370371               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 
    371372            END DO 
    372          END DO 
    373       END DO 
    374       DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    375          DO ji = fs_2, fs_jpim1    ! vector opt. 
    376             zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke 
    377          END DO 
    378       END DO 
    379       DO jk = 3, jpkm1 
    380          DO jj = 2, jpjm1 
    381             DO ji = fs_2, fs_jpim1    ! vector opt. 
     373            ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
     374            zd_lw(ji,jj,mikt(ji,jj)+1) = en(ji,jj,mikt(ji,jj)+1) - zd_lw(ji,jj,mikt(ji,jj)+1) * en(ji,jj,mikt(ji,jj))    ! Surface boudary conditions on tke 
     375            ! 
     376            DO jk = mikt(ji,jj)+2, jpkm1 
    382377               zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1) 
    383378            END DO 
    384          END DO 
    385       END DO 
    386       DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    387          DO ji = fs_2, fs_jpim1    ! vector opt. 
     379            ! 
     380            ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    388381            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 
    389          END DO 
    390       END DO 
    391       DO jk = jpk-2, 2, -1 
    392          DO jj = 2, jpjm1 
    393             DO ji = fs_2, fs_jpim1    ! vector opt. 
     382            ! 
     383            DO jk = jpk-2, mikt(ji,jj)+1, -1 
    394384               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 
    395385            END DO 
    396          END DO 
    397       END DO 
    398       DO jk = 2, jpkm1                             ! set the minimum value of tke 
    399          DO jj = 2, jpjm1 
    400             DO ji = fs_2, fs_jpim1   ! vector opt. 
     386            ! 
     387            DO jk = mikt(ji,jj), jpkm1                             ! set the minimum value of tke 
    401388               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk) 
    402389            END DO 
    … …  
    412399               DO ji = fs_2, fs_jpim1   ! vector opt. 
    413400                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    414                      &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) 
     401                     &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) * tmask(ji,jj,1) 
    415402               END DO 
    416403            END DO 
    … …  
    421408               jk = nmln(ji,jj) 
    422409               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    423                   &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) 
     410                  &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk) * tmask(ji,jj,1) 
    424411            END DO 
    425412         END DO 
    … …  
    433420                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj) 
    434421                  zty2 = vtau(ji  ,jj-1) + vtau(ji,jj) 
    435                   ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )    ! module of the mean stress  
     422                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress  
    436423                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean  
    437424                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications... 
    438425                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   & 
    439                      &                                        * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) 
     426                     &                                        * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * tmask(ji,jj,1) 
    440427               END DO 
    441428            END DO 
    … …  
    506493      ! 
    507494      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 
    508          zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
    509          zmxlm(:,:,1) = MAX(  rn_mxl0,  zraug * taum(:,:)  ) 
    510       ELSE                          ! surface set to the minimum value 
    511          zmxlm(:,:,1) = rn_mxl0 
     495         DO jj = 2, jpjm1 
     496            DO ji = fs_2, fs_jpim1 
     497               IF (mikt(ji,jj) .GT. 1) THEN 
     498                  zmxlm(ji,jj,mikt(ji,jj)) = rmxl_min 
     499               ELSE 
     500                  zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 
     501                  zmxlm(ji,jj,mikt(ji,jj)) = MAX( rn_mxl0, zraug * taum(ji,jj) ) 
     502               END IF 
     503            END DO 
     504         END DO 
     505      ELSE  
     506         DO jj = 2, jpjm1 
     507            DO ji = fs_2, fs_jpim1                         ! surface set to the minimum value 
     508               zmxlm(ji,jj,mikt(ji,jj)) = MAX( tmask(ji,jj,1) * rn_mxl0, rmxl_min) 
     509            END DO 
     510         END DO 
    512511      ENDIF 
    513512      zmxlm(:,:,jpk)  = rmxl_min     ! last level set to the interior minium value 
    514513      ! 
    515514!CDIR NOVERRCHK 
    516       DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2) 
    517 !CDIR NOVERRCHK 
    518          DO jj = 2, jpjm1 
    519 !CDIR NOVERRCHK 
    520             DO ji = fs_2, fs_jpim1   ! vector opt. 
     515      DO jj = 2, jpjm1 
     516!CDIR NOVERRCHK 
     517         DO ji = fs_2, fs_jpim1   ! vector opt. 
     518            !CDIR NOVERRCHK 
     519            DO jk = mikt(ji,jj)+1, jpkm1              ! interior value : l=sqrt(2*e/n^2) 
    521520               zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    522521               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  ) 
    523522            END DO 
     523            zmxld(ji,jj,mikt(ji,jj)) = zmxlm(ji,jj,mikt(ji,jj))   ! surface set to the minimum value  
    524524         END DO 
    525525      END DO 
    … …  
    533533      ! 
    534534      CASE ( 0 )           ! bounded by the distance to surface and bottom 
    535          DO jk = 2, jpkm1 
    536             DO jj = 2, jpjm1 
    537                DO ji = fs_2, fs_jpim1   ! vector opt. 
    538                   zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk),   & 
     535         DO jj = 2, jpjm1 
     536            DO ji = fs_2, fs_jpim1   ! vector opt. 
     537               DO jk = mikt(ji,jj)+1, jpkm1 
     538                  zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   & 
    539539                  &            fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 
    540540                  zmxlm(ji,jj,jk) = zemxl 
    … …  
    545545         ! 
    546546      CASE ( 1 )           ! bounded by the vertical scale factor 
    547          DO jk = 2, jpkm1 
    548             DO jj = 2, jpjm1 
    549                DO ji = fs_2, fs_jpim1   ! vector opt. 
     547         DO jj = 2, jpjm1 
     548            DO ji = fs_2, fs_jpim1   ! vector opt. 
     549               DO jk = mikt(ji,jj)+1, jpkm1 
    550550                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 
    551551                  zmxlm(ji,jj,jk) = zemxl 
    … …  
    556556         ! 
    557557      CASE ( 2 )           ! |dk[xml]| bounded by e3t : 
    558          DO jk = 2, jpkm1         ! from the surface to the bottom : 
    559             DO jj = 2, jpjm1 
    560                DO ji = fs_2, fs_jpim1   ! vector opt. 
     558         DO jj = 2, jpjm1 
     559            DO ji = fs_2, fs_jpim1   ! vector opt. 
     560               DO jk = mikt(ji,jj)+1, jpkm1         ! from the surface to the bottom : 
    561561                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    562562               END DO 
    563             END DO 
    564          END DO 
    565          DO jk = jpkm1, 2, -1     ! from the bottom to the surface : 
    566             DO jj = 2, jpjm1 
    567                DO ji = fs_2, fs_jpim1   ! vector opt. 
     563               DO jk = jpkm1, mikt(ji,jj)+1, -1     ! from the bottom to the surface : 
    568564                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    569565                  zmxlm(ji,jj,jk) = zemxl 
    … …  
    574570         ! 
    575571      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t : 
    576          DO jk = 2, jpkm1         ! from the surface to the bottom : lup 
    577             DO jj = 2, jpjm1 
    578                DO ji = fs_2, fs_jpim1   ! vector opt. 
     572         DO jj = 2, jpjm1 
     573            DO ji = fs_2, fs_jpim1   ! vector opt. 
     574               DO jk = mikt(ji,jj)+1, jpkm1         ! from the surface to the bottom : lup 
    579575                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 
    580576               END DO 
    581             END DO 
    582          END DO 
    583          DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown 
    584             DO jj = 2, jpjm1 
    585                DO ji = fs_2, fs_jpim1   ! vector opt. 
     577               DO jk = jpkm1, mikt(ji,jj)+1, -1     ! from the bottom to the surface : ldown 
    586578                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 
    587579               END DO 
    … …  
    628620      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged) 
    629621      ! 
    630       DO jk = 2, jpkm1            !* vertical eddy viscosity at u- and v-points 
     622      DO jj = 2, jpjm1 
     623         DO ji = fs_2, fs_jpim1   ! vector opt. 
     624            DO jk = miku(ji,jj)+1, jpkm1            !* vertical eddy viscosity at u- and v-points 
     625               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk) 
     626            END DO 
     627            DO jk = mikv(ji,jj)+1, jpkm1 
     628               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk) 
     629            END DO 
     630         END DO 
     631      END DO 
     632      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions 
     633      ! 
     634      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt 
    631635         DO jj = 2, jpjm1 
    632636            DO ji = fs_2, fs_jpim1   ! vector opt. 
    633                avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk) 
    634                avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk) 
    635             END DO 
    636          END DO 
    637       END DO 
    638       CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions 
    639       ! 
    640       IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt 
    641          DO jk = 2, jpkm1 
    642             DO jj = 2, jpjm1 
    643                DO ji = fs_2, fs_jpim1   ! vector opt. 
     637               DO jk = mikt(ji,jj)+1, jpkm1 
    644638                  zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 
    645639                  !                                          ! shear 
    … …  
    655649                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
    656650# if defined key_c1d 
    657                   e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number 
    658                   e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri 
     651                  e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)  ! c1d configuration : save masked Prandlt number 
     652                  e_ric(ji,jj,jk) = zri   * tmask(ji,jj,jk)  ! c1d config. : save Ri 
    659653# endif 
    660654              END DO 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90

    r4624 r4666  
    126126      zkz(:,:) = 0.e0               !* Associated potential energy consummed over the whole water column 
    127127      DO jk = 2, jpkm1 
    128          zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk)* tmask(:,:,jk) 
     128         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk)* tmask(:,:,jk) * tmask(:,:,jk-1) 
    129129      END DO 
    130130 
    … …  
    135135      END DO 
    136136 
    137       DO jk = 2, jpkm1              !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s 
    138          zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. )   !kz max = 300 cm2/s 
     137      DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 
     138         DO ji = 1, jpi 
     139            DO jk = mikt(ji,jj)+1, jpkm1     !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s 
     140               zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. )   !kz max = 300 cm2/s 
     141            END DO 
     142         END DO 
    139143      END DO 
    140144 
    … …  
    162166      !                          !   Update  mixing coefs  !                           
    163167      !                          ! ----------------------- ! 
    164       DO jk = 2, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
    165          avt(:,:,jk) = avt(:,:,jk) + zav_tide(:,:,jk) 
    166          avm(:,:,jk) = avm(:,:,jk) + zav_tide(:,:,jk) 
    167          DO jj = 2, jpjm1 
    168             DO ji = fs_2, fs_jpim1  ! vector opt. 
     168      DO jj = 1, jpj                !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 
     169         DO ji = 1, jpi 
     170            DO jk = mikt(ji,jj)+1, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
     171               avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk) 
     172               avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk) 
     173            END DO 
     174         END DO 
     175      END DO 
     176       
     177      DO jj = 2, jpjm1 
     178         DO ji = fs_2, fs_jpim1  ! vector opt. 
     179            DO jk = mikt(ji,jj)+1, jpkm1              !* update momentum & tracer diffusivity with tidal mixing 
    169180               avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj  ,jk) ) * umask(ji,jj,jk) 
    170181               avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji  ,jj+1,jk) ) * vmask(ji,jj,jk) 
    … …  
    237248      zsum2(:,:) = 0.e0 
    238249      DO jk= 2, jpk 
    239          zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk) 
    240          zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk)                 
     250         zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1) 
     251         zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)                
    241252      END DO 
    242253      DO jj = 1, jpj 
    … …  
    274285      zkz(:,:) = 0.e0               ! Associated potential energy consummed over the whole water column 
    275286      DO jk = 2, jpkm1 
    276          zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * tmask(:,:,jk) 
     287         zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1) 
    277288      END DO 
    278289 
    … …  
    284295 
    285296      DO jk = 2, jpkm1              ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zavt_itf bound by 300 cm2/s 
    286          zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. )   ! kz max = 120 cm2/s 
     297         zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. ) * tmask(:,:,jk) * tmask(:,:,jk-1)   ! kz max = 120 cm2/s 
    287298      END DO 
    288299 
    … …  
    377388      READ  ( numnam_cfg, namzdf_tmx, IOSTAT = ios, ERR = 902 ) 
    378389902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist', lwp ) 
    379       IF(lwm) WRITE ( numond, namzdf_tmx ) 
     390      WRITE ( numond, namzdf_tmx ) 
    380391 
    381392      IF(lwp) THEN                   ! Control print 
    … …  
    414425      ! only the energy available for mixing is taken into account, 
    415426      ! (mixing efficiency tidal dissipation efficiency) 
    416       en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * tmask(:,:,1) 
     427      en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * lmask(:,:) 
    417428 
    418429      ! Vertical structure (az_tmx) 
    … …  
    443454         ztpc = 0.e0 
    444455         zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:) 
    445          DO jk= 2, jpkm1 
    446             DO jj = 1, jpj 
    447                DO ji = 1, jpi 
     456         DO jj = 1, jpj 
     457            DO ji = 1, jpi 
     458               DO jk= mikt(ji,jj)+1, jpkm1 
    448459                  ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 
    449460               END DO 
    … …  
    459470         zav_tide(:,:,:) = MIN( zav_tide(:,:,:), 60.e-4 )    
    460471         zkz(:,:) = 0.e0 
    461          DO jk = 2, jpkm1 
    462          DO jj = 1, jpj 
    463             DO ji = 1, jpi 
     472         DO jj = 1, jpj 
     473            DO ji = 1, jpi 
     474               DO jk = mikt(ji,jj)+1, jpkm1 
    464475               zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zav_tide(ji,jj,jk)* tmask(ji,jj,jk) 
    465             END DO 
    466          END DO 
     476               END DO 
     477            END DO 
    467478         END DO 
    468479         ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz 
    … …  
    484495         WRITE(numout,*) '          Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) ) 
    485496 
    486          DO jk = 2, jpkm1 
    487             zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. )   !kz max = 300 cm2/s 
     497         DO jj = 1, jpj 
     498            DO ji = 1, jpi 
     499               DO jk = mikt(ji,jj)+1, jpkm1 
     500                  zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. )   !kz max = 300 cm2/s 
     501               END DO 
     502            END DO 
    488503         END DO 
    489504         ztpc = 0.e0 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/oce.F90

    r4354 r4666  
    4646   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gtsu, gtsv   !: horizontal gradient of T, S bottom u-point  
    4747   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gru , grv    !: horizontal gradient of rd at bottom u-point 
     48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aru , arv     
     49   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gzu , gzv     
     50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ge3ru, ge3rv   !: horizontal gradient of T, S and rd at top v-point   
     51 
     52   !! (ISF) interpolated gradient (only used for ice shelf case)  
     53   !! ---------------------  
     54   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gtui, gtvi   !: horizontal gradient of T, S and rd at top u-point  
     55   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   grui, grvi   !: horizontal gradient of T, S and rd at top v-point   
     56   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   arui, arvi   !: horizontal average  of rd          at top v-point   
     57   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   gzui, gzvi   !: horizontal gradient of z           at top v-point   
     58   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ge3rui, ge3rvi   !: horizontal gradient of T, S and rd at top v-point   
    4859 
    4960   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rke          !: kinetic energy 
    … …  
    90101         &     spgu  (jpi,jpj)   , spgv(jpi,jpj)   ,                       & 
    91102         &     gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts),                     & 
    92          &     gru(jpi,jpj)      , grv(jpi,jpj)                      , STAT=ierr(2) ) 
     103         &     aru(jpi,jpj)      , arv(jpi,jpj)      ,                     & 
     104         &     gzu(jpi,jpj)      , gzv(jpi,jpj)      ,                     & 
     105         &     gru(jpi,jpj)      , grv(jpi,jpj)      ,                     & 
     106         &     ge3ru(jpi,jpj)    , ge3rv(jpi,jpj)    ,                     & 
     107         &     gtui(jpi,jpj,jpts), gtvi(jpi,jpj,jpts),                     & 
     108         &     arui(jpi,jpj)     , arvi(jpi,jpj)     ,                     & 
     109         &     gzui(jpi,jpj)     , gzvi(jpi,jpj)     ,                     & 
     110         &     ge3rui(jpi,jpj)   , ge3rvi(jpi,jpj)   ,                     & 
     111         &     grui(jpi,jpj)     , grvi(jpi,jpj)     ,                           STAT=ierr(2) ) 
    93112         ! 
    94113      ALLOCATE( snwice_mass(jpi,jpj) , snwice_mass_b(jpi,jpj), snwice_fmass(jpi,jpj) , STAT=ierr(3) ) 

  • branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/step.F90

    r4624 r4666  
    3131   !!---------------------------------------------------------------------- 
    3232   USE step_oce         ! time stepping definition modules 
     33   USE iom 
    3334 
    3435   IMPLICIT NONE 
    … …  
    141142      IF( lk_ldfslp ) THEN                            ! slope of lateral mixing 
    142143                         CALL eos( tsb, rhd, gdept_0(:,:,:) )             ! before in situ density 
    143          IF( ln_zps )    CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient 
    144             &                                      rhd, gru , grv  )      ! of t, s, rd at the last ocean level 
     144         IF( ln_zps )    CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
     145            &                                      rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             ! 
     146            &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
    145147         IF( ln_traldf_grif ) THEN                           ! before slope for Griffies operator 
    146148                         CALL ldf_slp_grif( kstp ) 
    … …  
    170172          ! Note that the computation of vertical velocity above, hence "after" sea level 
    171173          ! is necessary to compute momentum advection for the rhs of barotropic loop: 
    172                                   CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 
    173           IF( ln_zps      )       CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &   ! zps: now hor. derivative 
    174                 &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level 
     174                            CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 
     175            IF( ln_zps )    CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
     176               &                                      rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             ! 
     177               &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
    175178 
    176179                                  ua(:,:,:) = 0.e0             ! set dynamics trends to zero 
    … …  
    245248                             CALL tra_nxt( kstp )                ! tracer fields at next time step 
    246249                             CALL eos    ( tsa, rhd, rhop, fsdept_n(:,:,:) )  ! Time-filtered in situ density for hpg computation 
    247          IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv,  &    ! zps: time filtered hor. derivative 
    248             &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level 
    249  
     250         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
     251            &                                      rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             ! 
     252            &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
    250253      ELSE                                                  ! centered hpg  (eos then time stepping) 
    251254         IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case 
    252255                                CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation 
    253             IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &    ! zps: now hor. derivative 
    254             &                                          rhd, gru , grv  )      ! of t, s, rd at the last ocean level 
     256            IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient 
     257            &                                      rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             ! 
     258            &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level 
    255259         ENDIF 
    256260         IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection 
Note: See TracChangeset for help on using the changeset viewer.