source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OCE_SRC/SBC/sbcisf.F90 @ 9580

Last change on this file since 9580 was 9570, checked in by nicolasmartin, 3 years ago

Global renaming for core routines (./NEMO)

  • Folders
    • LIM_SRC_3 → ICE_SRC
    • OPA_SRC → OCE_SRC
  • CPP key: key_lim3 → key_si3
  • Modules, (sub)routines and variables names
    • MPI: mpi_comm_opa → mpi_comm_oce, MPI_COMM_OPA → MPI_COMM_OCE, mpi_init_opa → mpi_init_oce
    • AGRIF: agrif_opa_* → agrif_oce_*, agrif_lim3_* → agrif_si3_* and few more
    • TOP-PISCES: p.zlim → p.zice, namp.zlim → namp.zice
  • Comments
    • NEMO/OPA → NEMO/OCE
    • ESIM|LIM3 → SI3
  • Property svn:keywords set to Id
File size: 45.1 KB
Line 
1MODULE sbcisf
2   !!======================================================================
3   !!                       ***  MODULE  sbcisf  ***
4   !! Surface module :  update surface ocean boundary condition under ice
5   !!                   shelf
6   !!======================================================================
7   !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
8   !!            X.X  !  2006-02  (C. Wang   ) Original code bg03
9   !!            3.4  !  2013-03  (P. Mathiot) Merging + parametrization
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   sbc_isf       : update sbc under ice shelf
14   !!----------------------------------------------------------------------
15   USE oce            ! ocean dynamics and tracers
16   USE dom_oce        ! ocean space and time domain
17   USE phycst         ! physical constants
18   USE eosbn2         ! equation of state
19   USE sbc_oce        ! surface boundary condition: ocean fields
20   USE zdfdrg         ! vertical physics: top/bottom drag coef.
21   !
22   USE in_out_manager ! I/O manager
23   USE iom            ! I/O library
24   USE fldread        ! read input field at current time step
25   USE lbclnk         !
26   USE lib_fortran    ! glob_sum
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   sbc_isf, sbc_isf_init, sbc_isf_div, sbc_isf_alloc  ! routine called in sbcmod and divhor
32
33   ! public in order to be able to output then
34
35   REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m]
36   INTEGER , PUBLIC ::   nn_isf                      !: flag to choose between explicit/param/specified 
37   INTEGER , PUBLIC ::   nn_isfblk                   !: flag to choose the bulk formulation to compute the ice shelf melting
38   INTEGER , PUBLIC ::   nn_gammablk                 !: flag to choose how the exchange coefficient is computed
39   REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient []
40   REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient []
41
42   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   misfkt   , misfkb        !: Level of ice shelf base
43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rzisf_tbl                !: depth of calving front (shallowest point) nn_isf ==2/3
44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rhisf_tbl, rhisf_tbl_0   !: thickness of tbl  [m]
45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   r1_hisf_tbl              !: 1/thickness of tbl
46   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ralpha                   !: proportion of bottom cell influenced by tbl
47   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   risfLeff                 !: effective length (Leff) BG03 nn_isf==2
48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ttbl, stbl, utbl, vtbl   !: top boundary layer variable at T point
49   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf                     !: net heat flux from ice shelf      [W/m2]
50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc     !: before and now T & S isf contents [K.m/s & PSU.m/s] 
51
52   LOGICAL, PUBLIC ::   l_isfcpl = .false.       !: isf recieved from oasis
53
54   REAL(wp), PUBLIC, SAVE ::   rcpi     = 2000.0_wp     !: specific heat of ice shelf             [J/kg/K]
55   REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    !: heat diffusivity through the ice-shelf [m2/s]
56   REAL(wp), PUBLIC, SAVE ::   rhoisf   = 920.0_wp      !: volumic mass of ice shelf              [kg/m3]
57   REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      !: air temperature on top of ice shelf    [C]
58   REAL(wp), PUBLIC, SAVE ::   rlfusisf = 0.334e6_wp    !: latent heat of fusion of ice shelf     [J/kg]
59
60!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
61   CHARACTER(len=100), PUBLIC           :: cn_dirisf  = './' !: Root directory for location of ssr files
62   TYPE(FLD_N)       , PUBLIC           :: sn_fwfisf         !: information about the isf melting file to be read
63   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf
64   TYPE(FLD_N)       , PUBLIC           :: sn_rnfisf         !: information about the isf melting param.   file to be read
65   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf           
66   TYPE(FLD_N)       , PUBLIC           :: sn_depmax_isf     !: information about the grounding line depth file to be read
67   TYPE(FLD_N)       , PUBLIC           :: sn_depmin_isf     !: information about the calving   line depth file to be read
68   TYPE(FLD_N)       , PUBLIC           :: sn_Leff_isf       !: information about the effective length     file to be read
69   
70   !!----------------------------------------------------------------------
71   !! NEMO/OCE 4.0 , LOCEAN-IPSL (2017)
72   !! $Id$
73   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
74   !!----------------------------------------------------------------------
75CONTAINS
76 
77  SUBROUTINE sbc_isf( kt )
78      !!---------------------------------------------------------------------
79      !!                  ***  ROUTINE sbc_isf  ***
80      !!
81      !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf
82      !!              melting and freezing
83      !!
84      !! ** Method  :  4 parameterizations are available according to nn_isf
85      !!               nn_isf = 1 : Realistic ice_shelf formulation
86      !!                        2 : Beckmann & Goose parameterization
87      !!                        3 : Specified runoff in deptht (Mathiot & al. )
88      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf
89      !!----------------------------------------------------------------------
90      INTEGER, INTENT(in) ::   kt   ! ocean time step
91      !
92      INTEGER ::   ji, jj, jk   ! loop index
93      INTEGER ::   ikt, ikb     ! local integers
94      REAL(wp), DIMENSION(jpi,jpj) ::   zt_frz, zdep   ! freezing temperature (zt_frz) at depth (zdep)
95      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zqhcisf2d
96      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zfwfisf3d, zqhcisf3d, zqlatisf3d
97      !!---------------------------------------------------------------------
98      !
99      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN    ! compute salt and heat flux
100         !
101         SELECT CASE ( nn_isf )
102         CASE ( 1 )    ! realistic ice shelf formulation
103            ! compute T/S/U/V for the top boundary layer
104            CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')
105            CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T')
106            CALL sbc_isf_tbl(un(:,:,:)        ,utbl(:,:),'U')
107            CALL sbc_isf_tbl(vn(:,:,:)        ,vtbl(:,:),'V')
108            ! iom print
109            CALL iom_put('ttbl',ttbl(:,:))
110            CALL iom_put('stbl',stbl(:,:))
111            CALL iom_put('utbl',utbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:))
112            CALL iom_put('vtbl',vtbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:))
113            ! compute fwf and heat flux
114            ! compute fwf and heat flux
115            IF( .NOT.l_isfcpl ) THEN    ;   CALL sbc_isf_cav (kt)
116            ELSE                        ;   qisf(:,:)  = fwfisf(:,:) * rlfusisf  ! heat        flux
117            ENDIF
118            !
119         CASE ( 2 )    ! Beckmann and Goosse parametrisation
120            stbl(:,:)   = soce
121            CALL sbc_isf_bg03(kt)
122            !
123         CASE ( 3 )    ! specified runoff in depth (Mathiot et al., XXXX in preparation)
124            ! specified runoff in depth (Mathiot et al., XXXX in preparation)
125            IF( .NOT.l_isfcpl ) THEN
126               CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
127               fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting)
128            ENDIF
129            qisf(:,:)   = fwfisf(:,:) * rlfusisf             ! heat flux
130            stbl(:,:)   = soce
131            !
132         CASE ( 4 )    ! specified fwf and heat flux forcing beneath the ice shelf
133            !          ! specified fwf and heat flux forcing beneath the ice shelf
134            IF( .NOT.l_isfcpl ) THEN
135               CALL fld_read ( kt, nn_fsbc, sf_fwfisf   )
136               !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
137               fwfisf(:,:) = -sf_fwfisf(1)%fnow(:,:,1)            ! fwf
138            ENDIF
139            qisf(:,:)   = fwfisf(:,:) * rlfusisf               ! heat flux
140            stbl(:,:)   = soce
141            !
142         END SELECT
143
144         ! compute tsc due to isf
145         ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU.
146         ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rau0).
147         ! compute freezing point beneath ice shelf (or top cell if nn_isf = 3)
148         DO jj = 1,jpj
149            DO ji = 1,jpi
150               zdep(ji,jj)=gdepw_n(ji,jj,misfkt(ji,jj))
151            END DO
152         END DO
153         CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) )
154         
155         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rau0 !
156         risf_tsc(:,:,jp_sal) = 0.0_wp
157
158         ! lbclnk
159         CALL lbc_lnk_multi( risf_tsc(:,:,jp_tem), 'T', 1., risf_tsc(:,:,jp_sal), 'T', 1., fwfisf,'T', 1., qisf, 'T', 1.)
160         ! output
161         IF( iom_use('iceshelf_cea') )   CALL iom_put( 'iceshelf_cea', -fwfisf(:,:)                      )   ! isf mass flux
162         IF( iom_use('hflx_isf_cea') )   CALL iom_put( 'hflx_isf_cea', risf_tsc(:,:,jp_tem) * rau0 * rcp )   ! isf sensible+latent heat (W/m2)
163         IF( iom_use('qlatisf' ) )       CALL iom_put( 'qlatisf'     , qisf(:,:)                         )   ! isf latent heat
164         IF( iom_use('fwfisf'  ) )       CALL iom_put( 'fwfisf'      , fwfisf(:,:)                       )   ! isf mass flux (opposite sign)
165
166         ! Diagnostics
167         IF( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN
168            ALLOCATE( zfwfisf3d(jpi,jpj,jpk) , zqhcisf3d(jpi,jpj,jpk) , zqlatisf3d(jpi,jpj,jpk) )
169            ALLOCATE( zqhcisf2d(jpi,jpj) )
170            !
171            zfwfisf3d (:,:,:) = 0._wp                         ! 3d ice shelf melting (kg/m2/s)
172            zqhcisf3d (:,:,:) = 0._wp                         ! 3d heat content flux (W/m2)
173            zqlatisf3d(:,:,:) = 0._wp                         ! 3d ice shelf melting latent heat flux (W/m2)
174            zqhcisf2d (:,:)   = fwfisf(:,:) * zt_frz * rcp    ! 2d heat content flux (W/m2)
175            !
176            DO jj = 1,jpj
177               DO ji = 1,jpi
178                  ikt = misfkt(ji,jj)
179                  ikb = misfkb(ji,jj)
180                  DO jk = ikt, ikb - 1
181                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
182                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
183                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
184                  END DO
185                  zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj)   & 
186                     &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk)
187                  zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj)   & 
188                     &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk)
189                  zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj)   & 
190                     &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk)
191               END DO
192            END DO
193            !
194            CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:))
195            CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:))
196            CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:))
197            CALL iom_put('qhcisf'   , zqhcisf2d (:,:  ))
198            !
199            DEALLOCATE( zfwfisf3d, zqhcisf3d, zqlatisf3d )
200            DEALLOCATE( zqhcisf2d )
201         ENDIF
202         !
203      ENDIF
204
205      IF( kt == nit000 ) THEN                         !   set the forcing field at nit000 - 1    !
206         IF( ln_rstart .AND.    &                     ! Restart: read in restart file
207            &   iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN
208            IF(lwp) WRITE(numout,*) '          nit000-1 isf tracer content forcing fields read in the restart file'
209            CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:)         , ldxios = lrxios )   ! before salt content isf_tsc trend
210            CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b' , risf_tsc_b(:,:,jp_sal), ldxios = lrxios )   ! before salt content isf_tsc trend
211            CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b' , risf_tsc_b(:,:,jp_tem), ldxios = lrxios )   ! before salt content isf_tsc trend
212         ELSE
213            fwfisf_b(:,:)    = fwfisf(:,:)
214            risf_tsc_b(:,:,:)= risf_tsc(:,:,:)
215         ENDIF
216      ENDIF
217      !
218      IF( lrst_oce ) THEN
219         IF(lwp) WRITE(numout,*)
220         IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ',   &
221            &                    'at it= ', kt,' date= ', ndastp
222         IF(lwp) WRITE(numout,*) '~~~~'
223         IF( lwxios ) CALL iom_swap(      cwxios_context          )
224         CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:)         , ldxios = lwxios )
225         CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem), ldxios = lwxios )
226         CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal), ldxios = lwxios )
227         IF( lwxios ) CALL iom_swap(      cxios_context          )
228      ENDIF
229      !
230   END SUBROUTINE sbc_isf
231
232
233   INTEGER FUNCTION sbc_isf_alloc()
234      !!----------------------------------------------------------------------
235      !!               ***  FUNCTION sbc_isf_rnf_alloc  ***
236      !!----------------------------------------------------------------------
237      sbc_isf_alloc = 0       ! set to zero if no array to be allocated
238      IF( .NOT. ALLOCATED( qisf ) ) THEN
239         ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj)   , &
240               &    rhisf_tbl(jpi,jpj)    , r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj)  , &
241               &    ttbl(jpi,jpj)         , stbl(jpi,jpj)       , utbl(jpi,jpj)       , &
242               &    vtbl(jpi, jpj)        , risfLeff(jpi,jpj)   , rhisf_tbl_0(jpi,jpj), &
243               &    ralpha(jpi,jpj)       , misfkt(jpi,jpj)     , misfkb(jpi,jpj)     , &
244               &    STAT= sbc_isf_alloc )
245         !
246         IF( lk_mpp             )   CALL mpp_sum ( sbc_isf_alloc )
247         IF( sbc_isf_alloc /= 0 )   CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.')
248         !
249      ENDIF
250   END FUNCTION
251
252
253  SUBROUTINE sbc_isf_init
254      !!---------------------------------------------------------------------
255      !!                  ***  ROUTINE sbc_isf_init  ***
256      !!
257      !! ** Purpose : Initialisation of variables for iceshelf fluxes formulation
258      !!
259      !! ** Method  :  4 parameterizations are available according to nn_isf
260      !!               nn_isf = 1 : Realistic ice_shelf formulation
261      !!                        2 : Beckmann & Goose parameterization
262      !!                        3 : Specified runoff in deptht (Mathiot & al. )
263      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf
264      !!----------------------------------------------------------------------
265      INTEGER               :: ji, jj, jk           ! loop index
266      INTEGER               :: ik                   ! current level index
267      INTEGER               :: ikt, ikb             ! top and bottom level of the isf boundary layer
268      INTEGER               :: inum, ierror
269      INTEGER               :: ios                  ! Local integer output status for namelist read
270      REAL(wp)              :: zhk
271      CHARACTER(len=256)    :: cvarzisf, cvarhisf   ! name for isf file
272      CHARACTER(LEN=32 )    :: cvarLeff             ! variable name for efficient Length scale
273      !!----------------------------------------------------------------------
274      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf, &
275                         & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
276      !!----------------------------------------------------------------------
277
278      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs
279      READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
280901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp )
281
282      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
283      READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
284902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
285      IF(lwm) WRITE ( numond, namsbc_isf )
286
287      IF(lwp) WRITE(numout,*)
288      IF(lwp) WRITE(numout,*) 'sbc_isf_init : heat flux of the ice shelf'
289      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
290      IF(lwp) WRITE(numout,*) '   Namelist namsbc_isf :'
291      IF(lwp) WRITE(numout,*) '      type ice shelf melting/freezing         nn_isf      = ', nn_isf
292      IF(lwp) WRITE(numout,*) '      bulk formulation (nn_isf=1 only)        nn_isfblk   = ', nn_isfblk
293      IF(lwp) WRITE(numout,*) '      thickness of the top boundary layer     rn_hisf_tbl = ', rn_hisf_tbl
294      IF(lwp) WRITE(numout,*) '      gamma formulation                       nn_gammablk = ', nn_gammablk 
295      IF(lwp) WRITE(numout,*) '      gammat coefficient                      rn_gammat0  = ', rn_gammat0 
296      IF(lwp) WRITE(numout,*) '      gammas coefficient                      rn_gammas0  = ', rn_gammas0 
297      IF(lwp) WRITE(numout,*) '      top drag coef. used (from namdrg_top)   rn_Cd0      = ', r_Cdmin_top 
298
299
300                           !  1 = presence of ISF    2 = bg03 parametrisation
301                           !  3 = rnf file for isf   4 = ISF fwf specified
302                           !  option 1 and 4 need ln_isfcav = .true. (domzgr)
303      !
304      ! Allocate public variable
305      IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
306      !
307      ! initialisation
308      qisf    (:,:)    = 0._wp   ;   fwfisf  (:,:) = 0._wp
309      risf_tsc(:,:,:)  = 0._wp   ;   fwfisf_b(:,:) = 0._wp
310      !
311      ! define isf tbl tickness, top and bottom indice
312      SELECT CASE ( nn_isf )
313      CASE ( 1 ) 
314         IF(lwp) WRITE(numout,*)
315         IF(lwp) WRITE(numout,*) '      ==>>>   presence of under iceshelf seas (nn_isf = 1)'
316         rhisf_tbl(:,:) = rn_hisf_tbl
317         misfkt   (:,:) = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
318         !
319      CASE ( 2 , 3 )
320         IF( .NOT.l_isfcpl ) THEN
321             ALLOCATE( sf_rnfisf(1), STAT=ierror )
322             ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
323             CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
324          ENDIF
325          !  read effective lenght (BG03)
326          IF( nn_isf == 2 ) THEN
327            IF(lwp) WRITE(numout,*)
328            IF(lwp) WRITE(numout,*) '      ==>>>   bg03 parametrisation (nn_isf = 2)'
329            CALL iom_open( sn_Leff_isf%clname, inum )
330            cvarLeff = TRIM(sn_Leff_isf%clvar)
331            CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
332            CALL iom_close(inum)
333            !
334            risfLeff = risfLeff*1000.0_wp           !: convertion in m
335         ELSE
336            IF(lwp) WRITE(numout,*)
337            IF(lwp) WRITE(numout,*) '      ==>>>   rnf file for isf (nn_isf = 3)'
338         ENDIF
339         ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
340         CALL iom_open( sn_depmax_isf%clname, inum )
341         cvarhisf = TRIM(sn_depmax_isf%clvar)
342         CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
343         CALL iom_close(inum)
344         !
345         CALL iom_open( sn_depmin_isf%clname, inum )
346         cvarzisf = TRIM(sn_depmin_isf%clvar)
347         CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
348         CALL iom_close(inum)
349         !
350         rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer
351
352         !! compute first level of the top boundary layer
353         DO ji = 1, jpi
354            DO jj = 1, jpj
355                ik = 2
356!!gm potential bug: use gdepw_0 not _n
357                DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO
358                misfkt(ji,jj) = ik-1
359            END DO
360         END DO
361         !
362      CASE ( 4 ) 
363         IF(lwp) WRITE(numout,*)
364         IF(lwp) WRITE(numout,*) '      ==>>>   specified fresh water flux in ISF (nn_isf = 4)'
365         ! as in nn_isf == 1
366         rhisf_tbl(:,:) = rn_hisf_tbl
367         misfkt   (:,:) = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
368         !
369         ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
370         IF( .NOT.l_isfcpl ) THEN
371           ALLOCATE( sf_fwfisf(1), STAT=ierror )
372           ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
373           CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
374         ENDIF
375         !
376      CASE DEFAULT
377         CALL ctl_stop( 'sbc_isf_init: wrong value of nn_isf' )
378      END SELECT
379         
380      rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
381
382      ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
383      DO jj = 1,jpj
384         DO ji = 1,jpi
385            ikt = misfkt(ji,jj)
386            ikb = misfkt(ji,jj)
387            ! thickness of boundary layer at least the top level thickness
388            rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
389
390            ! determine the deepest level influenced by the boundary layer
391            DO jk = ikt+1, mbkt(ji,jj)
392               IF( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) )   ikb = jk
393            END DO
394            rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
395            misfkb(ji,jj) = ikb                                                   ! last wet level of the tbl
396            r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
397
398            zhk           = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1
399            ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
400         END DO
401      END DO
402
403      IF( lwxios ) THEN
404          CALL iom_set_rstw_var_active('fwf_isf_b')
405          CALL iom_set_rstw_var_active('isf_hc_b')
406          CALL iom_set_rstw_var_active('isf_sc_b')
407      ENDIF
408
409
410  END SUBROUTINE sbc_isf_init
411
412
413  SUBROUTINE sbc_isf_bg03(kt)
414      !!---------------------------------------------------------------------
415      !!                  ***  ROUTINE sbc_isf_bg03  ***
416      !!
417      !! ** Purpose : add net heat and fresh water flux from ice shelf melting
418      !!          into the adjacent ocean
419      !!
420      !! ** Method  :   See reference
421      !!
422      !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
423      !!         interaction for climate models", Ocean Modelling 5(2003) 157-170.
424      !!         (hereafter BG)
425      !! History :  06-02  (C. Wang) Original code
426      !!----------------------------------------------------------------------
427      INTEGER, INTENT ( in ) :: kt
428      !
429      INTEGER  :: ji, jj, jk ! dummy loop index
430      INTEGER  :: ik         ! current level
431      REAL(wp) :: zt_sum     ! sum of the temperature between 200m and 600m
432      REAL(wp) :: zt_ave     ! averaged temperature between 200m and 600m
433      REAL(wp) :: zt_frz     ! freezing point temperature at depth z
434      REAL(wp) :: zpress     ! pressure to compute the freezing point in depth
435      !!----------------------------------------------------------------------
436      !
437      DO ji = 1, jpi
438         DO jj = 1, jpj
439            ik = misfkt(ji,jj)
440            !! Initialize arrays to 0 (each step)
441            zt_sum = 0.e0_wp
442            IF ( ik > 1 ) THEN
443               ! 1. -----------the average temperature between 200m and 600m ---------------------
444               DO jk = misfkt(ji,jj),misfkb(ji,jj)
445                  ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!)
446                  ! after verif with UNESCO, wrong sign in BG eq. 2
447                  ! Calculate freezing temperature
448                  CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 
449                  zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp
450               END DO
451               zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
452               ! 2. ------------Net heat flux and fresh water flux due to the ice shelf
453               ! For those corresponding to zonal boundary   
454               qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
455                           & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk)
456             
457               fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf          !fresh water flux kg/(m2s)                 
458               fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
459               !add to salinity trend
460            ELSE
461               qisf(ji,jj) = 0._wp   ;   fwfisf(ji,jj) = 0._wp
462            END IF
463         END DO
464      END DO
465      !
466  END SUBROUTINE sbc_isf_bg03
467
468
469  SUBROUTINE sbc_isf_cav( kt )
470      !!---------------------------------------------------------------------
471      !!                     ***  ROUTINE sbc_isf_cav  ***
472      !!
473      !! ** Purpose :   handle surface boundary condition under ice shelf
474      !!
475      !! ** Method  : -
476      !!
477      !! ** Action  :   utau, vtau : remain unchanged
478      !!                taum, wndm : remain unchanged
479      !!                qns        : update heat flux below ice shelf
480      !!                emp, emps  : update freshwater flux below ice shelf
481      !!---------------------------------------------------------------------
482      INTEGER, INTENT(in) ::   kt   ! ocean time step
483      !
484      INTEGER  ::   ji, jj     ! dummy loop indices
485      INTEGER  ::   nit
486      LOGICAL  ::   lit
487      REAL(wp) ::   zlamb1, zlamb2, zlamb3
488      REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7
489      REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac
490      REAL(wp) ::   zeps = 1.e-20_wp       
491      REAL(wp) ::   zerr
492      REAL(wp), DIMENSION(jpi,jpj) ::   zfrz
493      REAL(wp), DIMENSION(jpi,jpj) ::   zgammat, zgammas 
494      REAL(wp), DIMENSION(jpi,jpj) ::   zfwflx, zhtflx, zhtflx_b
495      !!---------------------------------------------------------------------
496      !
497      ! coeficient for linearisation of potential tfreez
498      ! Crude approximation for pressure (but commonly used)
499      zlamb1 =-0.0573_wp
500      zlamb2 = 0.0832_wp
501      zlamb3 =-7.53e-08_wp * grav * rau0
502      !
503      ! initialisation
504      zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0
505      zhtflx (:,:) = 0.0_wp     ; zhtflx_b(:,:) = 0.0_wp   
506      zfwflx (:,:) = 0.0_wp
507
508      ! compute ice shelf melting
509      nit = 1 ; lit = .TRUE.
510      DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine
511         SELECT CASE ( nn_isfblk )
512         CASE ( 1 )   !  ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006)
513            ! Calculate freezing temperature
514            CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) )
515
516            ! compute gammat every where (2d)
517            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
518           
519            ! compute upward heat flux zhtflx and upward water flux zwflx
520            DO jj = 1, jpj
521               DO ji = 1, jpi
522                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj))
523                  zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf
524               END DO
525            END DO
526
527            ! Compute heat flux and upward fresh water flux
528            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
529            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
530
531         CASE ( 2 )  ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
532            ! compute gammat every where (2d)
533            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
534
535            ! compute upward heat flux zhtflx and upward water flux zwflx
536            ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015)
537            DO jj = 1, jpj
538               DO ji = 1, jpi
539                  ! compute coeficient to solve the 2nd order equation
540                  zeps1 = rcp*rau0*zgammat(ji,jj)
541                  zeps2 = rlfusisf*rau0*zgammas(ji,jj)
542                  zeps3 = rhoisf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps)
543                  zeps4 = zlamb2+zlamb3*risfdep(ji,jj)
544                  zeps6 = zeps4-ttbl(ji,jj)
545                  zeps7 = zeps4-tsurf
546                  zaqe  = zlamb1 * (zeps1 + zeps3)
547                  zaqer = 0.5_wp/MIN(zaqe,-zeps)
548                  zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2
549                  zcqe  = zeps2*stbl(ji,jj)
550                  zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe               
551
552                  ! Presumably zdis can never be negative because gammas is very small compared to gammat
553                  ! compute s freeze
554                  zsfrz=(-zbqe-SQRT(zdis))*zaqer
555                  IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer
556
557                  ! compute t freeze (eq. 22)
558                  zfrz(ji,jj)=zeps4+zlamb1*zsfrz
559 
560                  ! zfwflx is upward water flux
561                  ! zhtflx is upward heat flux (out of ocean)
562                  ! compute the upward water and heat flux (eq. 28 and eq. 29)
563                  zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps)
564                  zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) ) 
565               END DO
566            END DO
567
568            ! compute heat and water flux
569            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
570            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
571
572         END SELECT
573
574         ! define if we need to iterate (nn_gammablk 0/1 do not need iteration)
575         IF ( nn_gammablk <  2 ) THEN ; lit = .FALSE.
576         ELSE                           
577            ! check total number of iteration
578            IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
579            ELSE                 ; nit = nit + 1
580            END IF
581
582            ! compute error between 2 iterations
583            ! if needed save gammat and compute zhtflx_b for next iteration
584            zerr = MAXVAL(ABS(zhtflx-zhtflx_b))
585            IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE.
586            ELSE                        ; zhtflx_b(:,:) = zhtflx(:,:)
587            END IF
588         END IF
589      END DO
590      !
591      CALL iom_put('isfgammat', zgammat)
592      CALL iom_put('isfgammas', zgammas)
593      !
594   END SUBROUTINE sbc_isf_cav
595
596
597   SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf )
598      !!----------------------------------------------------------------------
599      !! ** Purpose    : compute the coefficient echange for heat flux
600      !!
601      !! ** Method     : gamma assume constant or depends of u* and stability
602      !!
603      !! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
604      !!                Jenkins et al., 2010, JPO, p2298-2312
605      !!---------------------------------------------------------------------
606      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pgt   , pgs      !
607      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pqhisf, pqwisf   !
608      !
609      INTEGER  :: ji, jj                     ! loop index
610      INTEGER  :: ikt                        ! local integer
611      REAL(wp) :: zdku, zdkv                 ! U, V shear
612      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number
613      REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
614      REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
615      REAL(wp) :: zhmax                      ! limitation of mol
616      REAL(wp) :: zetastar                   ! stability parameter
617      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence
618      REAL(wp) :: zcoef                      ! temporary coef
619      REAL(wp) :: zdep
620      REAL(wp) :: zeps = 1.0e-20_wp   
621      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant
622      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
623      REAL(wp), DIMENSION(2) :: zts, zab
624      REAL(wp), DIMENSION(jpi,jpj) :: zustar   ! U, V at T point and friction velocity
625      !!---------------------------------------------------------------------
626      !
627      SELECT CASE ( nn_gammablk )
628      CASE ( 0 ) ! gamma is constant (specified in namelist)
629         !! ISOMIP formulation (Hunter et al, 2006)
630         pgt(:,:) = rn_gammat0
631         pgs(:,:) = rn_gammas0
632
633      CASE ( 1 ) ! gamma is assume to be proportional to u*
634         !! Jenkins et al., 2010, JPO, p2298-2312
635         !! Adopted by Asay-Davis et al. (2015)
636         !! compute ustar (eq. 24)
637!!gm  NB  use pCdU here so that it will incorporate local boost of Cd0 and log layer case :
638!!         zustar(:,:) = SQRT( rCdU_top(:,:) * SQRT(utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
639!! or better :  compute ustar in zdfdrg  and use it here as well as in TKE, GLS and Co
640!!
641!!     ===>>>>    GM  to be done this chrismas
642!!         
643!!gm end
644         zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
645
646         !! Compute gammats
647         pgt(:,:) = zustar(:,:) * rn_gammat0
648         pgs(:,:) = zustar(:,:) * rn_gammas0
649     
650      CASE ( 2 ) ! gamma depends of stability of boundary layer
651         !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
652         !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
653         !! compute ustar
654         zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
655
656         !! compute Pr and Sc number (can be improved)
657         zPr =   13.8_wp
658         zSc = 2432.0_wp
659
660         !! compute gamma mole
661         zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp
662         zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp
663
664         !! compute gamma
665         DO ji = 2, jpi
666            DO jj = 2, jpj
667               ikt = mikt(ji,jj)
668
669               IF( zustar(ji,jj) == 0._wp ) THEN           ! only for kt = 1 I think
670                  pgt = rn_gammat0
671                  pgs = rn_gammas0
672               ELSE
673                  !! compute Rc number (as done in zdfric.F90)
674!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation
675!!gm moreover, use Max(rn2,0) to take care of static instabilities....
676                  zcoef = 0.5_wp / e3w_n(ji,jj,ikt)
677                  !                                            ! shear of horizontal velocity
678                  zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  &
679                     &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
680                  zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  &
681                     &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
682                  !                                            ! richardson number (minimum value set to zero)
683                  zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps )
684
685                  !! compute bouyancy
686                  zts(jp_tem) = ttbl(ji,jj)
687                  zts(jp_sal) = stbl(ji,jj)
688                  zdep        = gdepw_n(ji,jj,ikt)
689                  !
690                  CALL eos_rab( zts, zdep, zab )
691                  !
692                  !! compute length scale
693                  zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
694
695                  !! compute Monin Obukov Length
696                  ! Maximum boundary layer depth
697                  zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp
698                  ! Compute Monin obukhov length scale at the surface and Ekman depth:
699                  zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
700                  zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
701
702                  !! compute eta* (stability parameter)
703                  zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp)))
704
705                  !! compute the sublayer thickness
706                  zhnu = 5 * znu / zustar(ji,jj)
707
708                  !! compute gamma turb
709                  zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) &
710                  &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
711
712                  !! compute gammats
713                  pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
714                  pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
715               END IF
716            END DO
717         END DO
718         CALL lbc_lnk_multi( pgt, 'T', 1., pgs, 'T', 1.)
719      END SELECT
720      !
721   END SUBROUTINE sbc_isf_gammats
722
723
724   SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin )
725      !!----------------------------------------------------------------------
726      !!                  ***  SUBROUTINE sbc_isf_tbl  ***
727      !!
728      !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
729      !!
730      !!----------------------------------------------------------------------
731      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) :: pvarin
732      REAL(wp), DIMENSION(:,:)  , INTENT(  out) :: pvarout
733      CHARACTER(len=1),           INTENT(in   ) :: cd_ptin ! point of variable in/out
734      !
735      INTEGER ::   ji, jj, jk                ! loop index
736      INTEGER ::   ikt, ikb                    ! top and bottom index of the tbl
737      REAL(wp) ::   ze3, zhk
738      REAL(wp), DIMENSION(jpi,jpj) :: zhisf_tbl ! thickness of the tbl
739      !!----------------------------------------------------------------------
740     
741      ! initialisation
742      pvarout(:,:)=0._wp
743   
744      SELECT CASE ( cd_ptin )
745      CASE ( 'U' ) ! compute U in the top boundary layer at T- point
746         DO jj = 1,jpj
747            DO ji = 1,jpi
748               ikt = miku(ji,jj) ; ikb = miku(ji,jj)
749               ! thickness of boundary layer at least the top level thickness
750               zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u_n(ji,jj,ikt) )
751
752               ! determine the deepest level influenced by the boundary layer
753               DO jk = ikt+1, mbku(ji,jj)
754                  IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
755               END DO
756               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
757
758               ! level fully include in the ice shelf boundary layer
759               DO jk = ikt, ikb - 1
760                  ze3 = e3u_n(ji,jj,jk)
761                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
762               END DO
763
764               ! level partially include in ice shelf boundary layer
765               zhk = SUM( e3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
766               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
767            END DO
768         END DO
769         DO jj = 2, jpj
770            DO ji = 2, jpi
771!!gm a wet-point only average should be used here !!!
772               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj))
773            END DO
774         END DO
775         CALL lbc_lnk(pvarout,'T',-1.)
776     
777      CASE ( 'V' ) ! compute V in the top boundary layer at T- point
778         DO jj = 1,jpj
779            DO ji = 1,jpi
780               ikt = mikv(ji,jj) ; ikb = mikv(ji,jj)
781               ! thickness of boundary layer at least the top level thickness
782               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v_n(ji,jj,ikt))
783
784               ! determine the deepest level influenced by the boundary layer
785               DO jk = ikt+1, mbkv(ji,jj)
786                  IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
787               END DO
788               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
789
790               ! level fully include in the ice shelf boundary layer
791               DO jk = ikt, ikb - 1
792                  ze3 = e3v_n(ji,jj,jk)
793                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
794               END DO
795
796               ! level partially include in ice shelf boundary layer
797               zhk = SUM( e3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
798               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
799            END DO
800         END DO
801         DO jj = 2, jpj
802            DO ji = 2, jpi
803!!gm a wet-point only average should be used here !!!
804               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1))
805            END DO
806         END DO
807         CALL lbc_lnk(pvarout,'T',-1.)
808
809      CASE ( 'T' ) ! compute T in the top boundary layer at T- point
810         DO jj = 1,jpj
811            DO ji = 1,jpi
812               ikt = misfkt(ji,jj)
813               ikb = misfkb(ji,jj)
814
815               ! level fully include in the ice shelf boundary layer
816               DO jk = ikt, ikb - 1
817                  ze3 = e3t_n(ji,jj,jk)
818                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
819               END DO
820
821               ! level partially include in ice shelf boundary layer
822               zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
823               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
824            END DO
825         END DO
826      END SELECT
827      !
828      ! mask mean tbl value
829      pvarout(:,:) = pvarout(:,:) * ssmask(:,:)
830      !
831   END SUBROUTINE sbc_isf_tbl
832     
833
834   SUBROUTINE sbc_isf_div( phdivn )
835      !!----------------------------------------------------------------------
836      !!                  ***  SUBROUTINE sbc_isf_div  ***
837      !!       
838      !! ** Purpose :   update the horizontal divergence with the runoff inflow
839      !!
840      !! ** Method  :   
841      !!                CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the
842      !!                          divergence and expressed in m/s
843      !!
844      !! ** Action  :   phdivn   decreased by the runoff inflow
845      !!----------------------------------------------------------------------
846      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdivn   ! horizontal divergence
847      !
848      INTEGER  ::   ji, jj, jk   ! dummy loop indices
849      INTEGER  ::   ikt, ikb 
850      REAL(wp) ::   zhk
851      REAL(wp) ::   zfact     ! local scalar
852      !!----------------------------------------------------------------------
853      !
854      zfact   = 0.5_wp
855      !
856      IF(.NOT.ln_linssh ) THEN     ! need to re compute level distribution of isf fresh water
857         DO jj = 1,jpj
858            DO ji = 1,jpi
859               ikt = misfkt(ji,jj)
860               ikb = misfkt(ji,jj)
861               ! thickness of boundary layer at least the top level thickness
862               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
863
864               ! determine the deepest level influenced by the boundary layer
865               DO jk = ikt, mbkt(ji,jj)
866                  IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
867               END DO
868               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
869               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
870               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
871
872               zhk           = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
873               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
874            END DO
875         END DO
876      END IF 
877      !
878      !==   ice shelf melting distributed over several levels   ==!
879      DO jj = 1,jpj
880         DO ji = 1,jpi
881               ikt = misfkt(ji,jj)
882               ikb = misfkb(ji,jj)
883               ! level fully include in the ice shelf boundary layer
884               DO jk = ikt, ikb - 1
885                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
886                    &              * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
887               END DO
888               ! level partially include in ice shelf boundary layer
889               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
890                    &            + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
891         END DO
892      END DO
893      !
894   END SUBROUTINE sbc_isf_div
895
896   !!======================================================================
897END MODULE sbcisf
Note: See TracBrowser for help on using the repository browser.