source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/SBC/sbcisf.F90 @ 10023

Last change on this file since 10023 was 10023, checked in by gm, 2 years ago

#1911 (ENHANCE-04): RK3 branch - step II.2 bug correction in dynnxt + domvvl_RK3 creation

  • Property svn:keywords set to Id
File size: 45.4 KB
Line 
1MODULE sbcisf
2   !!======================================================================
3   !!                       ***  MODULE  sbcisf  ***
4   !! Surface module :  update surface ocean boundary condition under ice
5   !!                   shelf
6   !!======================================================================
7   !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
8   !!            X.X  !  2006-02  (C. Wang   ) Original code bg03
9   !!            3.4  !  2013-03  (P. Mathiot) Merging + parametrization
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   sbc_isf       : update sbc under ice shelf
14   !!----------------------------------------------------------------------
15   USE oce            ! ocean dynamics and tracers
16   USE dom_oce        ! ocean space and time domain
17   USE phycst         ! physical constants
18   USE eosbn2         ! equation of state
19   USE sbc_oce        ! surface boundary condition: ocean fields
20   USE zdfdrg         ! vertical physics: top/bottom drag coef.
21   !
22   USE in_out_manager ! I/O manager
23   USE iom            ! I/O library
24   USE fldread        ! read input field at current time step
25   USE lbclnk         !
26   USE lib_fortran    ! glob_sum
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   sbc_isf, sbc_isf_init, sbc_isf_div, sbc_isf_alloc  ! routine called in sbcmod and divhor
32
33   ! public in order to be able to output then
34
35   REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m]
36   INTEGER , PUBLIC ::   nn_isf                      !: flag to choose between explicit/param/specified 
37   INTEGER , PUBLIC ::   nn_isfblk                   !: flag to choose the bulk formulation to compute the ice shelf melting
38   INTEGER , PUBLIC ::   nn_gammablk                 !: flag to choose how the exchange coefficient is computed
39   REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient []
40   REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient []
41
42   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   misfkt   , misfkb        !: Level of ice shelf base
43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rzisf_tbl                !: depth of calving front (shallowest point) nn_isf ==2/3
44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rhisf_tbl, rhisf_tbl_0   !: thickness of tbl  [m]
45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   r1_hisf_tbl              !: 1/thickness of tbl
46   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ralpha                   !: proportion of bottom cell influenced by tbl
47   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   risfLeff                 !: effective length (Leff) BG03 nn_isf==2
48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ttbl, stbl, utbl, vtbl   !: top boundary layer variable at T point
49   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf                     !: net heat flux from ice shelf      [W/m2]
50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc     !: before and now T & S isf contents [K.m/s & PSU.m/s] 
51
52   LOGICAL, PUBLIC ::   l_isfcpl = .false.       !: isf recieved from oasis
53
54   REAL(wp)        , SAVE ::   rcp_isf  = 2000.0_wp     !: specific heat of ice shelf             [J/kg/K]
55   REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    !: heat diffusivity through the ice-shelf [m2/s]
56   REAL(wp), PUBLIC, SAVE ::   rho_isf  = 920.0_wp      !: volumic mass of ice shelf              [kg/m3]
57   REAL(wp), PUBLIC, SAVE ::   tsurf    = -20.0_wp      !: air temperature on top of ice shelf    [C]
58
59!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
60   CHARACTER(len=100), PUBLIC           :: cn_dirisf  = './' !: Root directory for location of ssr files
61   TYPE(FLD_N)       , PUBLIC           :: sn_fwfisf         !: information about the isf melting file to be read
62   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_fwfisf
63   TYPE(FLD_N)       , PUBLIC           :: sn_rnfisf         !: information about the isf melting param.   file to be read
64   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_rnfisf           
65   TYPE(FLD_N)       , PUBLIC           :: sn_depmax_isf     !: information about the grounding line depth file to be read
66   TYPE(FLD_N)       , PUBLIC           :: sn_depmin_isf     !: information about the calving   line depth file to be read
67   TYPE(FLD_N)       , PUBLIC           :: sn_Leff_isf       !: information about the effective length     file to be read
68   
69   !!----------------------------------------------------------------------
70   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
71   !! $Id$
72   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
73   !!----------------------------------------------------------------------
74CONTAINS
75 
76  SUBROUTINE sbc_isf( kt )
77      !!---------------------------------------------------------------------
78      !!                  ***  ROUTINE sbc_isf  ***
79      !!
80      !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf
81      !!              melting and freezing
82      !!
83      !! ** Method  :  4 parameterizations are available according to nn_isf
84      !!               nn_isf = 1 : Realistic ice_shelf formulation
85      !!                        2 : Beckmann & Goose parameterization
86      !!                        3 : Specified runoff in deptht (Mathiot & al. )
87      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf
88      !!----------------------------------------------------------------------
89      INTEGER, INTENT(in) ::   kt   ! ocean time step
90      !
91      INTEGER ::   ji, jj, jk   ! loop index
92      INTEGER ::   ikt, ikb     ! local integers
93      REAL(wp), DIMENSION(jpi,jpj) ::   zt_frz, zdep   ! freezing temperature (zt_frz) at depth (zdep)
94      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zqhcisf2d
95      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zfwfisf3d, zqhcisf3d, zqlatisf3d
96      !!---------------------------------------------------------------------
97      !
98      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN    ! compute salt and heat flux
99         !
100         SELECT CASE ( nn_isf )
101         CASE ( 1 )    ! realistic ice shelf formulation
102            ! compute T/S/U/V for the top boundary layer
103            CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')
104            CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T')
105            CALL sbc_isf_tbl(un(:,:,:)        ,utbl(:,:),'U')
106            CALL sbc_isf_tbl(vn(:,:,:)        ,vtbl(:,:),'V')
107            ! iom print
108            CALL iom_put('ttbl',ttbl(:,:))
109            CALL iom_put('stbl',stbl(:,:))
110            CALL iom_put('utbl',utbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:))
111            CALL iom_put('vtbl',vtbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:))
112            ! compute fwf and heat flux
113            ! compute fwf and heat flux
114            IF( .NOT.l_isfcpl ) THEN    ;   CALL sbc_isf_cav (kt)
115            ELSE                        ;   qisf(:,:)  = fwfisf(:,:) * rLfus    ! heat flux
116            ENDIF
117            !
118         CASE ( 2 )    ! Beckmann and Goosse parametrisation
119            stbl(:,:)   = soce
120            CALL sbc_isf_bg03(kt)
121            !
122         CASE ( 3 )    ! specified runoff in depth (Mathiot et al., XXXX in preparation)
123            ! specified runoff in depth (Mathiot et al., XXXX in preparation)
124            IF( .NOT.l_isfcpl ) THEN
125               CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
126               fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting)
127            ENDIF
128            qisf(:,:)   = fwfisf(:,:) * rLfus                   ! heat flux
129            stbl(:,:)   = soce
130            !
131         CASE ( 4 )    ! specified fwf and heat flux forcing beneath the ice shelf
132            !          ! specified fwf and heat flux forcing beneath the ice shelf
133            IF( .NOT.l_isfcpl ) THEN
134               CALL fld_read ( kt, nn_fsbc, sf_fwfisf   )
135               !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
136               fwfisf(:,:) = -sf_fwfisf(1)%fnow(:,:,1)            ! fwf
137            ENDIF
138            qisf(:,:)   = fwfisf(:,:) * rLfus                     ! heat flux
139            stbl(:,:)   = soce
140            !
141         END SELECT
142
143         ! compute tsc due to isf
144         ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU.
145         ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / rho0).
146         ! compute freezing point beneath ice shelf (or top cell if nn_isf = 3)
147         DO jj = 1,jpj
148            DO ji = 1,jpi
149               zdep(ji,jj)=gdepw_n(ji,jj,misfkt(ji,jj))
150            END DO
151         END DO
152         CALL eos_fzp( stbl(:,:), zt_frz(:,:), zdep(:,:) )
153         
154         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rho0_rcp - fwfisf(:,:) * zt_frz(:,:) * r1_rho0 !
155         risf_tsc(:,:,jp_sal) = 0.0_wp
156
157         ! lbclnk
158         CALL lbc_lnk_multi( risf_tsc(:,:,jp_tem), 'T', 1., risf_tsc(:,:,jp_sal), 'T', 1., fwfisf,'T', 1., qisf, 'T', 1.)
159         ! output
160         IF( iom_use('iceshelf_cea') )   CALL iom_put( 'iceshelf_cea', -fwfisf(:,:)                      )   ! isf mass flux
161         IF( iom_use('hflx_isf_cea') )   CALL iom_put( 'hflx_isf_cea', risf_tsc(:,:,jp_tem) * rho0 * rcp )   ! isf sensible+latent heat (W/m2)
162         IF( iom_use('qlatisf' ) )       CALL iom_put( 'qlatisf'     , qisf(:,:)                         )   ! isf latent heat
163         IF( iom_use('fwfisf'  ) )       CALL iom_put( 'fwfisf'      , fwfisf(:,:)                       )   ! isf mass flux (opposite sign)
164
165         ! Diagnostics
166         IF( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN
167            ALLOCATE( zfwfisf3d(jpi,jpj,jpk) , zqhcisf3d(jpi,jpj,jpk) , zqlatisf3d(jpi,jpj,jpk) )
168            ALLOCATE( zqhcisf2d(jpi,jpj) )
169            !
170            zfwfisf3d (:,:,:) = 0._wp                         ! 3d ice shelf melting (kg/m2/s)
171            zqhcisf3d (:,:,:) = 0._wp                         ! 3d heat content flux (W/m2)
172            zqlatisf3d(:,:,:) = 0._wp                         ! 3d ice shelf melting latent heat flux (W/m2)
173            zqhcisf2d (:,:)   = fwfisf(:,:) * zt_frz * rcp    ! 2d heat content flux (W/m2)
174            !
175            DO jj = 1,jpj
176               DO ji = 1,jpi
177                  ikt = misfkt(ji,jj)
178                  ikb = misfkb(ji,jj)
179                  DO jk = ikt, ikb - 1
180                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
181                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
182                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
183                  END DO
184                  zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj)   & 
185                     &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk)
186                  zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj)   & 
187                     &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk)
188                  zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj)   & 
189                     &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk)
190               END DO
191            END DO
192            !
193            CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:))
194            CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:))
195            CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:))
196            CALL iom_put('qhcisf'   , zqhcisf2d (:,:  ))
197            !
198            DEALLOCATE( zfwfisf3d, zqhcisf3d, zqlatisf3d )
199            DEALLOCATE( zqhcisf2d )
200         ENDIF
201         !
202      ENDIF
203
204      IF( kt == nit000 ) THEN                         !   set the forcing field at nit000 - 1    !
205         IF( ln_rstart .AND.    &                     ! Restart: read in restart file
206            &   iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN
207            IF(lwp) WRITE(numout,*) '          nit000-1 isf tracer content forcing fields read in the restart file'
208            CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:)         , ldxios = lrxios )   ! before salt content isf_tsc trend
209            CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b' , risf_tsc_b(:,:,jp_sal), ldxios = lrxios )   ! before salt content isf_tsc trend
210            CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b' , risf_tsc_b(:,:,jp_tem), ldxios = lrxios )   ! before salt content isf_tsc trend
211         ELSE
212            fwfisf_b(:,:)    = fwfisf(:,:)
213            risf_tsc_b(:,:,:)= risf_tsc(:,:,:)
214         ENDIF
215      ENDIF
216      !
217      IF( lrst_oce ) THEN
218         IF(lwp) WRITE(numout,*)
219         IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ',   &
220            &                    'at it= ', kt,' date= ', ndastp
221         IF(lwp) WRITE(numout,*) '~~~~'
222         IF( lwxios ) CALL iom_swap(      cwxios_context          )
223         CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:)         , ldxios = lwxios )
224         CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem), ldxios = lwxios )
225         CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal), ldxios = lwxios )
226         IF( lwxios ) CALL iom_swap(      cxios_context          )
227      ENDIF
228      !
229   END SUBROUTINE sbc_isf
230
231
232   INTEGER FUNCTION sbc_isf_alloc()
233      !!----------------------------------------------------------------------
234      !!               ***  FUNCTION sbc_isf_rnf_alloc  ***
235      !!----------------------------------------------------------------------
236      sbc_isf_alloc = 0       ! set to zero if no array to be allocated
237      IF( .NOT. ALLOCATED( qisf ) ) THEN
238         ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj)   , &
239               &    rhisf_tbl(jpi,jpj)    , r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj)  , &
240               &    ttbl(jpi,jpj)         , stbl(jpi,jpj)       , utbl(jpi,jpj)       , &
241               &    vtbl(jpi, jpj)        , risfLeff(jpi,jpj)   , rhisf_tbl_0(jpi,jpj), &
242               &    ralpha(jpi,jpj)       , misfkt(jpi,jpj)     , misfkb(jpi,jpj)     , &
243               &    STAT= sbc_isf_alloc )
244         !
245         IF( lk_mpp             )   CALL mpp_sum ( sbc_isf_alloc )
246         IF( sbc_isf_alloc /= 0 )   CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.')
247         !
248      ENDIF
249   END FUNCTION
250
251
252  SUBROUTINE sbc_isf_init
253      !!---------------------------------------------------------------------
254      !!                  ***  ROUTINE sbc_isf_init  ***
255      !!
256      !! ** Purpose : Initialisation of variables for iceshelf fluxes formulation
257      !!
258      !! ** Method  :  4 parameterizations are available according to nn_isf
259      !!               nn_isf = 1 : Realistic ice_shelf formulation
260      !!                        2 : Beckmann & Goose parameterization
261      !!                        3 : Specified runoff in deptht (Mathiot & al. )
262      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf
263      !!----------------------------------------------------------------------
264      INTEGER               :: ji, jj, jk           ! loop index
265      INTEGER               :: ik                   ! current level index
266      INTEGER               :: ikt, ikb             ! top and bottom level of the isf boundary layer
267      INTEGER               :: inum, ierror
268      INTEGER               :: ios                  ! Local integer output status for namelist read
269      REAL(wp)              :: zhk
270      CHARACTER(len=256)    :: cvarzisf, cvarhisf   ! name for isf file
271      CHARACTER(LEN=32 )    :: cvarLeff             ! variable name for efficient Length scale
272      !!----------------------------------------------------------------------
273      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf, &
274                         & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
275      !!----------------------------------------------------------------------
276
277      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs
278      READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
279901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp )
280
281      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
282      READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
283902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
284      IF(lwm) WRITE ( numond, namsbc_isf )
285
286      IF(lwp) WRITE(numout,*)
287      IF(lwp) WRITE(numout,*) 'sbc_isf_init : heat flux of the ice shelf'
288      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
289      IF(lwp) WRITE(numout,*) '   Namelist namsbc_isf :'
290      IF(lwp) WRITE(numout,*) '      type ice shelf melting/freezing         nn_isf      = ', nn_isf
291      IF(lwp) WRITE(numout,*) '      bulk formulation (nn_isf=1 only)        nn_isfblk   = ', nn_isfblk
292      IF(lwp) WRITE(numout,*) '      thickness of the top boundary layer     rn_hisf_tbl = ', rn_hisf_tbl
293      IF(lwp) WRITE(numout,*) '      gamma formulation                       nn_gammablk = ', nn_gammablk 
294      IF(lwp) WRITE(numout,*) '      gammat coefficient                      rn_gammat0  = ', rn_gammat0 
295      IF(lwp) WRITE(numout,*) '      gammas coefficient                      rn_gammas0  = ', rn_gammas0 
296      IF(lwp) WRITE(numout,*) '      top drag coef. used (from namdrg_top)   rn_Cd0      = ', r_Cdmin_top 
297
298
299                           !  1 = presence of ISF    2 = bg03 parametrisation
300                           !  3 = rnf file for isf   4 = ISF fwf specified
301                           !  option 1 and 4 need ln_isfcav = .true. (domzgr)
302      !
303      ! Allocate public variable
304      IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
305      !
306      ! initialisation
307      qisf    (:,:)    = 0._wp   ;   fwfisf  (:,:) = 0._wp
308      risf_tsc(:,:,:)  = 0._wp   ;   fwfisf_b(:,:) = 0._wp
309
310      SELECT CASE ( nn_isf )      ! define isf tbl tickness, top and bottom indice
311      !
312      CASE ( 1 ) 
313         IF(lwp) WRITE(numout,*)
314         IF(lwp) WRITE(numout,*) '      ==>>>   presence of under iceshelf seas (nn_isf = 1)'
315         rhisf_tbl(:,:) = rn_hisf_tbl
316         misfkt   (:,:) = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
317         !
318      CASE ( 2 , 3 )
319         IF( .NOT.l_isfcpl ) THEN
320             ALLOCATE( sf_rnfisf(1), STAT=ierror )
321             ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
322             CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
323          ENDIF
324          !  read effective lenght (BG03)
325          IF( nn_isf == 2 ) THEN
326            IF(lwp) WRITE(numout,*)
327            IF(lwp) WRITE(numout,*) '      ==>>>   bg03 parametrisation (nn_isf = 2)'
328            CALL iom_open( sn_Leff_isf%clname, inum )
329            cvarLeff = TRIM(sn_Leff_isf%clvar)
330            CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
331            CALL iom_close(inum)
332            !
333            risfLeff = risfLeff*1000.0_wp           !: convertion in m
334         ELSE
335            IF(lwp) WRITE(numout,*)
336            IF(lwp) WRITE(numout,*) '      ==>>>   rnf file for isf (nn_isf = 3)'
337         ENDIF
338         ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
339         CALL iom_open( sn_depmax_isf%clname, inum )
340         cvarhisf = TRIM(sn_depmax_isf%clvar)
341         CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
342         CALL iom_close(inum)
343         !
344         CALL iom_open( sn_depmin_isf%clname, inum )
345         cvarzisf = TRIM(sn_depmin_isf%clvar)
346         CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
347         CALL iom_close(inum)
348         !
349         rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer
350
351         !! compute first level of the top boundary layer
352         DO ji = 1, jpi
353            DO jj = 1, jpj
354                ik = 2
355!!gm potential bug: use gdepw_0 not _n
356                DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO
357                misfkt(ji,jj) = ik-1
358            END DO
359         END DO
360         !
361      CASE ( 4 ) 
362         IF(lwp) WRITE(numout,*)
363         IF(lwp) WRITE(numout,*) '      ==>>>   specified fresh water flux in ISF (nn_isf = 4)'
364         ! as in nn_isf == 1
365         rhisf_tbl(:,:) = rn_hisf_tbl
366         misfkt   (:,:) = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
367         !
368         ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
369         IF( .NOT.l_isfcpl ) THEN
370           ALLOCATE( sf_fwfisf(1), STAT=ierror )
371           ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
372           CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
373         ENDIF
374         !
375      CASE DEFAULT
376         CALL ctl_stop( 'sbc_isf_init: wrong value of nn_isf' )
377      END SELECT
378         
379      rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
380
381      ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
382      DO jj = 1,jpj
383         DO ji = 1,jpi
384            ikt = misfkt(ji,jj)
385            ikb = misfkt(ji,jj)
386            ! thickness of boundary layer at least the top level thickness
387            rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
388
389            ! determine the deepest level influenced by the boundary layer
390            DO jk = ikt+1, mbkt(ji,jj)
391               IF( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) )   ikb = jk
392            END DO
393            rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
394            misfkb(ji,jj) = ikb                                                   ! last wet level of the tbl
395            r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
396
397            zhk           = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1
398            ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
399         END DO
400      END DO
401
402      IF( lwxios ) THEN
403          CALL iom_set_rstw_var_active('fwf_isf_b')
404          CALL iom_set_rstw_var_active('isf_hc_b')
405          CALL iom_set_rstw_var_active('isf_sc_b')
406      ENDIF
407
408
409  END SUBROUTINE sbc_isf_init
410
411
412  SUBROUTINE sbc_isf_bg03(kt)
413      !!---------------------------------------------------------------------
414      !!                  ***  ROUTINE sbc_isf_bg03  ***
415      !!
416      !! ** Purpose : add net heat and fresh water flux from ice shelf melting
417      !!          into the adjacent ocean
418      !!
419      !! ** Method  :   See reference
420      !!
421      !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
422      !!         interaction for climate models", Ocean Modelling 5(2003) 157-170.
423      !!         (hereafter BG)
424      !! History :  06-02  (C. Wang) Original code
425      !!----------------------------------------------------------------------
426      INTEGER, INTENT ( in ) :: kt
427      !
428      INTEGER  :: ji, jj, jk ! dummy loop index
429      INTEGER  :: ik         ! current level
430      REAL(wp) :: zt_sum     ! sum of the temperature between 200m and 600m
431      REAL(wp) :: zt_ave     ! averaged temperature between 200m and 600m
432      REAL(wp) :: zt_frz     ! freezing point temperature at depth z
433      REAL(wp) :: zpress     ! pressure to compute the freezing point in depth
434      !!----------------------------------------------------------------------
435      !
436      DO ji = 1, jpi
437         DO jj = 1, jpj
438            ik = misfkt(ji,jj)
439            !! Initialize arrays to 0 (each step)
440            zt_sum = 0.e0_wp
441            IF ( ik > 1 ) THEN
442               ! 1. -----------the average temperature between 200m and 600m ---------------------
443               DO jk = misfkt(ji,jj),misfkb(ji,jj)
444                  ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!)
445                  ! after verif with UNESCO, wrong sign in BG eq. 2
446                  ! Calculate freezing temperature
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) = - rho0 * 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) * r1_Lfus                        ! 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 * rho0
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 * rho0
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(:,:), ht_isf(:,:) )
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*rho0*(ttbl(ji,jj)-zfrz(ji,jj))
528                  zfwflx(ji,jj) = - zhtflx(ji,jj) * r1_Lfus
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 = rho0_rcp   * zgammat(ji,jj)
546                  zeps2 = rho0*rLfus * zgammas(ji,jj)
547                  zeps3 = rho_isf*rcp_isf*rkappa / MAX( ht_isf(ji,jj) , zeps )
548                  zeps4 = zlamb2 + zlamb3 * ht_isf(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) = rho0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps)
569                  zhtflx(ji,jj) = zgammat(ji,jj) * rho0 * 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)
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( 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      !!----------------------------------------------------------------------
745     
746      ! initialisation
747      pvarout(:,:)=0._wp
748   
749      SELECT CASE ( cd_ptin )
750      CASE ( 'U' ) ! compute U in the top boundary layer at T- point
751         DO jj = 1,jpj
752            DO ji = 1,jpi
753               ikt = miku(ji,jj) ; ikb = miku(ji,jj)
754               ! thickness of boundary layer at least the top level thickness
755               zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u_n(ji,jj,ikt) )
756
757               ! determine the deepest level influenced by the boundary layer
758               DO jk = ikt+1, mbku(ji,jj)
759                  IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
760               END DO
761               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
762
763               ! level fully include in the ice shelf boundary layer
764               DO jk = ikt, ikb - 1
765                  ze3 = e3u_n(ji,jj,jk)
766                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
767               END DO
768
769               ! level partially include in ice shelf boundary layer
770               zhk = SUM( e3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
771               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
772            END DO
773         END DO
774         DO jj = 2, jpj
775            DO ji = 2, jpi
776!!gm a wet-point only average should be used here !!!
777               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj))
778            END DO
779         END DO
780         CALL lbc_lnk(pvarout,'T',-1.)
781     
782      CASE ( 'V' ) ! compute V in the top boundary layer at T- point
783         DO jj = 1,jpj
784            DO ji = 1,jpi
785               ikt = mikv(ji,jj) ; ikb = mikv(ji,jj)
786               ! thickness of boundary layer at least the top level thickness
787               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v_n(ji,jj,ikt))
788
789               ! determine the deepest level influenced by the boundary layer
790               DO jk = ikt+1, mbkv(ji,jj)
791                  IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
792               END DO
793               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
794
795               ! level fully include in the ice shelf boundary layer
796               DO jk = ikt, ikb - 1
797                  ze3 = e3v_n(ji,jj,jk)
798                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
799               END DO
800
801               ! level partially include in ice shelf boundary layer
802               zhk = SUM( e3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
803               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
804            END DO
805         END DO
806         DO jj = 2, jpj
807            DO ji = 2, jpi
808!!gm a wet-point only average should be used here !!!
809               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1))
810            END DO
811         END DO
812         CALL lbc_lnk(pvarout,'T',-1.)
813
814      CASE ( 'T' ) ! compute T in the top boundary layer at T- point
815         DO jj = 1,jpj
816            DO ji = 1,jpi
817               ikt = misfkt(ji,jj)
818               ikb = misfkb(ji,jj)
819
820               ! level fully include in the ice shelf boundary layer
821               DO jk = ikt, ikb - 1
822                  ze3 = e3t_n(ji,jj,jk)
823                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
824               END DO
825
826               ! level partially include in ice shelf boundary layer
827               zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
828               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
829            END DO
830         END DO
831      END SELECT
832      !
833      ! mask mean tbl value
834      pvarout(:,:) = pvarout(:,:) * ssmask(:,:)
835      !
836   END SUBROUTINE sbc_isf_tbl
837     
838
839   SUBROUTINE sbc_isf_div( phdivn )
840      !!----------------------------------------------------------------------
841      !!                  ***  SUBROUTINE sbc_isf_div  ***
842      !!       
843      !! ** Purpose :   update the horizontal divergence with the runoff inflow
844      !!
845      !! ** Method  :   
846      !!                CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the
847      !!                          divergence and expressed in m/s
848      !!
849      !! ** Action  :   phdivn   decreased by the runoff inflow
850      !!----------------------------------------------------------------------
851      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdivn   ! horizontal divergence
852      !
853      INTEGER  ::   ji, jj, jk   ! dummy loop indices
854      INTEGER  ::   ikt, ikb 
855      REAL(wp) ::   zhk
856      REAL(wp) ::   zfact     ! local scalar
857      !!----------------------------------------------------------------------
858      !
859      zfact   = 0.5_wp
860      !
861      IF(.NOT.ln_linssh ) THEN     ! need to re compute level distribution of isf fresh water
862         DO jj = 1,jpj
863            DO ji = 1,jpi
864               ikt = misfkt(ji,jj)
865               ikb = misfkt(ji,jj)
866               ! thickness of boundary layer at least the top level thickness
867               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
868
869               ! determine the deepest level influenced by the boundary layer
870               DO jk = ikt, mbkt(ji,jj)
871                  IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
872               END DO
873               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
874               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
875               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
876
877               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
878               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
879            END DO
880         END DO
881      END IF 
882      !
883      !==   ice shelf melting distributed over several levels   ==!
884      DO jj = 1,jpj
885         DO ji = 1,jpi
886               ikt = misfkt(ji,jj)
887               ikb = misfkb(ji,jj)
888               ! level fully include in the ice shelf boundary layer
889               DO jk = ikt, ikb - 1
890                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
891                    &              * r1_hisf_tbl(ji,jj) * r1_rho0 * zfact
892               END DO
893               ! level partially include in ice shelf boundary layer
894               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
895                    &            + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rho0 * zfact * ralpha(ji,jj) 
896         END DO
897      END DO
898      !
899   END SUBROUTINE sbc_isf_div
900
901   !!======================================================================
902END MODULE sbcisf
Note: See TracBrowser for help on using the repository browser.