source: NEMO/branches/2019/dev_r11756_SI3restart_XIOS/src/OCE/SBC/sbcisf.F90 @ 11837

Last change on this file since 11837 was 11837, checked in by andmirek, 18 months ago

ticket #2323 read SI3 restart with XIOS

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