source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 9124

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

dev_merge_2017: ln_timing instead of nn_timing + restricted timing to nemo_init and routine called by step in OPA_SRC

  • Property svn:keywords set to Id
File size: 43.6 KB
Line 
1MODULE sbcisf
2   !!======================================================================
3   !!                       ***  MODULE  sbcisf  ***
4   !! Surface module :  update surface ocean boundary condition under ice
5   !!                   shelf
6   !!======================================================================
7   !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
8   !!            X.X  !  2006-02  (C. Wang   ) Original code bg03
9   !!            3.4  !  2013-03  (P. Mathiot) Merging + parametrization
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   sbc_isf       : update sbc under ice shelf
14   !!----------------------------------------------------------------------
15   USE oce            ! ocean dynamics and tracers
16   USE dom_oce        ! ocean space and time domain
17   USE phycst         ! physical constants
18   USE eosbn2         ! equation of state
19   USE sbc_oce        ! surface boundary condition: ocean fields
20   USE zdfdrg         ! vertical physics: top/bottom drag coef.
21   !
22   USE in_out_manager ! I/O manager
23   USE iom            ! I/O library
24   USE fldread        ! read input field at current time step
25   USE lbclnk         !
26   USE lib_fortran    ! glob_sum
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   sbc_isf, sbc_isf_init, sbc_isf_div, sbc_isf_alloc  ! routine called in sbcmod and divhor
32
33   ! public in order to be able to output then
34
35   REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m]
36   INTEGER , PUBLIC ::   nn_isf                      !: flag to choose between explicit/param/specified 
37   INTEGER , PUBLIC ::   nn_isfblk                   !: flag to choose the bulk formulation to compute the ice shelf melting
38   INTEGER , PUBLIC ::   nn_gammablk                 !: flag to choose how the exchange coefficient is computed
39   REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient []
40   REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient []
41
42   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   misfkt   , misfkb        !: Level of ice shelf base
43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rzisf_tbl                !: depth of calving front (shallowest point) nn_isf ==2/3
44   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rhisf_tbl, rhisf_tbl_0   !: thickness of tbl  [m]
45   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   r1_hisf_tbl              !: 1/thickness of tbl
46   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ralpha                   !: proportion of bottom cell influenced by tbl
47   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   risfLeff                 !: effective length (Leff) BG03 nn_isf==2
48   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ttbl, stbl, utbl, vtbl   !: top boundary layer variable at T point
49   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf                     !: net heat flux from ice shelf      [W/m2]
50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc     !: before and now T & S isf contents [K.m/s & PSU.m/s] 
51
52   LOGICAL, PUBLIC ::   l_isfcpl = .false.       !: isf recieved from oasis
53
54   REAL(wp), PUBLIC, SAVE ::   rcpi     = 2000.0_wp     !: specific heat of ice shelf             [J/kg/K]
55   REAL(wp), PUBLIC, SAVE ::   rkappa   = 1.54e-6_wp    !: heat diffusivity through the ice-shelf [m2/s]
56   REAL(wp), PUBLIC, SAVE ::   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/OPA 4.0 , LOCEAN-IPSL (2017)
72   !! $Id$
73   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
74   !!----------------------------------------------------------------------
75CONTAINS
76 
77  SUBROUTINE sbc_isf( kt )
78      !!---------------------------------------------------------------------
79      !!                  ***  ROUTINE sbc_isf  ***
80      !!
81      !! ** Purpose : Compute Salt and Heat fluxes related to ice_shelf
82      !!              melting and freezing
83      !!
84      !! ** Method  :  4 parameterizations are available according to nn_isf
85      !!               nn_isf = 1 : Realistic ice_shelf formulation
86      !!                        2 : Beckmann & Goose parameterization
87      !!                        3 : Specified runoff in deptht (Mathiot & al. )
88      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf
89      !!----------------------------------------------------------------------
90      INTEGER, INTENT(in) ::   kt   ! ocean time step
91      !
92      INTEGER ::   ji, jj, jk   ! loop index
93      INTEGER ::   ikt, ikb     ! local integers
94      REAL(wp), DIMENSION(jpi,jpj) ::   zt_frz, zdep   ! freezing temperature (zt_frz) at depth (zdep)
95      REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   zqhcisf2d
96      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zfwfisf3d, zqhcisf3d, zqlatisf3d
97      !!---------------------------------------------------------------------
98      !
99      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN    ! compute salt and heat flux
100         !
101         SELECT CASE ( nn_isf )
102         CASE ( 1 )    ! realistic ice shelf formulation
103            ! compute T/S/U/V for the top boundary layer
104            CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')
105            CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T')
106            CALL sbc_isf_tbl(un(:,:,:)        ,utbl(:,:),'U')
107            CALL sbc_isf_tbl(vn(:,:,:)        ,vtbl(:,:),'V')
108            ! iom print
109            CALL iom_put('ttbl',ttbl(:,:))
110            CALL iom_put('stbl',stbl(:,:))
111            CALL iom_put('utbl',utbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:))
112            CALL iom_put('vtbl',vtbl(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:))
113            ! compute fwf and heat flux
114            ! compute fwf and heat flux
115            IF( .NOT.l_isfcpl ) THEN    ;   CALL sbc_isf_cav (kt)
116            ELSE                        ;   qisf(:,:)  = fwfisf(:,:) * rlfusisf  ! heat        flux
117            ENDIF
118            !
119         CASE ( 2 )    ! Beckmann and Goosse parametrisation
120            stbl(:,:)   = soce
121            CALL sbc_isf_bg03(kt)
122            !
123         CASE ( 3 )    ! specified runoff in depth (Mathiot et al., XXXX in preparation)
124            ! specified runoff in depth (Mathiot et al., XXXX in preparation)
125            IF( .NOT.l_isfcpl ) THEN
126               CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
127               fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting)
128            ENDIF
129            qisf(:,:)   = fwfisf(:,:) * rlfusisf             ! heat flux
130            stbl(:,:)   = soce
131            !
132         CASE ( 4 )    ! specified fwf and heat flux forcing beneath the ice shelf
133            !          ! specified fwf and heat flux forcing beneath the ice shelf
134            IF( .NOT.l_isfcpl ) THEN
135               CALL fld_read ( kt, nn_fsbc, sf_fwfisf   )
136               !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
137               fwfisf(:,:) = -sf_fwfisf(1)%fnow(:,:,1)            ! fwf
138            ENDIF
139            qisf(:,:)   = fwfisf(:,:) * rlfusisf               ! heat flux
140            stbl(:,:)   = soce
141            !
142         END SELECT
143
144         ! compute tsc due to isf
145         ! isf melting implemented as a volume flux and we assume that melt water is at 0 PSU.
146         ! WARNING water add at temp = 0C, need to add a correction term (fwfisf * tfreez / 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( risf_tsc(:,:,jp_tem), 'T', 1., risf_tsc(:,:,jp_sal), 'T', 1., fwfisf,'T', 1., qisf, 'T', 1.)
160         ! output
161         IF( iom_use('iceshelf_cea') )   CALL iom_put( 'iceshelf_cea', -fwfisf(:,:)                      )   ! isf mass flux
162         IF( iom_use('hflx_isf_cea') )   CALL iom_put( 'hflx_isf_cea', risf_tsc(:,:,jp_tem) * rau0 * rcp )   ! isf sensible+latent heat (W/m2)
163         IF( iom_use('qlatisf' ) )       CALL iom_put( 'qlatisf'     , qisf(:,:)                         )   ! isf latent heat
164         IF( iom_use('fwfisf'  ) )       CALL iom_put( 'fwfisf'      , fwfisf(:,:)                       )   ! isf mass flux (opposite sign)
165
166         ! Diagnostics
167         IF( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN
168            ALLOCATE( zfwfisf3d(jpi,jpj,jpk) , zqhcisf3d(jpi,jpj,jpk) , zqlatisf3d(jpi,jpj,jpk) )
169            ALLOCATE( zqhcisf2d(jpi,jpj) )
170            !
171            zfwfisf3d (:,:,:) = 0._wp                         ! 3d ice shelf melting (kg/m2/s)
172            zqhcisf3d (:,:,:) = 0._wp                         ! 3d heat content flux (W/m2)
173            zqlatisf3d(:,:,:) = 0._wp                         ! 3d ice shelf melting latent heat flux (W/m2)
174            zqhcisf2d (:,:)   = fwfisf(:,:) * zt_frz * rcp    ! 2d heat content flux (W/m2)
175            !
176            DO jj = 1,jpj
177               DO ji = 1,jpi
178                  ikt = misfkt(ji,jj)
179                  ikb = misfkb(ji,jj)
180                  DO jk = ikt, ikb - 1
181                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
182                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
183                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * e3t_n(ji,jj,jk)
184                  END DO
185                  zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj)   & 
186                     &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk)
187                  zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj)   & 
188                     &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk)
189                  zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj)   & 
190                     &                                                                   * ralpha(ji,jj) * e3t_n(ji,jj,jk)
191               END DO
192            END DO
193            !
194            CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:))
195            CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:))
196            CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:))
197            CALL iom_put('qhcisf'   , zqhcisf2d (:,:  ))
198            !
199            DEALLOCATE( zfwfisf3d, zqhcisf3d, zqlatisf3d )
200            DEALLOCATE( zqhcisf2d )
201         ENDIF
202         !
203      ENDIF
204
205      IF( kt == nit000 ) THEN                         !   set the forcing field at nit000 - 1    !
206         IF( ln_rstart .AND.    &                     ! Restart: read in restart file
207            &   iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN
208            IF(lwp) WRITE(numout,*) '          nit000-1 isf tracer content forcing fields read in the restart file'
209            CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:) )   ! before salt content isf_tsc trend
210            CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b', risf_tsc_b(:,:,jp_sal) )   ! before salt content isf_tsc trend
211            CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b', risf_tsc_b(:,:,jp_tem) )   ! before salt content isf_tsc trend
212         ELSE
213            fwfisf_b(:,:)    = fwfisf(:,:)
214            risf_tsc_b(:,:,:)= risf_tsc(:,:,:)
215         ENDIF
216      ENDIF
217      !
218      IF( lrst_oce ) THEN
219         IF(lwp) WRITE(numout,*)
220         IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ',   &
221            &                    'at it= ', kt,' date= ', ndastp
222         IF(lwp) WRITE(numout,*) '~~~~'
223         CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:) )
224         CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem) )
225         CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal) )
226      ENDIF
227      !
228   END SUBROUTINE sbc_isf
229
230
231   INTEGER FUNCTION sbc_isf_alloc()
232      !!----------------------------------------------------------------------
233      !!               ***  FUNCTION sbc_isf_rnf_alloc  ***
234      !!----------------------------------------------------------------------
235      sbc_isf_alloc = 0       ! set to zero if no array to be allocated
236      IF( .NOT. ALLOCATED( qisf ) ) THEN
237         ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj)   , &
238               &    rhisf_tbl(jpi,jpj)    , r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj)  , &
239               &    ttbl(jpi,jpj)         , stbl(jpi,jpj)       , utbl(jpi,jpj)       , &
240               &    vtbl(jpi, jpj)        , risfLeff(jpi,jpj)   , rhisf_tbl_0(jpi,jpj), &
241               &    ralpha(jpi,jpj)       , misfkt(jpi,jpj)     , misfkb(jpi,jpj)     , &
242               &    STAT= sbc_isf_alloc )
243         !
244         IF( lk_mpp             )   CALL mpp_sum ( sbc_isf_alloc )
245         IF( sbc_isf_alloc /= 0 )   CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.')
246         !
247      ENDIF
248   END FUNCTION
249
250
251  SUBROUTINE sbc_isf_init
252      !!---------------------------------------------------------------------
253      !!                  ***  ROUTINE sbc_isf_init  ***
254      !!
255      !! ** Purpose : Initialisation of variables for iceshelf fluxes formulation
256      !!
257      !! ** Method  :  4 parameterizations are available according to nn_isf
258      !!               nn_isf = 1 : Realistic ice_shelf formulation
259      !!                        2 : Beckmann & Goose parameterization
260      !!                        3 : Specified runoff in deptht (Mathiot & al. )
261      !!                        4 : specified fwf and heat flux forcing beneath the ice shelf
262      !!----------------------------------------------------------------------
263      INTEGER               :: ji, jj, jk           ! loop index
264      INTEGER               :: ik                   ! current level index
265      INTEGER               :: ikt, ikb             ! top and bottom level of the isf boundary layer
266      INTEGER               :: inum, ierror
267      INTEGER               :: ios                  ! Local integer output status for namelist read
268      REAL(wp)              :: zhk
269      CHARACTER(len=256)    :: cvarzisf, cvarhisf   ! name for isf file
270      CHARACTER(LEN=32 )    :: cvarLeff             ! variable name for efficient Length scale
271      !!----------------------------------------------------------------------
272      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, rn_gammat0, rn_gammas0, nn_gammablk, nn_isf, &
273                         & sn_fwfisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
274      !!----------------------------------------------------------------------
275
276      REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs
277      READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
278901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp )
279
280      REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
281      READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
282902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
283      IF(lwm) WRITE ( numond, namsbc_isf )
284
285      IF ( lwp ) WRITE(numout,*)
286      IF ( lwp ) WRITE(numout,*) 'sbc_isf_init : heat flux of the ice shelf'
287      IF ( lwp ) WRITE(numout,*) '~~~~~~~~~~~'
288      IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf
289      IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk
290      IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl
291      IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk 
292      IF ( lwp ) WRITE(numout,*) '        rn_gammat0  = ', rn_gammat0 
293      IF ( lwp ) WRITE(numout,*) '        rn_gammas0  = ', rn_gammas0 
294      IF ( lwp ) WRITE(numout,*) '        rn_Cd0      = ', r_Cdmin_top 
295      !
296      ! Allocate public variable
297      IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
298      !
299      ! initialisation
300      qisf    (:,:)    = 0._wp   ;   fwfisf  (:,:) = 0._wp
301      risf_tsc(:,:,:)  = 0._wp   ;   fwfisf_b(:,:) = 0._wp
302      !
303      ! define isf tbl tickness, top and bottom indice
304      SELECT CASE ( nn_isf )
305      CASE ( 1 ) 
306         rhisf_tbl(:,:) = rn_hisf_tbl
307         misfkt   (:,:) = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
308
309      CASE ( 2 , 3 )
310         IF( .NOT.l_isfcpl ) THEN
311             ALLOCATE( sf_rnfisf(1), STAT=ierror )
312             ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
313             CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
314          ENDIF
315          !  read effective lenght (BG03)
316          IF (nn_isf == 2) THEN
317            CALL iom_open( sn_Leff_isf%clname, inum )
318            cvarLeff = TRIM(sn_Leff_isf%clvar)
319            CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
320            CALL iom_close(inum)
321            !
322            risfLeff = risfLeff*1000.0_wp           !: convertion in m
323          END IF
324         ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
325         CALL iom_open( sn_depmax_isf%clname, inum )
326         cvarhisf = TRIM(sn_depmax_isf%clvar)
327         CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
328         CALL iom_close(inum)
329         !
330         CALL iom_open( sn_depmin_isf%clname, inum )
331         cvarzisf = TRIM(sn_depmin_isf%clvar)
332         CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
333         CALL iom_close(inum)
334         !
335         rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer
336
337         !! compute first level of the top boundary layer
338         DO ji = 1, jpi
339            DO jj = 1, jpj
340                ik = 2
341!!gm potential bug: use gdepw_0 not _n
342                DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO
343                misfkt(ji,jj) = ik-1
344            END DO
345         END DO
346
347      CASE ( 4 ) 
348         ! as in nn_isf == 1
349         rhisf_tbl(:,:) = rn_hisf_tbl
350         misfkt   (:,:) = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
351         
352         ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
353         IF( .NOT.l_isfcpl ) THEN
354           ALLOCATE( sf_fwfisf(1), STAT=ierror )
355           ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
356           CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
357         ENDIF
358
359      END SELECT
360         
361      rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
362
363      ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
364      DO jj = 1,jpj
365         DO ji = 1,jpi
366            ikt = misfkt(ji,jj)
367            ikb = misfkt(ji,jj)
368            ! thickness of boundary layer at least the top level thickness
369            rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
370
371            ! determine the deepest level influenced by the boundary layer
372            DO jk = ikt+1, mbkt(ji,jj)
373               IF( (SUM(e3t_n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) )   ikb = jk
374            END DO
375            rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.
376            misfkb(ji,jj) = ikb                                                   ! last wet level of the tbl
377            r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
378
379            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
380            ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
381         END DO
382      END DO
383
384  END SUBROUTINE sbc_isf_init
385
386
387  SUBROUTINE sbc_isf_bg03(kt)
388      !!---------------------------------------------------------------------
389      !!                  ***  ROUTINE sbc_isf_bg03  ***
390      !!
391      !! ** Purpose : add net heat and fresh water flux from ice shelf melting
392      !!          into the adjacent ocean
393      !!
394      !! ** Method  :   See reference
395      !!
396      !! ** Reference : Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
397      !!         interaction for climate models", Ocean Modelling 5(2003) 157-170.
398      !!         (hereafter BG)
399      !! History :  06-02  (C. Wang) Original code
400      !!----------------------------------------------------------------------
401      INTEGER, INTENT ( in ) :: kt
402      !
403      INTEGER  :: ji, jj, jk ! dummy loop index
404      INTEGER  :: ik         ! current level
405      REAL(wp) :: zt_sum     ! sum of the temperature between 200m and 600m
406      REAL(wp) :: zt_ave     ! averaged temperature between 200m and 600m
407      REAL(wp) :: zt_frz     ! freezing point temperature at depth z
408      REAL(wp) :: zpress     ! pressure to compute the freezing point in depth
409      !!----------------------------------------------------------------------
410      !
411      DO ji = 1, jpi
412         DO jj = 1, jpj
413            ik = misfkt(ji,jj)
414            !! Initialize arrays to 0 (each step)
415            zt_sum = 0.e0_wp
416            IF ( ik > 1 ) THEN
417               ! 1. -----------the average temperature between 200m and 600m ---------------------
418               DO jk = misfkt(ji,jj),misfkb(ji,jj)
419                  ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!)
420                  ! after verif with UNESCO, wrong sign in BG eq. 2
421                  ! Calculate freezing temperature
422                  CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 
423                  zt_sum = zt_sum + (tsn(ji,jj,jk,jp_tem)-zt_frz) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk)  ! sum temp
424               END DO
425               zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
426               ! 2. ------------Net heat flux and fresh water flux due to the ice shelf
427               ! For those corresponding to zonal boundary   
428               qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
429                           & * r1_e1e2t(ji,jj) * tmask(ji,jj,jk)
430             
431               fwfisf(ji,jj) = qisf(ji,jj) / rlfusisf          !fresh water flux kg/(m2s)                 
432               fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
433               !add to salinity trend
434            ELSE
435               qisf(ji,jj) = 0._wp   ;   fwfisf(ji,jj) = 0._wp
436            END IF
437         END DO
438      END DO
439      !
440  END SUBROUTINE sbc_isf_bg03
441
442
443  SUBROUTINE sbc_isf_cav( kt )
444      !!---------------------------------------------------------------------
445      !!                     ***  ROUTINE sbc_isf_cav  ***
446      !!
447      !! ** Purpose :   handle surface boundary condition under ice shelf
448      !!
449      !! ** Method  : -
450      !!
451      !! ** Action  :   utau, vtau : remain unchanged
452      !!                taum, wndm : remain unchanged
453      !!                qns        : update heat flux below ice shelf
454      !!                emp, emps  : update freshwater flux below ice shelf
455      !!---------------------------------------------------------------------
456      INTEGER, INTENT(in) ::   kt   ! ocean time step
457      !
458      INTEGER  ::   ji, jj     ! dummy loop indices
459      INTEGER  ::   nit
460      LOGICAL  ::   lit
461      REAL(wp) ::   zlamb1, zlamb2, zlamb3
462      REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7
463      REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac
464      REAL(wp) ::   zeps = 1.e-20_wp       
465      REAL(wp) ::   zerr
466      REAL(wp), DIMENSION(jpi,jpj) ::   zfrz
467      REAL(wp), DIMENSION(jpi,jpj) ::   zgammat, zgammas 
468      REAL(wp), DIMENSION(jpi,jpj) ::   zfwflx, zhtflx, zhtflx_b
469      !!---------------------------------------------------------------------
470      !
471      ! coeficient for linearisation of potential tfreez
472      ! Crude approximation for pressure (but commonly used)
473      zlamb1 =-0.0573_wp
474      zlamb2 = 0.0832_wp
475      zlamb3 =-7.53e-08_wp * grav * rau0
476      !
477      ! initialisation
478      zgammat(:,:) = rn_gammat0 ; zgammas (:,:) = rn_gammas0
479      zhtflx (:,:) = 0.0_wp     ; zhtflx_b(:,:) = 0.0_wp   
480      zfwflx (:,:) = 0.0_wp
481
482      ! compute ice shelf melting
483      nit = 1 ; lit = .TRUE.
484      DO WHILE ( lit )    ! maybe just a constant number of iteration as in blk_core is fine
485         SELECT CASE ( nn_isfblk )
486         CASE ( 1 )   !  ISOMIP formulation (2 equations) for volume flux (Hunter et al., 2006)
487            ! Calculate freezing temperature
488            CALL eos_fzp( stbl(:,:), zfrz(:,:), risfdep(:,:) )
489
490            ! compute gammat every where (2d)
491            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
492           
493            ! compute upward heat flux zhtflx and upward water flux zwflx
494            DO jj = 1, jpj
495               DO ji = 1, jpi
496                  zhtflx(ji,jj) =   zgammat(ji,jj)*rcp*rau0*(ttbl(ji,jj)-zfrz(ji,jj))
497                  zfwflx(ji,jj) = - zhtflx(ji,jj)/rlfusisf
498               END DO
499            END DO
500
501            ! Compute heat flux and upward fresh water flux
502            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
503            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
504
505         CASE ( 2 )  ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015)
506            ! compute gammat every where (2d)
507            CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx)
508
509            ! compute upward heat flux zhtflx and upward water flux zwflx
510            ! Resolution of a 2d equation from equation 21, 22 and 23 to find Sb (Asay-Davis et al., 2015)
511            DO jj = 1, jpj
512               DO ji = 1, jpi
513                  ! compute coeficient to solve the 2nd order equation
514                  zeps1 = rcp*rau0*zgammat(ji,jj)
515                  zeps2 = rlfusisf*rau0*zgammas(ji,jj)
516                  zeps3 = rhoisf*rcpi*rkappa/MAX(risfdep(ji,jj),zeps)
517                  zeps4 = zlamb2+zlamb3*risfdep(ji,jj)
518                  zeps6 = zeps4-ttbl(ji,jj)
519                  zeps7 = zeps4-tsurf
520                  zaqe  = zlamb1 * (zeps1 + zeps3)
521                  zaqer = 0.5_wp/MIN(zaqe,-zeps)
522                  zbqe  = zeps1*zeps6+zeps3*zeps7-zeps2
523                  zcqe  = zeps2*stbl(ji,jj)
524                  zdis  = zbqe*zbqe-4.0_wp*zaqe*zcqe               
525
526                  ! Presumably zdis can never be negative because gammas is very small compared to gammat
527                  ! compute s freeze
528                  zsfrz=(-zbqe-SQRT(zdis))*zaqer
529                  IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer
530
531                  ! compute t freeze (eq. 22)
532                  zfrz(ji,jj)=zeps4+zlamb1*zsfrz
533 
534                  ! zfwflx is upward water flux
535                  ! zhtflx is upward heat flux (out of ocean)
536                  ! compute the upward water and heat flux (eq. 28 and eq. 29)
537                  zfwflx(ji,jj) = rau0 * zgammas(ji,jj) * (zsfrz-stbl(ji,jj)) / MAX(zsfrz,zeps)
538                  zhtflx(ji,jj) = zgammat(ji,jj) * rau0 * rcp * (ttbl(ji,jj) - zfrz(ji,jj) ) 
539               END DO
540            END DO
541
542            ! compute heat and water flux
543            qisf  (:,:) = - zhtflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
544            fwfisf(:,:) =   zfwflx(:,:) * (1._wp - tmask(:,:,1)) * ssmask(:,:)
545
546         END SELECT
547
548         ! define if we need to iterate (nn_gammablk 0/1 do not need iteration)
549         IF ( nn_gammablk <  2 ) THEN ; lit = .FALSE.
550         ELSE                           
551            ! check total number of iteration
552            IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
553            ELSE                 ; nit = nit + 1
554            END IF
555
556            ! compute error between 2 iterations
557            ! if needed save gammat and compute zhtflx_b for next iteration
558            zerr = MAXVAL(ABS(zhtflx-zhtflx_b))
559            IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE.
560            ELSE                        ; zhtflx_b(:,:) = zhtflx(:,:)
561            END IF
562         END IF
563      END DO
564      !
565      CALL iom_put('isfgammat', zgammat)
566      CALL iom_put('isfgammas', zgammas)
567      !
568   END SUBROUTINE sbc_isf_cav
569
570
571   SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf )
572      !!----------------------------------------------------------------------
573      !! ** Purpose    : compute the coefficient echange for heat flux
574      !!
575      !! ** Method     : gamma assume constant or depends of u* and stability
576      !!
577      !! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
578      !!                Jenkins et al., 2010, JPO, p2298-2312
579      !!---------------------------------------------------------------------
580      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pgt   , pgs      !
581      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pqhisf, pqwisf   !
582      !
583      INTEGER  :: ji, jj                     ! loop index
584      INTEGER  :: ikt                        ! local integer
585      REAL(wp) :: zdku, zdkv                 ! U, V shear
586      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number
587      REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
588      REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
589      REAL(wp) :: zhmax                      ! limitation of mol
590      REAL(wp) :: zetastar                   ! stability parameter
591      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence
592      REAL(wp) :: zcoef                      ! temporary coef
593      REAL(wp) :: zdep
594      REAL(wp) :: zeps = 1.0e-20_wp   
595      REAL(wp), PARAMETER :: zxsiN = 0.052_wp   ! dimensionless constant
596      REAL(wp), PARAMETER :: znu   = 1.95e-6_wp ! kinamatic viscosity of sea water (m2.s-1)
597      REAL(wp), DIMENSION(2) :: zts, zab
598      REAL(wp), DIMENSION(jpi,jpj) :: zustar   ! U, V at T point and friction velocity
599      !!---------------------------------------------------------------------
600      !
601      SELECT CASE ( nn_gammablk )
602      CASE ( 0 ) ! gamma is constant (specified in namelist)
603         !! ISOMIP formulation (Hunter et al, 2006)
604         pgt(:,:) = rn_gammat0
605         pgs(:,:) = rn_gammas0
606
607      CASE ( 1 ) ! gamma is assume to be proportional to u*
608         !! Jenkins et al., 2010, JPO, p2298-2312
609         !! Adopted by Asay-Davis et al. (2015)
610         !! compute ustar (eq. 24)
611!!gm  NB  use pCdU here so that it will incorporate local boost of Cd0 and log layer case :
612!!         zustar(:,:) = SQRT( rCdU_top(:,:) * SQRT(utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
613!! or better :  compute ustar in zdfdrg  and use it here as well as in TKE, GLS and Co
614!!
615!!     ===>>>>    GM  to be done this chrismas
616!!         
617!!gm end
618         zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
619
620         !! Compute gammats
621         pgt(:,:) = zustar(:,:) * rn_gammat0
622         pgs(:,:) = zustar(:,:) * rn_gammas0
623     
624      CASE ( 2 ) ! gamma depends of stability of boundary layer
625         !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
626         !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
627         !! compute ustar
628         zustar(:,:) = SQRT( r_Cdmin_top * (utbl(:,:) * utbl(:,:) + vtbl(:,:) * vtbl(:,:) + r_ke0_top) )
629
630         !! compute Pr and Sc number (can be improved)
631         zPr =   13.8_wp
632         zSc = 2432.0_wp
633
634         !! compute gamma mole
635         zgmolet = 12.5_wp * zPr ** (2.0/3.0) - 6.0_wp
636         zgmoles = 12.5_wp * zSc ** (2.0/3.0) - 6.0_wp
637
638         !! compute gamma
639         DO ji = 2, jpi
640            DO jj = 2, jpj
641               ikt = mikt(ji,jj)
642
643               IF( zustar(ji,jj) == 0._wp ) THEN           ! only for kt = 1 I think
644                  pgt = rn_gammat0
645                  pgs = rn_gammas0
646               ELSE
647                  !! compute Rc number (as done in zdfric.F90)
648!!gm better to do it like in the new zdfric.F90   i.e. avm weighted Ri computation
649!!gm moreover, use Max(rn2,0) to take care of static instabilities....
650                  zcoef = 0.5_wp / e3w_n(ji,jj,ikt)
651                  !                                            ! shear of horizontal velocity
652                  zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )  &
653                     &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
654                  zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )  &
655                     &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
656                  !                                            ! richardson number (minimum value set to zero)
657                  zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps )
658
659                  !! compute bouyancy
660                  zts(jp_tem) = ttbl(ji,jj)
661                  zts(jp_sal) = stbl(ji,jj)
662                  zdep        = gdepw_n(ji,jj,ikt)
663                  !
664                  CALL eos_rab( zts, zdep, zab )
665                  !
666                  !! compute length scale
667                  zbuofdep = grav * ( zab(jp_tem) * pqhisf(ji,jj) - zab(jp_sal) * pqwisf(ji,jj) )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
668
669                  !! compute Monin Obukov Length
670                  ! Maximum boundary layer depth
671                  zhmax = gdept_n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp
672                  ! Compute Monin obukhov length scale at the surface and Ekman depth:
673                  zmob   = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps))
674                  zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
675
676                  !! compute eta* (stability parameter)
677                  zetastar = 1._wp / ( SQRT(1._wp + MAX(zxsiN * zustar(ji,jj) / ( ABS(ff_f(ji,jj)) * zmols * zRc ), 0._wp)))
678
679                  !! compute the sublayer thickness
680                  zhnu = 5 * znu / zustar(ji,jj)
681
682                  !! compute gamma turb
683                  zgturb = 1._wp / vkarmn * LOG(zustar(ji,jj) * zxsiN * zetastar * zetastar / ( ABS(ff_f(ji,jj)) * zhnu )) &
684                  &      + 1._wp / ( 2 * zxsiN * zetastar ) - 1._wp / vkarmn
685
686                  !! compute gammats
687                  pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet)
688                  pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles)
689               END IF
690            END DO
691         END DO
692         CALL lbc_lnk_multi( pgt, 'T', 1., pgs, 'T', 1.)
693      END SELECT
694      !
695   END SUBROUTINE sbc_isf_gammats
696
697
698   SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin )
699      !!----------------------------------------------------------------------
700      !!                  ***  SUBROUTINE sbc_isf_tbl  ***
701      !!
702      !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
703      !!
704      !!----------------------------------------------------------------------
705      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) :: pvarin
706      REAL(wp), DIMENSION(:,:)  , INTENT(  out) :: pvarout
707      CHARACTER(len=1),           INTENT(in   ) :: cd_ptin ! point of variable in/out
708      !
709      INTEGER ::   ji, jj, jk                ! loop index
710      INTEGER ::   ikt, ikb                    ! top and bottom index of the tbl
711      REAL(wp) ::   ze3, zhk
712      REAL(wp), DIMENSION(jpi,jpj) :: zhisf_tbl ! thickness of the tbl
713      !!----------------------------------------------------------------------
714     
715      ! initialisation
716      pvarout(:,:)=0._wp
717   
718      SELECT CASE ( cd_ptin )
719      CASE ( 'U' ) ! compute U in the top boundary layer at T- point
720         DO jj = 1,jpj
721            DO ji = 1,jpi
722               ikt = miku(ji,jj) ; ikb = miku(ji,jj)
723               ! thickness of boundary layer at least the top level thickness
724               zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u_n(ji,jj,ikt) )
725
726               ! determine the deepest level influenced by the boundary layer
727               DO jk = ikt+1, mbku(ji,jj)
728                  IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk
729               END DO
730               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
731
732               ! level fully include in the ice shelf boundary layer
733               DO jk = ikt, ikb - 1
734                  ze3 = e3u_n(ji,jj,jk)
735                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
736               END DO
737
738               ! level partially include in ice shelf boundary layer
739               zhk = SUM( e3u_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
740               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
741            END DO
742         END DO
743         DO jj = 2, jpj
744            DO ji = 2, jpi
745!!gm a wet-point only average should be used here !!!
746               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji-1,jj))
747            END DO
748         END DO
749         CALL lbc_lnk(pvarout,'T',-1.)
750     
751      CASE ( 'V' ) ! compute V in the top boundary layer at T- point
752         DO jj = 1,jpj
753            DO ji = 1,jpi
754               ikt = mikv(ji,jj) ; ikb = mikv(ji,jj)
755               ! thickness of boundary layer at least the top level thickness
756               zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v_n(ji,jj,ikt))
757
758               ! determine the deepest level influenced by the boundary layer
759               DO jk = ikt+1, mbkv(ji,jj)
760                  IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk
761               END DO
762               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
763
764               ! level fully include in the ice shelf boundary layer
765               DO jk = ikt, ikb - 1
766                  ze3 = e3v_n(ji,jj,jk)
767                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3
768               END DO
769
770               ! level partially include in ice shelf boundary layer
771               zhk = SUM( e3v_n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)
772               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
773            END DO
774         END DO
775         DO jj = 2, jpj
776            DO ji = 2, jpi
777!!gm a wet-point only average should be used here !!!
778               pvarout(ji,jj) = 0.5_wp * (pvarout(ji,jj) + pvarout(ji,jj-1))
779            END DO
780         END DO
781         CALL lbc_lnk(pvarout,'T',-1.)
782
783      CASE ( 'T' ) ! compute T in the top boundary layer at T- point
784         DO jj = 1,jpj
785            DO ji = 1,jpi
786               ikt = misfkt(ji,jj)
787               ikb = misfkb(ji,jj)
788
789               ! level fully include in the ice shelf boundary layer
790               DO jk = ikt, ikb - 1
791                  ze3 = e3t_n(ji,jj,jk)
792                  pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
793               END DO
794
795               ! level partially include in ice shelf boundary layer
796               zhk = SUM( e3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
797               pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk)
798            END DO
799         END DO
800      END SELECT
801      !
802      ! mask mean tbl value
803      pvarout(:,:) = pvarout(:,:) * ssmask(:,:)
804      !
805   END SUBROUTINE sbc_isf_tbl
806     
807
808   SUBROUTINE sbc_isf_div( phdivn )
809      !!----------------------------------------------------------------------
810      !!                  ***  SUBROUTINE sbc_isf_div  ***
811      !!       
812      !! ** Purpose :   update the horizontal divergence with the runoff inflow
813      !!
814      !! ** Method  :   
815      !!                CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the
816      !!                          divergence and expressed in m/s
817      !!
818      !! ** Action  :   phdivn   decreased by the runoff inflow
819      !!----------------------------------------------------------------------
820      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   phdivn   ! horizontal divergence
821      !
822      INTEGER  ::   ji, jj, jk   ! dummy loop indices
823      INTEGER  ::   ikt, ikb 
824      REAL(wp) ::   zhk
825      REAL(wp) ::   zfact     ! local scalar
826      !!----------------------------------------------------------------------
827      !
828      zfact   = 0.5_wp
829      !
830      IF(.NOT.ln_linssh ) THEN     ! need to re compute level distribution of isf fresh water
831         DO jj = 1,jpj
832            DO ji = 1,jpi
833               ikt = misfkt(ji,jj)
834               ikb = misfkt(ji,jj)
835               ! thickness of boundary layer at least the top level thickness
836               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t_n(ji,jj,ikt))
837
838               ! determine the deepest level influenced by the boundary layer
839               DO jk = ikt, mbkt(ji,jj)
840                  IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
841               END DO
842               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
843               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
844               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
845
846               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
847               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t_n(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
848            END DO
849         END DO
850      END IF 
851      !
852      !==   ice shelf melting distributed over several levels   ==!
853      DO jj = 1,jpj
854         DO ji = 1,jpi
855               ikt = misfkt(ji,jj)
856               ikb = misfkb(ji,jj)
857               ! level fully include in the ice shelf boundary layer
858               DO jk = ikt, ikb - 1
859                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
860                    &              * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
861               END DO
862               ! level partially include in ice shelf boundary layer
863               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
864                    &            + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
865         END DO
866      END DO
867      !
868   END SUBROUTINE sbc_isf_div
869
870   !!======================================================================
871END MODULE sbcisf
Note: See TracBrowser for help on using the repository browser.