source: NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/SBC/sbcisf.F90 @ 9923

Last change on this file since 9923 was 9923, checked in by gm, 3 years ago

#1911 (ENHANCE-04): step I.2: dev_r9838_ENHANCE04_MLF

  • Property svn:keywords set to Id
File size: 45.4 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 ::   rho_isf  = 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 , NEMO Consortium (2018)
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 / rho0).
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_rho0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rho0 !
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) * rho0 * 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) = - rho0 * 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      IF ( l_useCT ) THEN   ! linearisation from Jourdain et al. (2017)
500         zlamb1 =-0.0564_wp
501         zlamb2 = 0.0773_wp
502         zlamb3 =-7.8633e-8 * grav * rho0
503      ELSE                  ! linearisation from table 4 (Asay-Davis et al., 2015)
504         zlamb1 =-0.0573_wp
505         zlamb2 = 0.0832_wp
506         zlamb3 =-7.53e-8 * grav * rho0
507      ENDIF
508      !
509      ! initialisation
510      zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0
511      zhtflx (:,:) = 0.0_wp     ; zhtflx_b(:,:) = 0.0_wp   
512      zfwflx (:,:) = 0.0_wp
513
514      ! compute ice shelf melting
515      nit = 1 ; lit = .TRUE.
516      DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine
517         SELECT CASE ( nn_isfblk )
518         CASE ( 1 )   !  ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006)
519            ! Calculate freezing temperature
520            CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) )
521
522            ! compute gammat every where (2d)
523            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
524           
525            ! compute upward heat flux zhtflx and upward water flux zwflx
526            DO jj = 1, jpj
527               DO ji = 1, jpi
528                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rho0*(ttbl(ji,jj)-zfrz(ji,jj))
529                  zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf
530               END DO
531            END DO
532
533            ! Compute heat flux and upward fresh water flux
534            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
535            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
536
537         CASE ( 2 )  ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
538            ! compute gammat every where (2d)
539            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
540
541            ! compute upward heat flux zhtflx and upward water flux zwflx
542            ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015)
543            DO jj = 1, jpj
544               DO ji = 1, jpi
545                  ! compute coeficient to solve the 2nd order equation
546                  zeps1 = rcp*rho0*zgammat(ji,jj)
547                  zeps2 = rlfusisf*rho0*zgammas(ji,jj)
548                  zeps3 = rho_isf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps)
549                  zeps4 = zlamb2+zlamb3*risfdep(ji,jj)
550                  zeps6 = zeps4-ttbl(ji,jj)
551                  zeps7 = zeps4-tsurf
552                  zaqe  = zlamb1 * (zeps1 + zeps3)
553                  zaqer = 0.5_wp/MIN(zaqe,-zeps)
554                  zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2
555                  zcqe  = zeps2*stbl(ji,jj)
556                  zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe               
557
558                  ! Presumably zdis can never be negative because gammas is very small compared to gammat
559                  ! compute s freeze
560                  zsfrz=(-zbqe-SQRT(zdis))*zaqer
561                  IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer
562
563                  ! compute t freeze (eq. 22)
564                  zfrz(ji,jj)=zeps4+zlamb1*zsfrz
565 
566                  ! zfwflx is upward water flux
567                  ! zhtflx is upward heat flux (out of ocean)
568                  ! compute the upward water and heat flux (eq. 28 and eq. 29)
569                  zfwflx(ji,jj) = rho0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps)
570                  zhtflx(ji,jj) = zgammat(ji,jj) * rho0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) ) 
571               END DO
572            END DO
573
574            ! compute heat and water flux
575            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
576            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
577
578         END SELECT
579
580         ! define if we need to iterate (nn_gammablk 0/1 do not need iteration)
581         IF ( nn_gammablk <  2 ) THEN ; lit = .FALSE.
582         ELSE                           
583            ! check total number of iteration
584            IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
585            ELSE                 ; nit = nit + 1
586            END IF
587
588            ! compute error between 2 iterations
589            ! if needed save gammat and compute zhtflx_b for next iteration
590            zerr = MAXVAL(ABS(zhtflx-zhtflx_b))
591            IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE.
592            ELSE                        ; zhtflx_b(:,:) = zhtflx(:,:)
593            END IF
594         END IF
595      END DO
596      !
597      CALL iom_put('isfgammat', zgammat)
598      CALL iom_put('isfgammas', zgammas)
599      !
600   END SUBROUTINE sbc_isf_cav
601
602
603   SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf )
604      !!----------------------------------------------------------------------
605      !! ** Purpose    : compute the coefficient echange for heat flux
606      !!
607      !! ** Method     : gamma assume constant or depends of u* and stability
608      !!
609      !! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
610      !!                Jenkins et al., 2010, JPO, p2298-2312
611      !!---------------------------------------------------------------------
612      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pgt   , pgs      !
613      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pqhisf, pqwisf   !
614      !
615      INTEGER  :: ji, jj                     ! loop index
616      INTEGER  :: ikt                        ! local integer
617      REAL(wp) :: zdku, zdkv                 ! U, V shear
618      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number
619      REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
620      REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
621      REAL(wp) :: zhmax                      ! limitation of mol
622      REAL(wp) :: zetastar                   ! stability parameter
623      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence
624      REAL(wp) :: zcoef                      ! temporary coef
625      REAL(wp) :: zdep
626      REAL(wp) :: zeps = 1.0e-20_wp   
627      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant
628      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
629      REAL(wp), DIMENSION(2) :: zts, zab
630      REAL(wp), DIMENSION(jpi,jpj) :: zustar   ! U, V at T point and friction velocity
631      !!---------------------------------------------------------------------
632      !
633      SELECT CASE ( nn_gammablk )
634      CASE ( 0 ) ! gamma is constant (specified in namelist)
635         !! ISOMIP formulation (Hunter et al, 2006)
636         pgt(:,:) = rn_gammat0
637         pgs(:,:) = rn_gammas0
638
639      CASE ( 1 ) ! gamma is assume to be proportional to u*
640         !! Jenkins et al., 2010, JPO, p2298-2312
641         !! Adopted by Asay-Davis et al. (2015)
642         !! compute ustar (eq. 24)
643!!gm  NB  use pCdU here so that it will incorporate local boost of Cd0 and log layer case :
644!!         zustar(:,:) = SQRT( rCdU_top(:,:) * SQRT(utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
645!! or better :  compute ustar in zdfdrg  and use it here as well as in TKE, GLS and Co
646!!
647!!     ===>>>>    GM  to be done this chrismas
648!!         
649!!gm end
650         zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
651
652         !! Compute gammats
653         pgt(:,:) = zustar(:,:) * rn_gammat0
654         pgs(:,:) = zustar(:,:) * rn_gammas0
655     
656      CASE ( 2 ) ! gamma depends of stability of boundary layer
657         !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
658         !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
659         !! compute ustar
660         zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
661
662         !! compute Pr and Sc number (can be improved)
663         zPr =   13.8_wp
664         zSc = 2432.0_wp
665
666         !! compute gamma mole
667         zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp
668         zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp
669
670         !! compute gamma
671         DO ji = 2, jpi
672            DO jj = 2, jpj
673               ikt = mikt(ji,jj)
674
675               IF( zustar(ji,jj) == 0._wp ) THEN           ! only for kt = 1 I think
676                  pgt = rn_gammat0
677                  pgs = rn_gammas0
678               ELSE
679                  !! compute Rc number (as done in zdfric.F90)
680!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation
681!!gm moreover, use Max(rn2,0) to take care of static instabilities....
682                  zcoef = 0.5_wp / e3w_n(ji,jj,ikt)
683                  !                                            ! shear of horizontal velocity
684                  zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  &
685                     &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
686                  zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  &
687                     &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
688                  !                                            ! richardson number (minimum value set to zero)
689                  zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps )
690
691                  !! compute bouyancy
692                  zts(jp_tem) = ttbl(ji,jj)
693                  zts(jp_sal) = stbl(ji,jj)
694                  zdep        = gdepw_n(ji,jj,ikt)
695                  !
696                  CALL eos_rab( zts, zdep, zab )
697                  !
698                  !! compute length scale
699                  zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
700
701                  !! compute Monin Obukov Length
702                  ! Maximum boundary layer depth
703                  zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp
704                  ! Compute Monin obukhov length scale at the surface and Ekman depth:
705                  zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
706                  zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
707
708                  !! compute eta* (stability parameter)
709                  zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp)))
710
711                  !! compute the sublayer thickness
712                  zhnu = 5 * znu / zustar(ji,jj)
713
714                  !! compute gamma turb
715                  zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) &
716                  &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
717
718                  !! compute gammats
719                  pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
720                  pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
721               END IF
722            END DO
723         END DO
724         CALL lbc_lnk_multi( pgt, 'T', 1., pgs, 'T', 1.)
725      END SELECT
726      !
727   END SUBROUTINE sbc_isf_gammats
728
729
730   SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin )
731      !!----------------------------------------------------------------------
732      !!                  ***  SUBROUTINE sbc_isf_tbl  ***
733      !!
734      !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
735      !!
736      !!----------------------------------------------------------------------
737      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) :: pvarin
738      REAL(wp), DIMENSION(:,:)  , INTENT(  out) :: pvarout
739      CHARACTER(len=1),           INTENT(in   ) :: cd_ptin ! point of variable in/out
740      !
741      INTEGER ::   ji, jj, jk                ! loop index
742      INTEGER ::   ikt, ikb                    ! top and bottom index of the tbl
743      REAL(wp) ::   ze3, zhk
744      REAL(wp), DIMENSION(jpi,jpj) :: zhisf_tbl ! thickness of the tbl
745      !!----------------------------------------------------------------------
746     
747      ! initialisation
748      pvarout(:,:)=0._wp
749   
750      SELECT CASE ( cd_ptin )
751      CASE ( 'U' ) ! compute U in the top boundary layer at T- point
752         DO jj = 1,jpj
753            DO ji = 1,jpi
754               ikt = miku(ji,jj) ; ikb = miku(ji,jj)
755               ! thickness of boundary layer at least the top level thickness
756               zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u_n(ji,jj,ikt) )
757
758               ! determine the deepest level influenced by the boundary layer
759               DO jk = ikt+1, mbku(ji,jj)
760                  IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
761               END DO
762               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
763
764               ! level fully include in the ice shelf boundary layer
765               DO jk = ikt, ikb - 1
766                  ze3 = e3u_n(ji,jj,jk)
767                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
768               END DO
769
770               ! level partially include in ice shelf boundary layer
771               zhk = SUM( e3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
772               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
773            END DO
774         END DO
775         DO jj = 2, jpj
776            DO ji = 2, jpi
777!!gm a wet-point only average should be used here !!!
778               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj))
779            END DO
780         END DO
781         CALL lbc_lnk(pvarout,'T',-1.)
782     
783      CASE ( 'V' ) ! compute V in the top boundary layer at T- point
784         DO jj = 1,jpj
785            DO ji = 1,jpi
786               ikt = mikv(ji,jj) ; ikb = mikv(ji,jj)
787               ! thickness of boundary layer at least the top level thickness
788               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v_n(ji,jj,ikt))
789
790               ! determine the deepest level influenced by the boundary layer
791               DO jk = ikt+1, mbkv(ji,jj)
792                  IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
793               END DO
794               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
795
796               ! level fully include in the ice shelf boundary layer
797               DO jk = ikt, ikb - 1
798                  ze3 = e3v_n(ji,jj,jk)
799                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
800               END DO
801
802               ! level partially include in ice shelf boundary layer
803               zhk = SUM( e3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
804               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
805            END DO
806         END DO
807         DO jj = 2, jpj
808            DO ji = 2, jpi
809!!gm a wet-point only average should be used here !!!
810               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1))
811            END DO
812         END DO
813         CALL lbc_lnk(pvarout,'T',-1.)
814
815      CASE ( 'T' ) ! compute T in the top boundary layer at T- point
816         DO jj = 1,jpj
817            DO ji = 1,jpi
818               ikt = misfkt(ji,jj)
819               ikb = misfkb(ji,jj)
820
821               ! level fully include in the ice shelf boundary layer
822               DO jk = ikt, ikb - 1
823                  ze3 = e3t_n(ji,jj,jk)
824                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
825               END DO
826
827               ! level partially include in ice shelf boundary layer
828               zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
829               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
830            END DO
831         END DO
832      END SELECT
833      !
834      ! mask mean tbl value
835      pvarout(:,:) = pvarout(:,:) * ssmask(:,:)
836      !
837   END SUBROUTINE sbc_isf_tbl
838     
839
840   SUBROUTINE sbc_isf_div( phdivn )
841      !!----------------------------------------------------------------------
842      !!                  ***  SUBROUTINE sbc_isf_div  ***
843      !!       
844      !! ** Purpose :   update the horizontal divergence with the runoff inflow
845      !!
846      !! ** Method  :   
847      !!                CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the
848      !!                          divergence and expressed in m/s
849      !!
850      !! ** Action  :   phdivn   decreased by the runoff inflow
851      !!----------------------------------------------------------------------
852      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdivn   ! horizontal divergence
853      !
854      INTEGER  ::   ji, jj, jk   ! dummy loop indices
855      INTEGER  ::   ikt, ikb 
856      REAL(wp) ::   zhk
857      REAL(wp) ::   zfact     ! local scalar
858      !!----------------------------------------------------------------------
859      !
860      zfact   = 0.5_wp
861      !
862      IF(.NOT.ln_linssh ) THEN     ! need to re compute level distribution of isf fresh water
863         DO jj = 1,jpj
864            DO ji = 1,jpi
865               ikt = misfkt(ji,jj)
866               ikb = misfkt(ji,jj)
867               ! thickness of boundary layer at least the top level thickness
868               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
869
870               ! determine the deepest level influenced by the boundary layer
871               DO jk = ikt, mbkt(ji,jj)
872                  IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
873               END DO
874               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
875               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
876               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
877
878               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
879               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
880            END DO
881         END DO
882      END IF 
883      !
884      !==   ice shelf melting distributed over several levels   ==!
885      DO jj = 1,jpj
886         DO ji = 1,jpi
887               ikt = misfkt(ji,jj)
888               ikb = misfkb(ji,jj)
889               ! level fully include in the ice shelf boundary layer
890               DO jk = ikt, ikb - 1
891                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
892                    &              * r1_hisf_tbl(ji,jj) * r1_rho0 * zfact
893               END DO
894               ! level partially include in ice shelf boundary layer
895               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
896                    &            + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rho0 * zfact * ralpha(ji,jj) 
897         END DO
898      END DO
899      !
900   END SUBROUTINE sbc_isf_div
901
902   !!======================================================================
903END MODULE sbcisf
Note: See TracBrowser for help on using the repository browser.