New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcisf.F90 in branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 7982

Last change on this file since 7982 was 7982, checked in by mathiot, 7 years ago

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