source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 4726

Last change on this file since 4726 was 4726, checked in by mathiot, 6 years ago

ISF branch: change name of 2 variables (icedep ⇒ risfdep and lmask ⇒ ssmask), cosmetic changes and add ldfslp key

File size: 44.0 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
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 lbclnk          !
21   USE iom             ! I/O manager library
22   USE in_out_manager  ! I/O manager
23   USE wrk_nemo        ! Memory allocation
24   USE timing          ! Timing
25   USE lib_fortran     ! glob_sum
26   USE zdfbfr
27   USE fldread         ! read input field at current time step
28
29
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   sbc_isf, sbc_isf_div, sbc_isf_alloc  ! routine called in sbcmod and divcur
35
36   ! public in order to be able to output then
37
38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc   
39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fwfisf_b, fwfisf  !: evaporation damping   [kg/m2/s]
40   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf            !: net heat flux from ice shelf
41   REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m]
42   LOGICAL , PUBLIC ::   ln_divisf                   !: flag to correct divergence
43   INTEGER , PUBLIC ::   nn_isfblk                   !:
44   INTEGER , PUBLIC ::   nn_gammablk                 !:
45   LOGICAL , PUBLIC ::   ln_conserve                 !:
46   REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient
47   REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient
48   REAL(wp), PUBLIC ::   rdivisf                     !: flag to test if fwf apply on divergence
49
50   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rzisf_tbl              !:depth of ice shelf base  ????
51   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rhisf_tbl, rhisf_tbl_0 !:depth of ice shelf base  ????
52   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  r1_hisf_tbl            !:1/depth of ice shelf base  ????
53   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ralpha                 !:proportion of bottom cell influenced by boundary layer
54   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  risfLeff               !:effective length (Leff) BG03 nn_isf==2 ?
55   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point
56   INTEGER(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
57                                                                                          !: (first wet level and last level include in the tbl)
58
59   REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ?
60   REAL(wp), PUBLIC, SAVE ::   kappa  = 1.54e-6_wp    ! phycst ?
61   REAL(wp), PUBLIC, SAVE ::   rhoisf = 920.0_wp      ! phycst ?
62   REAL(wp), PUBLIC, SAVE ::   tsurf  = -20.0_wp      ! phycst ?
63   REAL(wp), PUBLIC, SAVE ::   lfusisf= 0.334e6_wp    ! phycst ?
64
65!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
66   CHARACTER(len=100), PUBLIC ::   cn_dirisf  = './'    !: Root directory for location of ssr files
67   TYPE(FLD_N)       , PUBLIC ::   sn_qisf, sn_fwfisf     !: information about the runoff file to be read
68   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qisf, sf_fwfisf
69   TYPE(FLD_N)       , PUBLIC ::   sn_rnfisf              !: information about the runoff file to be read
70   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnfisf           
71   TYPE(FLD_N)       , PUBLIC ::   sn_depmax_isf, sn_depmin_isf, sn_Leff_isf     !: information about the runoff file to be read
72   
73   !! * Substitutions
74#  include "domzgr_substitute.h90"
75   !!----------------------------------------------------------------------
76   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
77   !! $Id: sbcice_if.F90 1730 2009-11-16 14:34:19Z smasson $
78   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
79   !!----------------------------------------------------------------------
80
81CONTAINS
82 
83  SUBROUTINE sbc_isf(kt)
84    INTEGER, INTENT(in)          ::   kt         ! ocean time step
85    INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
86    INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
87    REAL(wp)                     ::   rmin
88    REAL(wp)                     ::   zhk
89    CHARACTER(len=256)           ::   cfisf, cvarzisf, cvarhisf   ! name for isf file
90    CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file
91    CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale
92    INTEGER           ::   ios           ! Local integer output status for namelist read
93      !
94      !!---------------------------------------------------------------------
95      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, ln_divisf, ln_conserve, rn_gammat0, rn_gammas0, nn_gammablk, &
96                         & sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
97      !
98      !
99      !                                         ! ====================== !
100      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
101         !                                      ! ====================== !
102         REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs
103         READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
104901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp )
105
106         REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
107         READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
108902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
109         IF(lwm) WRITE ( numond, namsbc_isf )
110
111
112         IF ( lwp ) WRITE(numout,*)
113         IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf'
114         IF ( lwp ) WRITE(numout,*) '~~~~~~~~~'
115         IF ( lwp ) WRITE(numout,*) 'sbcisf :' 
116         IF ( lwp ) WRITE(numout,*) '~~~~~~~~'
117         IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf
118         IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk
119         IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl
120         IF ( lwp ) WRITE(numout,*) '        ln_divisf   = ', ln_divisf 
121         IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk 
122         IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2 
123         IF (ln_divisf) THEN       ! keep it in the namelist ??? used true anyway as for runoff ? (PM)
124            rdivisf = 1._wp
125         ELSE
126            rdivisf = 0._wp
127         END IF
128         !
129         ! Allocate public variable
130         IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
131         !
132         ! initialisation
133         qisf(:,:)        = 0._wp  ; fwfisf(:,:) = 0._wp
134         risf_tsc(:,:,:)  = 0._wp
135         !
136         ! define isf tbl tickness, top and bottom indice
137         IF      (nn_isf == 1) THEN
138            rhisf_tbl(:,:) = rn_hisf_tbl
139            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
140         ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN
141            ALLOCATE( sf_rnfisf(1), STAT=ierror )
142            ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
143            CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
144
145            !: read effective lenght (BG03)
146            IF (nn_isf == 2) THEN
147               ! Read Data and save some integral values
148               CALL iom_open( sn_Leff_isf%clname, inum )
149               cvarLeff  = 'soLeff'               !: variable name for Efficient Length scale
150               CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
151               CALL iom_close(inum)
152               !
153               risfLeff = risfLeff*1000           !: convertion in m
154            END IF
155
156           ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
157            CALL iom_open( sn_depmax_isf%clname, inum )
158            cvarhisf = TRIM(sn_depmax_isf%clvar)
159            CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
160            CALL iom_close(inum)
161            !
162            CALL iom_open( sn_depmin_isf%clname, inum )
163            cvarzisf = TRIM(sn_depmin_isf%clvar)
164            CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
165            CALL iom_close(inum)
166            !
167            rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer
168
169           !! compute first level of the top boundary layer
170           DO ji = 1, jpi
171              DO jj = 1, jpj
172                  jk = 2
173                  DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdepw(ji,jj,jk) < rzisf_tbl(ji,jj) ) ;  jk = jk + 1 ;  END DO
174                  misfkt(ji,jj) = jk-1
175               END DO
176            END DO
177
178         ELSE IF ( nn_isf == 4 ) THEN
179            ! as in nn_isf == 1
180            rhisf_tbl(:,:) = rn_hisf_tbl
181            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
182           
183            ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
184            ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror )
185            ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
186            ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) )
187            CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
188            !CALL fld_fill( sf_qisf  , (/ sn_qisf   /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data'       , 'namsbc_isf' )
189         END IF
190         
191         rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
192
193         ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
194         DO jj = 1,jpj
195            DO ji = 1,jpi
196               ikt = misfkt(ji,jj)
197               ikb = misfkt(ji,jj)
198               ! thickness of boundary layer at least the top level thickness
199               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt))
200
201               ! determine the deepest level influenced by the boundary layer
202               ! test on tmask useless ?????
203               DO jk = ikt, mbkt(ji,jj)
204                  IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
205               END DO
206               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
207               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
208               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
209
210               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
211               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
212            END DO
213         END DO
214         
215      END IF
216
217      !                                            ! ---------------------------------------- !
218      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
219         !                                         ! ---------------------------------------- !
220         fwfisf_b  (:,:  ) = fwfisf  (:,:  )               ! Swap the ocean forcing fields except at nit000
221         risf_tsc_b(:,:,:) = risf_tsc(:,:,:)               ! where before fields are set at the end of the routine
222         !
223      ENDIF
224
225      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN
226
227
228         ! compute salf and heat flux
229         IF (nn_isf == 1) THEN
230            ! realistic ice shelf formulation
231            ! compute T/S/U/V for the top boundary layer
232            CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')
233            CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T')
234            CALL sbc_isf_tbl(un(:,:,:),utbl(:,:),'U')
235            CALL sbc_isf_tbl(vn(:,:,:),vtbl(:,:),'V')
236            ! iom print
237            CALL iom_put('ttbl',ttbl(:,:))
238            CALL iom_put('stbl',stbl(:,:))
239            CALL iom_put('utbl',utbl(:,:))
240            CALL iom_put('vtbl',vtbl(:,:))
241            ! compute fwf and heat flux
242            CALL sbc_isf_cav (kt)
243
244         ELSE IF (nn_isf == 2) THEN
245            ! Beckmann and Goosse parametrisation
246            stbl(:,:)   = soce
247            CALL sbc_isf_bg03(kt)
248
249         ELSE IF (nn_isf == 3) THEN
250            ! specified runoff in depth (Mathiot et al., XXXX in preparation)
251            CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
252            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting)
253            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
254            stbl(:,:)   = soce
255
256         ELSE IF (nn_isf == 4) THEN
257            ! specified fwf and heat flux forcing beneath the ice shelf
258            CALL fld_read ( kt, nn_fsbc, sf_fwfisf   )
259            !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
260            fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf
261            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
262            !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux
263            stbl(:,:)   = soce
264
265         END IF
266         ! compute tsc due to isf
267         ! WARNING water add at temp = 0C, correction term is added in trasbc, maybe better here but need a 3D variable).
268         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp !
269         
270         ! salt effect already take into account in vertical advection
271         risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0
272         
273         ! lbclnk
274         CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.)
275         CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.)
276         CALL lbc_lnk(fwfisf(:,:)   ,'T',1.)
277         CALL lbc_lnk(qisf(:,:)     ,'T',1.)
278
279         IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
280            IF( ln_rstart .AND.    &                     ! Restart: read in restart file
281                 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN
282               IF(lwp) WRITE(numout,*) '          nit000-1 isf tracer content forcing fields read in the restart file'
283               CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:) )   ! before salt content isf_tsc trend
284               CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b', risf_tsc_b(:,:,jp_sal) )   ! before salt content isf_tsc trend
285               CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b', risf_tsc_b(:,:,jp_tem) )   ! before salt content isf_tsc trend
286            ELSE
287               fwfisf_b(:,:)    = fwfisf(:,:)
288               risf_tsc_b(:,:,:)= risf_tsc(:,:,:)
289            END IF
290         ENDIF
291         !
292         ! output
293         CALL iom_put('qisf'  , qisf)
294         CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce )
295      END IF
296 
297  END SUBROUTINE sbc_isf
298
299  INTEGER FUNCTION sbc_isf_alloc()
300      !!----------------------------------------------------------------------
301      !!               ***  FUNCTION sbc_isf_rnf_alloc  ***
302      !!----------------------------------------------------------------------
303      sbc_isf_alloc = 0       ! set to zero if no array to be allocated
304      IF( .NOT. ALLOCATED( qisf ) ) THEN
305         ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj), fwfisf(jpi,jpj), &
306               &    fwfisf_b(jpi,jpj), misfkt(jpi,jpj), rhisf_tbl(jpi,jpj), r1_hisf_tbl(jpi,jpj),     &
307               &    rzisf_tbl(jpi,jpj), misfkb(jpi,jpj), ttbl(jpi,jpj), stbl(jpi,jpj), utbl(jpi,jpj), &
308               &    vtbl(jpi, jpj), risfLeff(jpi,jpj), rhisf_tbl_0(jpi,jpj), ralpha(jpi,jpj), STAT= sbc_isf_alloc )
309         !
310         IF( lk_mpp                  )   CALL mpp_sum ( sbc_isf_alloc )
311         IF( sbc_isf_alloc /= 0 )   CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.')
312         !
313      ENDIF
314  END FUNCTION
315
316  SUBROUTINE sbc_isf_bg03(kt)
317   !!==========================================================================
318   !!                 *** SUBROUTINE sbcisf_bg03  ***
319   !! add net heat and fresh water flux from ice shelf melting
320   !! into the adjacent ocean using the parameterisation by
321   !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
322   !!     interaction for climate models", Ocean Modelling 5(2003) 157-170.
323   !!  (hereafter BG)
324   !!==========================================================================
325   !!----------------------------------------------------------------------
326   !!   sbc_isf_bg03      : routine called from sbcmod
327   !!----------------------------------------------------------------------
328   !!
329   !! ** Purpose   :   Add heat and fresh water fluxes due to ice shelf melting
330   !! ** Reference :   Beckmann et Goosse, 2003, Ocean Modelling
331   !!
332   !! History :
333   !!      !  06-02  (C. Wang) Original code
334   !!----------------------------------------------------------------------
335
336    INTEGER, INTENT ( in ) :: kt
337
338    INTEGER :: ji, jj, jk, jish  !temporary integer
339    INTEGER :: ijkmin
340    INTEGER :: ii, ij, ik 
341    INTEGER :: inum
342
343    REAL(wp) :: zt_sum      ! sum of the temperature between 200m and 600m
344    REAL(wp) :: zt_ave      ! averaged temperature between 200m and 600m
345    REAL(wp) :: zt_frz      ! freezing point temperature at depth z
346    REAL(wp) :: zpress      ! pressure to compute the freezing point in depth
347   
348    !!----------------------------------------------------------------------
349    IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03')
350     !
351
352    ! This test is false only in the very first time step of a run (JMM ???- Initialy build to skip 1rst year of run )
353    DO ji = 1, jpi
354       DO jj = 1, jpj
355          ik = misfkt(ji,jj)
356          !! Initialize arrays to 0 (each step)
357          zt_sum = 0.e0_wp
358          IF ( ik .GT. 1 ) THEN
359    ! 3. -----------the average temperature between 200m and 600m ---------------------
360             DO jk = misfkt(ji,jj),misfkb(ji,jj)
361             ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!)
362             ! after verif with UNESCO, wrong sign in BG eq. 2
363             ! Calculate freezing temperature
364                zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04 
365                zt_frz = tfreez1D(tsb(ji,jj,ik,jp_sal), zpress) 
366                zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik)  ! sum temp
367             ENDDO
368             zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
369   
370    ! 4. ------------Net heat flux and fresh water flux due to the ice shelf
371          ! For those corresponding to zonal boundary   
372             qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
373                         & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,ik) 
374             
375             fwfisf(ji,jj) = qisf(ji,jj) / lfusisf          !fresh water flux kg/(m2s)                 
376             fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
377             !add to salinity trend
378          ELSE
379             qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp
380          END IF
381       ENDDO
382    ENDDO
383    !
384    IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03')
385  END SUBROUTINE sbc_isf_bg03
386
387   SUBROUTINE sbc_isf_cav( kt )
388      !!---------------------------------------------------------------------
389      !!                     ***  ROUTINE sbc_isf_cav  ***
390      !!
391      !! ** Purpose :   handle surface boundary condition under ice shelf
392      !!
393      !! ** Method  : -
394      !!
395      !! ** Action  :   utau, vtau : remain unchanged
396      !!                taum, wndm : remain unchanged
397      !!                qns        : update heat flux below ice shelf
398      !!                emp, emps  : update freshwater flux below ice shelf
399      !!---------------------------------------------------------------------
400      INTEGER, INTENT(in)          ::   kt         ! ocean time step
401      !
402      LOGICAL :: ln_isomip = .true.
403      REAL(wp), DIMENSION(:,:), POINTER       ::   zfrz,zpress,zti
404      REAL(wp), DIMENSION(:,:), POINTER, SAVE ::   zgammat2d, zgammas2d 
405      !REAL(wp), DIMENSION(:,:), POINTER ::   zqisf, zfwfisf
406      REAL(wp) ::   zlamb1, zlamb2, zlamb3
407      REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7
408      REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac
409      REAL(wp) ::   zfwflx, zhtflx, zhtflx_b
410      REAL(wp) ::   zgammat, zgammas
411      REAL(wp) ::   zeps   =  -1.e-20_wp        !==   Local constant initialization   ==!
412      INTEGER  ::   ji, jj     ! dummy loop indices
413      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
414      INTEGER  ::   ierror     ! return error code
415      LOGICAL  ::   lit=.TRUE.
416      INTEGER  ::   nit
417      !!---------------------------------------------------------------------
418      !
419      ! coeficient for linearisation of tfreez
420      zlamb1=-0.0575
421      zlamb2=0.0901
422      zlamb3=-7.61e-04
423      IF( nn_timing == 1 )  CALL timing_start('sbc_isf_cav')
424      !
425      CALL wrk_alloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
426
427      zcfac=0.0_wp 
428      IF (ln_conserve)  zcfac=1.0_wp
429      zpress(:,:)=0.0_wp
430      zgammat2d(:,:)=0.0_wp
431      zgammas2d(:,:)=0.0_wp
432      !
433      !
434!CDIR COLLAPSE
435      DO jj = 1, jpj
436         DO ji = 1, jpi
437            ! Crude approximation for pressure (but commonly used)
438            ! 1e-04 to convert from Pa to dBar
439            zpress(ji,jj)=grav*rau0*fsdepw(ji,jj,mikt(ji,jj))*1.e-04
440            !
441         END DO
442      END DO
443
444! Calculate in-situ temperature (ref to surface)
445      zti(:,:)=tinsitu( ttbl, stbl, zpress )
446! Calculate freezing temperature
447      zfrz(:,:)=tfreez( sss_m(:,:), zpress )
448
449     
450      zhtflx=0._wp ; zfwflx=0._wp
451      IF (nn_isfblk == 1) THEN
452         DO jj = 1, jpj
453            DO ji = 1, jpi
454               IF (mikt(ji,jj) > 1 ) THEN
455                  nit = 1; lit = .TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp
456                  DO WHILE ( lit )
457! compute gamma
458                     CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit)
459! zhtflx is upward heat flux (out of ocean)
460                     zhtflx = zgammat*rcp*rau0*(zti(ji,jj)-zfrz(ji,jj))
461! zwflx is upward water flux
462                     zfwflx = - zhtflx/lfusisf
463! test convergence and compute gammat
464                     IF ( (zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE.
465
466                     nit = nit + 1
467                     IF (nit .GE. 100) THEN
468                        !WRITE(numout,*) "sbcisf : too many iteration ... ", zhtflx, zhtflx_b,zgammat, rn_gammat0, rn_tfri2, nn_gammablk, ji,jj
469                        !WRITE(numout,*) "sbcisf : too many iteration ... ", (zhtflx - zhtflx_b)/zhtflx
470                        CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
471                     END IF
472! save gammat and compute zhtflx_b
473                     zgammat2d(ji,jj)=zgammat
474                     zhtflx_b = zhtflx
475                  END DO
476
477                  qisf(ji,jj) = - zhtflx
478! For genuine ISOMIP protocol this should probably be something like
479                  fwfisf(ji,jj) = zfwflx  * ( soce / MAX(stbl(ji,jj),zeps))
480               ELSE
481                  fwfisf(ji,jj) = 0._wp
482                  qisf(ji,jj)   = 0._wp
483               END IF
484            !
485            END DO
486         END DO
487
488      ELSE IF (nn_isfblk == 2 ) THEN
489
490! More complicated 3 equation thermodynamics as in MITgcm
491!CDIR COLLAPSE
492         DO jj = 2, jpj
493            DO ji = 2, jpi
494               IF (mikt(ji,jj) > 1 ) THEN
495                  nit=1; lit=.TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp; zhtflx=0._wp
496                  DO WHILE ( lit )
497                     CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit)
498
499                     zeps1=rcp*rau0*zgammat
500                     zeps2=lfusisf*rau0*zgammas
501                     zeps3=rhoisf*rcpi*kappa/risfdep(ji,jj)
502                     zeps4=zlamb2+zlamb3*risfdep(ji,jj)
503                     zeps6=zeps4-zti(ji,jj)
504                     zeps7=zeps4-tsurf
505                     zaqe=zlamb1 * (zeps1 + zeps3)
506                     zaqer=0.5/zaqe
507                     zbqe=zeps1*zeps6+zeps3*zeps7-zeps2
508                     zcqe=zeps2*stbl(ji,jj)
509                     zdis=zbqe*zbqe-4.0*zaqe*zcqe               
510! Presumably zdis can never be negative because gammas is very small compared to gammat
511                     zsfrz=(-zbqe-SQRT(zdis))*zaqer
512                     IF (zsfrz .lt. 0.0) zsfrz=(-zbqe+SQRT(zdis))*zaqer
513                     zfrz(ji,jj)=zeps4+zlamb1*zsfrz
514 
515! zfwflx is upward water flux
516                     zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / zsfrz )
517! zhtflx is upward heat flux (out of ocean)
518! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
519                     zhtflx = ( zgammat*rau0 - zcfac*zfwflx ) * rcp * (zti(ji,jj) - zfrz(ji,jj) ) 
520! zwflx is upward water flux
521! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
522                     zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj)
523! test convergence and compute gammat
524                     IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE.
525
526                     nit = nit + 1
527                     IF (nit .GE. 51) THEN
528                        WRITE(numout,*) "sbcisf : too many iteration ... ", zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea
529                        CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
530                     END IF
531! save gammat and compute zhtflx_b
532                     zgammat2d(ji,jj)=zgammat
533                     zgammas2d(ji,jj)=zgammas
534                     zhtflx_b = zhtflx
535
536                  END DO
537! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
538                  qisf(ji,jj) = - zhtflx 
539! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
540                  fwfisf(ji,jj) = zfwflx 
541               ELSE
542                  fwfisf(ji,jj) = 0._wp
543                  qisf(ji,jj)   = 0._wp
544               ENDIF
545               !
546            END DO
547         END DO
548      ENDIF
549      ! lbclnk
550      CALL lbc_lnk(zgammas2d(:,:),'T',1.)
551      CALL lbc_lnk(zgammat2d(:,:),'T',1.)
552      ! output
553      CALL iom_put('isfgammat', zgammat2d)
554      CALL iom_put('isfgammas', zgammas2d)
555         !
556      !CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zqisf, zfwfisf  )
557      CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
558      !
559      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav')
560
561   END SUBROUTINE sbc_isf_cav
562
563   SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit )
564      !!----------------------------------------------------------------------
565      !! ** Purpose    : compute the coefficient echange for heat flux
566      !!
567      !! ** Method     : gamma assume constant or depends of u* and stability
568      !!
569      !! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
570      !!                Jenkins et al., 2010, JPO, p2298-2312
571      !!---------------------------------------------------------------------
572      REAL(wp), INTENT(inout) :: gt, gs, zqhisf, zqwisf
573      INTEGER , INTENT(in)    :: ji,jj
574      LOGICAL , INTENT(inout) :: lit
575
576      INTEGER  :: ikt                 ! loop index
577      REAL(wp) :: zut, zvt, zustar           ! U, V at T point and friction velocity
578      REAL(wp) :: zdku, zdkv                 ! U, V shear
579      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number
580      REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
581      REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
582      REAL(wp) :: zhmax                      ! limitation of mol
583      REAL(wp) :: zetastar                   ! stability parameter
584      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence
585      REAL(wp) :: zcoef                      ! temporary coef
586      REAL(wp) :: zrhos, zalbet, zbeta, zthermal, zhalin
587      REAL(wp) :: zt, zs, zh
588      REAL(wp), PARAMETER :: zxsiN = 0.052   ! dimensionless constant
589      REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number
590      REAL(wp), PARAMETER :: znu   = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1)
591      REAL(wp) ::   rcs      = 1.0e-3_wp        ! conversion: mm/s ==> m/s
592      !!---------------------------------------------------------------------
593      !
594      IF( nn_gammablk == 0 ) THEN
595      !! gamma is constant (specified in namelist)
596         gt = rn_gammat0
597         gs = rn_gammas0
598         lit = .FALSE.
599      ELSE IF ( nn_gammablk == 1 ) THEN
600      !! gamma is assume to be proportional to u*
601      !! WARNING in case of Losh 2008 tbl parametrization,
602      !! you have to used the mean value of u in the boundary layer)
603      !! not yet coded
604      !! Jenkins et al., 2010, JPO, p2298-2312
605         ikt = mikt(ji,jj)
606      !! Compute U and V at T points
607   !      zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
608   !      zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
609          zut = utbl(ji,jj)
610          zvt = vtbl(ji,jj)
611
612      !! compute ustar
613         zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
614      !! Compute mean value over the TBL
615
616      !! Compute gammats
617         gt = zustar * rn_gammat0
618         gs = zustar * rn_gammas0
619         lit = .FALSE.
620      ELSE IF ( nn_gammablk == 2 ) THEN
621      !! gamma depends of stability of boundary layer
622      !! WARNING in case of Losh 2008 tbl parametrization,
623      !! you have to used the mean value of u in the boundary layer)
624      !! not yet coded
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               ikt = mikt(ji,jj)
628
629      !! Compute U and V at T points
630               zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
631               zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
632
633      !! compute ustar
634               zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
635               IF (zustar == 0._wp) THEN           ! only for kt = 1 I think
636                 gt = rn_gammat0
637                 gs = rn_gammas0
638               ELSE
639      !! compute Rc number (as done in zdfric.F90)
640               zcoef = 0.5 / fse3w(ji,jj,ikt)
641               !                                            ! shear of horizontal velocity
642               zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )   &
643                  &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
644               zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )   &
645                  &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
646               !                                            ! richardson number (minimum value set to zero)
647               zRc = rn2(ji,jj,ikt+1) / ( zdku*zdku + zdkv*zdkv + 1.e-20 )
648
649      !! compute Pr and Sc number (can be improved)
650               zPr =   13.8
651               zSc = 2432.0
652
653      !! compute gamma mole
654               zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0
655               zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0
656
657      !! compute bouyancy
658               IF( nn_eos < 1) THEN
659                  zt     = ttbl(ji,jj)
660                  zs     = stbl(ji,jj) - 35.0
661                  zh     = fsdepw(ji,jj,ikt)
662                  !  potential volumic mass
663                  zrhos  = rhop(ji,jj,ikt)
664                  zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt   &   ! ratio alpha/beta
665                     &                               - 0.203814e-03 ) * zt   &
666                     &                               + 0.170907e-01 ) * zt   &
667                     &   + 0.665157e-01                                      &
668                     &   +     ( - 0.678662e-05 * zs                         &
669                     &           - 0.846960e-04 * zt + 0.378110e-02 ) * zs   &
670                     &   +   ( ( - 0.302285e-13 * zh                         &
671                     &           - 0.251520e-11 * zs                         &
672                     &           + 0.512857e-12 * zt * zt           ) * zh   &
673                     &           - 0.164759e-06 * zs                         &
674                     &        +(   0.791325e-08 * zt - 0.933746e-06 ) * zt   &
675                     &                               + 0.380374e-04 ) * zh
676
677                  zbeta  = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt      &   ! beta
678                     &                            - 0.301985e-05 ) * zt      &
679                     &   + 0.785567e-03                                      &
680                     &   + (     0.515032e-08 * zs                           &
681                     &         + 0.788212e-08 * zt - 0.356603e-06 ) * zs     &
682                     &   +(  (   0.121551e-17 * zh                           &
683                     &         - 0.602281e-15 * zs                           &
684                     &         - 0.175379e-14 * zt + 0.176621e-12 ) * zh     &
685                     &                             + 0.408195e-10   * zs     &
686                     &     + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt     &
687                     &                             - 0.121555e-07 ) * zh
688
689                  zthermal = zbeta * zalbet / ( rcp * zrhos + epsln )
690                  zhalin   = zbeta * stbl(ji,jj) * rcs
691               ELSE
692                  zrhos    = rhop(ji,jj,ikt) + rau0 * ( 1. - tmask(ji,jj,ikt) )
693                  zthermal = rn_alpha / ( rcp * zrhos + epsln )
694                  zhalin   = rn_beta * stbl(ji,jj) * rcs
695               ENDIF
696      !! compute length scale
697               zbuofdep = grav * ( zthermal * zqhisf - zhalin * zqwisf )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
698
699      !! compute Monin Obukov Length
700               ! Maximum boundary layer depth
701               zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) -0.001
702               ! Compute Monin obukhov length scale at the surface and Ekman depth:
703               zmob   = zustar ** 3 / (vkarmn * (zbuofdep + epsln))
704               zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
705
706      !! compute eta* (stability parameter)
707               zetastar = 1 / ( SQRT(1 + MAX(zxsiN * zustar / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0)))
708
709      !! compute the sublayer thickness
710               zhnu = 5 * znu / zustar
711      !! compute gamma turb
712               zgturb = 1/vkarmn * LOG(zustar * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &
713               &      + 1 / ( 2 * zxsiN * zetastar ) - 1/vkarmn
714
715      !! compute gammats
716               gt = zustar / (zgturb + zgmolet)
717               gs = zustar / (zgturb + zgmoles)
718               END IF
719      END IF
720
721   END SUBROUTINE
722
723   SUBROUTINE sbc_isf_tbl( varin, varout, cptin )
724      !!----------------------------------------------------------------------
725      !!                  ***  SUBROUTINE sbc_isf_tbl  ***
726      !!
727      !! ** Purpose : compute mean T/S/U/V in the boundary layer
728      !!
729      !!----------------------------------------------------------------------
730      REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin
731      REAL(wp), DIMENSION(:,:)  , INTENT(out):: varout
732     
733      CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out
734
735      REAL(wp) :: ze3, zhk
736      REAL(wp), DIMENSION(:,:), POINTER :: zikt
737
738      INTEGER :: ji,jj,jk
739      INTEGER :: ikt,ikb
740      INTEGER, DIMENSION(:,:), POINTER :: mkt, mkb
741
742      CALL wrk_alloc( jpi,jpj, mkt, mkb  )
743      CALL wrk_alloc( jpi,jpj, zikt )
744
745      ! get first and last level of tbl
746      mkt(:,:) = misfkt(:,:)
747      mkb(:,:) = misfkb(:,:)
748
749      varout(:,:)=0._wp
750      DO jj = 2,jpj
751         DO ji = 2,jpi
752            IF (ssmask(ji,jj) == 1) THEN
753               ikt = mkt(ji,jj)
754               ikb = mkb(ji,jj)
755
756               ! level fully include in the ice shelf boundary layer
757               DO jk = ikt, ikb - 1
758                  ze3 = fse3t_n(ji,jj,jk)
759                  IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
760                  IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) &
761                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
762                  IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) &
763                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
764               END DO
765
766               ! level partially include in ice shelf boundary layer
767               zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
768               IF (cptin == 'T') varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
769               IF (cptin == 'U') varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk)
770               IF (cptin == 'V') varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk)
771            END IF
772         END DO
773      END DO
774
775      CALL wrk_dealloc( jpi,jpj, mkt, mkb )     
776      CALL wrk_dealloc( jpi,jpj, zikt ) 
777
778      IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.)
779      IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.)
780
781   END SUBROUTINE sbc_isf_tbl
782     
783
784   SUBROUTINE sbc_isf_div( phdivn )
785      !!----------------------------------------------------------------------
786      !!                  ***  SUBROUTINE sbc_isf_div  ***
787      !!       
788      !! ** Purpose :   update the horizontal divergence with the runoff inflow
789      !!
790      !! ** Method  :   
791      !!                CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the
792      !!                          divergence and expressed in m/s
793      !!
794      !! ** Action  :   phdivn   decreased by the runoff inflow
795      !!----------------------------------------------------------------------
796      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
797      !!
798      INTEGER(wp)  ::   ji, jj, jk   ! dummy loop indices
799      INTEGER(wp)  ::   ikt, ikb 
800      INTEGER(wp)  ::   nk_isf
801      REAL(wp)     ::   zhk, z1_hisf_tbl, zhisf_tbl
802      REAL(wp)     ::   zfact     ! local scalar
803      !!----------------------------------------------------------------------
804      !
805      zfact   = 0.5_wp
806      !
807      IF (lk_vvl) THEN     ! need to re compute level distribution of isf fresh water
808         DO jj = 1,jpj
809            DO ji = 1,jpi
810               ikt = misfkt(ji,jj)
811               ikb = misfkt(ji,jj)
812               ! thickness of boundary layer at least the top level thickness
813               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t(ji,jj,ikt))
814
815               ! determine the deepest level influenced by the boundary layer
816               ! test on tmask useless ?????
817               DO jk = ikt, mbkt(ji,jj)
818                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
819               END DO
820               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
821               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
822               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
823
824               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
825               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
826            END DO
827         END DO
828      END IF ! vvl case
829      !
830      DO jj = 1,jpj
831         DO ji = 1,jpi
832               ikt = misfkt(ji,jj)
833               ikb = misfkb(ji,jj)
834               ! level fully include in the ice shelf boundary layer
835               DO jk = ikt, ikb - 1
836                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
837               END DO
838               ! level partially include in ice shelf boundary layer
839               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
840            !==   ice shelf melting mass distributed over several levels   ==!
841         END DO
842      END DO
843      !
844   END SUBROUTINE sbc_isf_div
845                       
846   FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti )
847      !!----------------------------------------------------------------------
848      !!                 ***  ROUTINE eos_init  ***
849      !!
850      !! ** Purpose :   Compute the in-situ temperature [Celcius]
851      !!
852      !! ** Method  :   
853      !!
854      !! Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
855      !!----------------------------------------------------------------------
856      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptem   ! potential temperature [Celcius]
857      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu]
858      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ppress ! pressure             [dBar]
859      REAL(wp), DIMENSION(:,:), POINTER           ::   pti    ! in-situ temperature [Celcius]
860!      REAL(wp) :: fsatg
861!      REAL(wp) :: pfps, pfpt, pfphp
862      REAL(wp) :: zt, zs, zp, zh, zq, zxk
863      INTEGER  :: ji, jj            ! dummy loop indices
864      !
865      CALL wrk_alloc( jpi,jpj, pti  )
866      !
867      DO jj=1,jpj
868         DO ji=1,jpi
869            zh = ppress(ji,jj)
870! Theta1
871            zt = ptem(ji,jj)
872            zs = psal(ji,jj)
873            zp = 0.0
874            zxk= zh * fsatg( zs, zt, zp )
875            zt = zt + 0.5 * zxk
876            zq = zxk
877! Theta2
878            zp = zp + 0.5 * zh
879            zxk= zh*fsatg( zs, zt, zp )
880            zt = zt + 0.29289322 * ( zxk - zq )
881            zq = 0.58578644 * zxk + 0.121320344 * zq
882! Theta3
883            zxk= zh * fsatg( zs, zt, zp )
884            zt = zt + 1.707106781 * ( zxk - zq )
885            zq = 3.414213562 * zxk - 4.121320344 * zq
886! Theta4
887            zp = zp + 0.5 * zh
888            zxk= zh * fsatg( zs, zt, zp )
889            pti(ji,jj) = zt + ( zxk - 2.0 * zq ) / 6.0
890         END DO
891      END DO
892      !
893      CALL wrk_dealloc( jpi,jpj, pti  )
894      !
895   END FUNCTION tinsitu
896   !
897   FUNCTION fsatg( pfps, pfpt, pfphp )
898      !!----------------------------------------------------------------------
899      !!                 ***  FUNCTION fsatg  ***
900      !!
901      !! ** Purpose    :   Compute the Adiabatic laspse rate [Celcius].[decibar]^-1
902      !!
903      !! ** Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
904      !!
905      !! ** units      :   pressure        pfphp    decibars
906      !!                   temperature     pfpt     deg celsius (ipts-68)
907      !!                   salinity        pfps     (ipss-78)
908      !!                   adiabatic       fsatg    deg. c/decibar
909      !!----------------------------------------------------------------------
910      REAL(wp) :: pfps, pfpt, pfphp 
911      REAL(wp) :: fsatg
912      !
913      fsatg = (((-2.1687e-16*pfpt+1.8676e-14)*pfpt-4.6206e-13)*pfphp         &
914        &    +((2.7759e-12*pfpt-1.1351e-10)*(pfps-35.)+((-5.4481e-14*pfpt    &
915        &    +8.733e-12)*pfpt-6.7795e-10)*pfpt+1.8741e-8))*pfphp             &
916        &    +(-4.2393e-8*pfpt+1.8932e-6)*(pfps-35.)                         &
917        &    +((6.6228e-10*pfpt-6.836e-8)*pfpt+8.5258e-6)*pfpt+3.5803e-5
918      !
919    END FUNCTION fsatg
920    !!======================================================================
921END MODULE sbcisf
Note: See TracBrowser for help on using the repository browser.