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/branches/UKMO/NEMO_4.0.2_mirror/src/OCE/SBC – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.2_mirror/src/OCE/SBC/sbcisf.F90 @ 12658

Last change on this file since 12658 was 12658, checked in by cguiavarch, 4 years ago

UKMO/NEMO_4.0.2_mirror : Remove SVN keywords.

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