New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcisf.F90 in NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/SBC/sbcisf.F90 @ 15231

Last change on this file since 15231 was 15231, checked in by clem, 2 years ago

4.0-HEAD: small bug fixes following tickets #2644 #2679 #2688

  • 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 ::   rcpisf   = 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 , NEMO Consortium (2018)
72   !! $Id$
73   !! Software governed by the CeCILL license (see ./LICENSE)
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( 'sbcisf', 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         CALL mpp_sum ( 'sbcisf', sbc_isf_alloc )
247         IF( sbc_isf_alloc /= 0 )   CALL ctl_stop( 'STOP', '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' )
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' )
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                DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_0(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO
357                misfkt(ji,jj) = ik-1
358            END DO
359         END DO
360         !
361      CASE ( 4 ) 
362         IF(lwp) WRITE(numout,*)
363         IF(lwp) WRITE(numout,*) '      ==>>>   specified fresh water flux in ISF (nn_isf = 4)'
364         ! as in nn_isf == 1
365         rhisf_tbl(:,:) = rn_hisf_tbl
366         misfkt   (:,:) = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
367         !
368         ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
369         IF( .NOT.l_isfcpl ) THEN
370           ALLOCATE( sf_fwfisf(1), STAT=ierror )
371           ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
372           CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
373         ENDIF
374         !
375      CASE DEFAULT
376         CALL ctl_stop( 'sbc_isf_init: wrong value of nn_isf' )
377      END SELECT
378         
379      rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
380
381      ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
382      DO jj = 1,jpj
383         DO ji = 1,jpi
384            ikt = misfkt(ji,jj)
385            ikb = misfkt(ji,jj)
386            ! thickness of boundary layer at least the top level thickness
387            rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
388
389            ! determine the deepest level influenced by the boundary layer
390            DO jk = ikt+1, mbkt(ji,jj)
391               IF( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) )   ikb = jk
392            END DO
393            rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
394            misfkb(ji,jj) = ikb                                                   ! last wet level of the tbl
395            r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
396
397            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
398            ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
399         END DO
400      END DO
401
402      IF( lwxios ) THEN
403          CALL iom_set_rstw_var_active('fwf_isf_b')
404          CALL iom_set_rstw_var_active('isf_hc_b')
405          CALL iom_set_rstw_var_active('isf_sc_b')
406      ENDIF
407
408
409  END SUBROUTINE sbc_isf_init
410
411
412  SUBROUTINE sbc_isf_bg03(kt)
413      !!---------------------------------------------------------------------
414      !!                  ***  ROUTINE sbc_isf_bg03  ***
415      !!
416      !! ** Purpose : add net heat and fresh water flux from ice shelf melting
417      !!          into the adjacent ocean
418      !!
419      !! ** Method  :   See reference
420      !!
421      !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
422      !!         interaction for climate models", Ocean Modelling 5(2003) 157-170.
423      !!         (hereafter BG)
424      !! History :  06-02  (C. Wang) Original code
425      !!----------------------------------------------------------------------
426      INTEGER, INTENT ( in ) :: kt
427      !
428      INTEGER  :: ji, jj, jk ! dummy loop index
429      INTEGER  :: ik         ! current level
430      REAL(wp) :: zt_sum     ! sum of the temperature between 200m and 600m
431      REAL(wp) :: zt_ave     ! averaged temperature between 200m and 600m
432      REAL(wp) :: zt_frz     ! freezing point temperature at depth z
433      REAL(wp) :: zpress     ! pressure to compute the freezing point in depth
434      !!----------------------------------------------------------------------
435      !
436      DO ji = 1, jpi
437         DO jj = 1, jpj
438            ik = misfkt(ji,jj)
439            !! Initialize arrays to 0 (each step)
440            zt_sum = 0.e0_wp
441            IF ( ik > 1 ) THEN
442               ! 1. -----------the average temperature between 200m and 600m ---------------------
443               DO jk = misfkt(ji,jj),misfkb(ji,jj)
444                  ! Calculate freezing temperature
445                  zpress = grav*rau0*gdept_n(ji,jj,ik)*1.e-04
446                  CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 
447                  zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp
448               END DO
449               zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
450               ! 2. ------------Net heat flux and fresh water flux due to the ice shelf
451               ! For those corresponding to zonal boundary   
452               qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
453                           & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk)
454             
455               fwfisf(ji,jj) = qisf(ji,jj) / rLfusisf          !fresh water flux kg/(m2s)                 
456               fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
457               !add to salinity trend
458            ELSE
459               qisf(ji,jj) = 0._wp   ;   fwfisf(ji,jj) = 0._wp
460            END IF
461         END DO
462      END DO
463      !
464  END SUBROUTINE sbc_isf_bg03
465
466
467  SUBROUTINE sbc_isf_cav( kt )
468      !!---------------------------------------------------------------------
469      !!                     ***  ROUTINE sbc_isf_cav  ***
470      !!
471      !! ** Purpose :   handle surface boundary condition under ice shelf
472      !!
473      !! ** Method  : -
474      !!
475      !! ** Action  :   utau, vtau : remain unchanged
476      !!                taum, wndm : remain unchanged
477      !!                qns        : update heat flux below ice shelf
478      !!                emp, emps  : update freshwater flux below ice shelf
479      !!---------------------------------------------------------------------
480      INTEGER, INTENT(in) ::   kt   ! ocean time step
481      !
482      INTEGER  ::   ji, jj     ! dummy loop indices
483      INTEGER  ::   nit
484      LOGICAL  ::   lit
485      REAL(wp) ::   zlamb1, zlamb2, zlamb3
486      REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7
487      REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac
488      REAL(wp) ::   zeps = 1.e-20_wp       
489      REAL(wp) ::   zerr
490      REAL(wp), DIMENSION(jpi,jpj) ::   zfrz
491      REAL(wp), DIMENSION(jpi,jpj) ::   zgammat, zgammas 
492      REAL(wp), DIMENSION(jpi,jpj) ::   zfwflx, zhtflx, zhtflx_b
493      !!---------------------------------------------------------------------
494      !
495      ! coeficient for linearisation of potential tfreez
496      ! Crude approximation for pressure (but commonly used)
497      IF ( l_useCT ) THEN   ! linearisation from Jourdain et al. (2017)
498         zlamb1 =-0.0564_wp
499         zlamb2 = 0.0773_wp
500         zlamb3 =-7.8633e-8 * grav * rau0
501      ELSE                  ! linearisation from table 4 (Asay-Davis et al., 2015)
502         zlamb1 =-0.0573_wp
503         zlamb2 = 0.0832_wp
504         zlamb3 =-7.53e-8 * grav * rau0
505      ENDIF
506      !
507      ! initialisation
508      zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0
509      zhtflx (:,:) = 0.0_wp     ; zhtflx_b(:,:) = 0.0_wp   
510      zfwflx (:,:) = 0.0_wp
511
512      ! compute ice shelf melting
513      nit = 1 ; lit = .TRUE.
514      DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine
515         SELECT CASE ( nn_isfblk )
516         CASE ( 1 )   !  ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006)
517            ! Calculate freezing temperature
518            CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) )
519
520            ! compute gammat every where (2d)
521            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
522           
523            ! compute upward heat flux zhtflx and upward water flux zwflx
524            DO jj = 1, jpj
525               DO ji = 1, jpi
526                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj))
527                  zfwflx(ji,jj) = - zhtflx(ji,jj)/rLfusisf
528               END DO
529            END DO
530
531            ! Compute heat flux and upward fresh water flux
532            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
533            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
534
535         CASE ( 2 )  ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
536            ! compute gammat every where (2d)
537            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
538
539            ! compute upward heat flux zhtflx and upward water flux zwflx
540            ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015)
541            DO jj = 1, jpj
542               DO ji = 1, jpi
543                  ! compute coeficient to solve the 2nd order equation
544                  zeps1 = rcp*rau0*zgammat(ji,jj)
545                  zeps2 = rLfusisf*rau0*zgammas(ji,jj)
546                  zeps3 = rhoisf*rcpisf*rkappa/MAX(risfdep(ji,jj),zeps)
547                  zeps4 = zlamb2+zlamb3*risfdep(ji,jj)
548                  zeps6 = zeps4-ttbl(ji,jj)
549                  zeps7 = zeps4-tsurf
550                  zaqe  = zlamb1 * (zeps1 + zeps3)
551                  zaqer = 0.5_wp/MIN(zaqe,-zeps)
552                  zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2
553                  zcqe  = zeps2*stbl(ji,jj)
554                  zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe               
555
556                  ! Presumably zdis can never be negative because gammas is very small compared to gammat
557                  ! compute s freeze
558                  zsfrz=(-zbqe-SQRT(zdis))*zaqer
559                  IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer
560
561                  ! compute t freeze (eq. 22)
562                  zfrz(ji,jj)=zeps4+zlamb1*zsfrz
563 
564                  ! zfwflx is upward water flux
565                  ! zhtflx is upward heat flux (out of ocean)
566                  ! compute the upward water and heat flux (eq. 28 and eq. 29)
567                  zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps)
568                  zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) ) 
569               END DO
570            END DO
571
572            ! compute heat and water flux
573            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
574            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
575
576         END SELECT
577
578         ! define if we need to iterate (nn_gammablk 0/1 do not need iteration)
579         IF ( nn_gammablk <  2 ) THEN ; lit = .FALSE.
580         ELSE                           
581            ! check total number of iteration
582            IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
583            ELSE                 ; nit = nit + 1
584            END IF
585
586            ! compute error between 2 iterations
587            ! if needed save gammat and compute zhtflx_b for next iteration
588            zerr = MAXVAL(ABS(zhtflx-zhtflx_b))
589            IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE.
590            ELSE                        ; zhtflx_b(:,:) = zhtflx(:,:)
591            END IF
592         END IF
593      END DO
594      !
595      CALL iom_put('isfgammat', zgammat)
596      CALL iom_put('isfgammas', zgammas)
597      !
598   END SUBROUTINE sbc_isf_cav
599
600
601   SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf )
602      !!----------------------------------------------------------------------
603      !! ** Purpose    : compute the coefficient echange for heat flux
604      !!
605      !! ** Method     : gamma assume constant or depends of u* and stability
606      !!
607      !! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
608      !!                Jenkins et al., 2010, JPO, p2298-2312
609      !!---------------------------------------------------------------------
610      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pgt   , pgs      !
611      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pqhisf, pqwisf   !
612      !
613      INTEGER  :: ji, jj                     ! loop index
614      INTEGER  :: ikt                        ! local integer
615      REAL(wp) :: zdku, zdkv                 ! U, V shear
616      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number
617      REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
618      REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
619      REAL(wp) :: zhmax                      ! limitation of mol
620      REAL(wp) :: zetastar                   ! stability parameter
621      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence
622      REAL(wp) :: zcoef                      ! temporary coef
623      REAL(wp) :: zdep
624      REAL(wp) :: zeps = 1.0e-20_wp   
625      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant
626      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
627      REAL(wp), DIMENSION(2) :: zts, zab
628      REAL(wp), DIMENSION(jpi,jpj) :: zustar   ! U, V at T point and friction velocity
629      !!---------------------------------------------------------------------
630      !
631      SELECT CASE ( nn_gammablk )
632      CASE ( 0 ) ! gamma is constant (specified in namelist)
633         !! ISOMIP formulation (Hunter et al, 2006)
634         pgt(:,:) = rn_gammat0
635         pgs(:,:) = rn_gammas0
636
637      CASE ( 1 ) ! gamma is assume to be proportional to u*
638         !! Jenkins et al., 2010, JPO, p2298-2312
639         !! Adopted by Asay-Davis et al. (2015)
640         !! compute ustar (eq. 24)
641!!gm  NB  use pCdU here so that it will incorporate local boost of Cd0 and log layer case :
642!!         zustar(:,:) = SQRT( rCdU_top(:,:) * SQRT(utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
643!! or better :  compute ustar in zdfdrg  and use it here as well as in TKE, GLS and Co
644!!
645!!     ===>>>>    GM  to be done this chrismas
646!!         
647!!gm end
648         zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
649
650         !! Compute gammats
651         pgt(:,:) = zustar(:,:) * rn_gammat0
652         pgs(:,:) = zustar(:,:) * rn_gammas0
653     
654      CASE ( 2 ) ! gamma depends of stability of boundary layer
655         !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
656         !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
657         !! compute ustar
658         zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
659
660         !! compute Pr and Sc number (can be improved)
661         zPr =   13.8_wp
662         zSc = 2432.0_wp
663
664         !! compute gamma mole
665         zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp
666         zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp
667
668         !! compute gamma
669         DO ji = 2, jpi
670            DO jj = 2, jpj
671               ikt = mikt(ji,jj)
672
673               IF( zustar(ji,jj) == 0._wp ) THEN           ! only for kt = 1 I think
674                  pgt = rn_gammat0
675                  pgs = rn_gammas0
676               ELSE
677                  !! compute Rc number (as done in zdfric.F90)
678!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation
679!!gm moreover, use Max(rn2,0) to take care of static instabilities....
680                  zcoef = 0.5_wp / e3w_n(ji,jj,ikt+1)
681                  !                                            ! shear of horizontal velocity
682                  zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  &
683                     &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
684                  zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  &
685                     &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
686                  !                                            ! richardson number (minimum value set to zero)
687                  zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps )
688
689                  !! compute bouyancy
690                  zts(jp_tem) = ttbl(ji,jj)
691                  zts(jp_sal) = stbl(ji,jj)
692                  zdep        = gdepw_n(ji,jj,ikt)
693                  !
694                  CALL eos_rab( zts, zdep, zab )
695                  !
696                  !! compute length scale
697                  zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
698
699                  !! compute Monin Obukov Length
700                  ! Maximum boundary layer depth
701                  zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp
702                  ! Compute Monin obukhov length scale at the surface and Ekman depth:
703                  zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
704                  zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
705
706                  !! compute eta* (stability parameter)
707                  zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp)))
708
709                  !! compute the sublayer thickness
710                  zhnu = 5 * znu / zustar(ji,jj)
711
712                  !! compute gamma turb
713                  zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) &
714                  &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
715
716                  !! compute gammats
717                  pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
718                  pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
719               END IF
720            END DO
721         END DO
722         CALL lbc_lnk_multi( 'sbcisf', pgt, 'T', 1., pgs, 'T', 1.)
723      END SELECT
724      !
725   END SUBROUTINE sbc_isf_gammats
726
727
728   SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin )
729      !!----------------------------------------------------------------------
730      !!                  ***  SUBROUTINE sbc_isf_tbl  ***
731      !!
732      !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
733      !!
734      !!----------------------------------------------------------------------
735      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) :: pvarin
736      REAL(wp), DIMENSION(:,:)  , INTENT(  out) :: pvarout
737      CHARACTER(len=1),           INTENT(in   ) :: cd_ptin ! point of variable in/out
738      !
739      INTEGER ::   ji, jj, jk                   ! loop index
740      INTEGER ::   ikt, ikb                     ! top and bottom index of the tbl
741      REAL(wp) ::   ze3, zhk
742      REAL(wp), DIMENSION(jpi,jpj) :: zhisf_tbl ! thickness of the tbl
743      REAL(wp), DIMENSION(jpi,jpj) :: zvarout
744      !!----------------------------------------------------------------------
745     
746      ! initialisation
747      pvarout(:,:)=0._wp
748   
749      SELECT CASE ( cd_ptin )
750      CASE ( 'U' ) ! compute U in the top boundary layer at T- point
751         !
752         zvarout(:,:)=0._wp
753         !
754         DO jj = 1,jpj
755            DO ji = 1,jpi
756               ikt = miku(ji,jj) ; ikb = miku(ji,jj)
757               ! thickness of boundary layer at least the top level thickness
758               zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u_n(ji,jj,ikt) )
759
760               ! determine the deepest level influenced by the boundary layer
761               DO jk = ikt+1, mbku(ji,jj)
762                  IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
763               END DO
764               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
765
766               ! level fully include in the ice shelf boundary layer
767               DO jk = ikt, ikb - 1
768                  ze3 = e3u_n(ji,jj,jk)
769                  zvarout(ji,jj) = zvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
770               END DO
771
772               ! level partially include in ice shelf boundary layer
773               zhk = SUM( e3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
774               zvarout(ji,jj) = zvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
775            END DO
776         END DO
777         DO jj = 2, jpj
778            DO ji = 2, jpi
779!!gm a wet-point only average should be used here !!!
780               pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji-1,jj))
781            END DO
782         END DO
783         CALL lbc_lnk('sbcisf', pvarout,'T',-1.)
784     
785      CASE ( 'V' ) ! compute V in the top boundary layer at T- point
786         !
787         zvarout(:,:)=0._wp
788         !
789         DO jj = 1,jpj
790            DO ji = 1,jpi
791               ikt = mikv(ji,jj) ; ikb = mikv(ji,jj)
792               ! thickness of boundary layer at least the top level thickness
793               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v_n(ji,jj,ikt))
794
795               ! determine the deepest level influenced by the boundary layer
796               DO jk = ikt+1, mbkv(ji,jj)
797                  IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
798               END DO
799               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
800
801               ! level fully include in the ice shelf boundary layer
802               DO jk = ikt, ikb - 1
803                  ze3 = e3v_n(ji,jj,jk)
804                  zvarout(ji,jj) = zvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
805               END DO
806
807               ! level partially include in ice shelf boundary layer
808               zhk = SUM( e3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
809               zvarout(ji,jj) = zvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
810            END DO
811         END DO
812         DO jj = 2, jpj
813            DO ji = 2, jpi
814!!gm a wet-point only average should be used here !!!
815               pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji,jj-1))
816            END DO
817         END DO
818         CALL lbc_lnk('sbcisf', pvarout,'T',-1.)
819
820      CASE ( 'T' ) ! compute T in the top boundary layer at T- point
821         DO jj = 1,jpj
822            DO ji = 1,jpi
823               ikt = misfkt(ji,jj)
824               ikb = misfkb(ji,jj)
825
826               ! level fully include in the ice shelf boundary layer
827               DO jk = ikt, ikb - 1
828                  ze3 = e3t_n(ji,jj,jk)
829                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
830               END DO
831
832               ! level partially include in ice shelf boundary layer
833               zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
834               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
835            END DO
836         END DO
837      END SELECT
838      !
839      ! mask mean tbl value
840      pvarout(:,:) = pvarout(:,:) * ssmask(:,:)
841      !
842   END SUBROUTINE sbc_isf_tbl
843     
844
845   SUBROUTINE sbc_isf_div( phdivn )
846      !!----------------------------------------------------------------------
847      !!                  ***  SUBROUTINE sbc_isf_div  ***
848      !!       
849      !! ** Purpose :   update the horizontal divergence with the runoff inflow
850      !!
851      !! ** Method  :   
852      !!                CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the
853      !!                          divergence and expressed in m/s
854      !!
855      !! ** Action  :   phdivn   decreased by the runoff inflow
856      !!----------------------------------------------------------------------
857      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdivn   ! horizontal divergence
858      !
859      INTEGER  ::   ji, jj, jk   ! dummy loop indices
860      INTEGER  ::   ikt, ikb 
861      REAL(wp) ::   zhk
862      REAL(wp) ::   zfact     ! local scalar
863      !!----------------------------------------------------------------------
864      !
865      zfact   = 0.5_wp
866      !
867      IF(.NOT.ln_linssh ) THEN     ! need to re compute level distribution of isf fresh water
868         DO jj = 1,jpj
869            DO ji = 1,jpi
870               ikt = misfkt(ji,jj)
871               ikb = misfkt(ji,jj)
872               ! thickness of boundary layer at least the top level thickness
873               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
874
875               ! determine the deepest level influenced by the boundary layer
876               DO jk = ikt, mbkt(ji,jj)
877                  IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
878               END DO
879               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
880               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
881               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
882
883               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
884               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
885            END DO
886         END DO
887      END IF 
888      !
889      !==   ice shelf melting distributed over several levels   ==!
890      DO jj = 1,jpj
891         DO ji = 1,jpi
892               ikt = misfkt(ji,jj)
893               ikb = misfkb(ji,jj)
894               ! level fully include in the ice shelf boundary layer
895               DO jk = ikt, ikb - 1
896                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
897                    &              * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
898               END DO
899               ! level partially include in ice shelf boundary layer
900               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
901                    &            + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
902         END DO
903      END DO
904      !
905   END SUBROUTINE sbc_isf_div
906
907   !!======================================================================
908END MODULE sbcisf
Note: See TracBrowser for help on using the repository browser.