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 @ 7963

Last change on this file since 7963 was 7963, checked in by clem, 7 years ago

solve ticket #1889
  • Property svn:keywords set to Id
File size: 44.9 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(:,:) * stbl(:,:) * r1_rau0
171
172         ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now
173         fwfisf(:,:) = rdivisf * fwfisf(:,:)         
174 
175         ! lbclnk
176         CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.)
177         CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.)
178         CALL lbc_lnk(fwfisf(:,:)   ,'T',1.)
179         CALL lbc_lnk(qisf(:,:)     ,'T',1.)
180
181         ! output
182         IF( iom_use('iceshelf_cea') )   CALL iom_put( 'iceshelf_cea', -fwfisf(:,:)                      )   ! isf mass flux
183         IF( iom_use('hflx_isf_cea') )   CALL iom_put( 'hflx_isf_cea', risf_tsc(:,:,jp_tem) * rau0 * rcp )   ! isf sensible+latent heat (W/m2)
184         IF( iom_use('qlatisf' ) )       CALL iom_put( 'qlatisf'     , qisf(:,:)                         )   ! isf latent heat
185         IF( iom_use('fwfisf'  ) )       CALL iom_put( 'fwfisf'      , fwfisf(:,:)                       )   ! isf mass flux (opposite sign)
186
187         ! Diagnostics
188         IF( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN
189            !
190            CALL wrk_alloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d )
191            CALL wrk_alloc( jpi,jpj,     zqhcisf2d                        )
192            !
193            zfwfisf3d(:,:,:) = 0.0_wp                         ! 3d ice shelf melting (kg/m2/s)
194            zqhcisf3d(:,:,:) = 0.0_wp                         ! 3d heat content flux (W/m2)
195            zqlatisf3d(:,:,:)= 0.0_wp                         ! 3d ice shelf melting latent heat flux (W/m2)
196            zqhcisf2d(:,:)   = fwfisf(:,:) * zt_frz * rcp     ! 2d heat content flux (W/m2)
197            !
198            DO jj = 1,jpj
199               DO ji = 1,jpi
200                  ikt = misfkt(ji,jj)
201                  ikb = misfkb(ji,jj)
202                  DO jk = ikt, ikb - 1
203                     zhisf = r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
204                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf(ji,jj)    * zhisf
205                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * zhisf
206                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf(ji,jj)      * zhisf
207                  END DO
208                  jk = ikb
209                  zhisf = r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
210                  zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * zhisf * ralpha(ji,jj) 
211                  zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * zhisf * ralpha(ji,jj)
212                  zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * zhisf * ralpha(ji,jj)
213               END DO
214            END DO
215            !
216            CALL iom_put( 'fwfisf3d' , zfwfisf3d (:,:,:) )
217            CALL iom_put( 'qlatisf3d', zqlatisf3d(:,:,:) )
218            CALL iom_put( 'qhcisf3d' , zqhcisf3d (:,:,:) )
219            CALL iom_put( 'qhcisf'   , zqhcisf2d (:,:  ) )
220            !
221            CALL wrk_dealloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d )
222            CALL wrk_dealloc( jpi,jpj,     zqhcisf2d                        )
223            !
224         END IF
225         !
226      END IF
227      !
228      !
229      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
230         IF( ln_rstart .AND.    &                     ! Restart: read in restart file
231              & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN
232            IF(lwp) WRITE(numout,*) '          nit000-1 isf tracer content forcing fields read in the restart file'
233            CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:) ) ! before salt content isf_tsc trend
234            CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b', risf_tsc_b(:,:,jp_sal) )   ! before salt content isf_tsc trend
235            CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b', risf_tsc_b(:,:,jp_tem) )   ! before salt content isf_tsc trend
236         ELSE
237            fwfisf_b(:,:)    = fwfisf(:,:)
238            risf_tsc_b(:,:,:)= risf_tsc(:,:,:)
239         END IF
240      ENDIF
241      !
242      IF( lrst_oce ) THEN
243         IF(lwp) WRITE(numout,*)
244         IF(lwp) WRITE(numout,*) 'sbc : isf surface tracer content forcing fields written in ocean restart file ',   &
245            &                    'at it= ', kt,' date= ', ndastp
246         IF(lwp) WRITE(numout,*) '~~~~'
247         CALL iom_rstput( kt, nitrst, numrow, 'fwf_isf_b', fwfisf(:,:) )
248         CALL iom_rstput( kt, nitrst, numrow, 'isf_hc_b' , risf_tsc(:,:,jp_tem) )
249         CALL iom_rstput( kt, nitrst, numrow, 'isf_sc_b' , risf_tsc(:,:,jp_sal) )
250      ENDIF
251       !
252  END SUBROUTINE sbc_isf
253
254  SUBROUTINE sbc_isf_init
255
256    INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
257    INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
258    REAL(wp)                     ::   zhk
259    CHARACTER(len=256)           ::   cfisf , cvarzisf, cvarhisf   ! name for isf file
260    CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file
261    CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale
262    INTEGER           ::   ios           ! Local integer output status for namelist read
263
264      !
265      !!---------------------------------------------------------------------
266      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, ln_divisf, ln_conserve, rn_gammat0, rn_gammas0, nn_gammablk, &
267                         & sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
268      !
269      !
270         REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs
271         READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
272901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp )
273
274         REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
275         READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
276902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
277         IF(lwm) WRITE ( numond, namsbc_isf )
278
279
280         IF ( lwp ) WRITE(numout,*)
281         IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf'
282         IF ( lwp ) WRITE(numout,*) '~~~~~~~~~'
283         IF ( lwp ) WRITE(numout,*) 'sbcisf :' 
284         IF ( lwp ) WRITE(numout,*) '~~~~~~~~'
285         IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf
286         IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk
287         IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl
288         IF ( lwp ) WRITE(numout,*) '        ln_divisf   = ', ln_divisf 
289         IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk 
290         IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2 
291         IF (ln_divisf) THEN       ! keep it in the namelist ??? used true anyway as for runoff ? (PM)
292            rdivisf = 1._wp
293         ELSE
294            rdivisf = 0._wp
295         END IF
296         !
297         ! Allocate public variable
298         IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
299         !
300         ! initialisation
301         qisf(:,:)        = 0._wp  ; fwfisf(:,:) = 0._wp
302         risf_tsc(:,:,:)  = 0._wp
303         !
304         ! define isf tbl tickness, top and bottom indice
305         IF      (nn_isf == 1) THEN
306            rhisf_tbl(:,:) = rn_hisf_tbl
307            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
308         ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN
309            IF( .NOT.l_isfcpl ) THEN
310               ALLOCATE( sf_rnfisf(1), STAT=ierror )
311               ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
312               CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
313             ENDIF
314
315            !: read effective lenght (BG03)
316            IF (nn_isf == 2) THEN
317               ! Read Data and save some integral values
318               CALL iom_open( sn_Leff_isf%clname, inum )
319               cvarLeff  = 'soLeff'               !: variable name for Efficient Length scale
320               CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
321               CALL iom_close(inum)
322               !
323               risfLeff = risfLeff*1000           !: convertion in m
324            END IF
325
326           ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
327            CALL iom_open( sn_depmax_isf%clname, inum )
328            cvarhisf = TRIM(sn_depmax_isf%clvar)
329            CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
330            CALL iom_close(inum)
331            !
332            CALL iom_open( sn_depmin_isf%clname, inum )
333            cvarzisf = TRIM(sn_depmin_isf%clvar)
334            CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
335            CALL iom_close(inum)
336            !
337            rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer
338
339           !! compute first level of the top boundary layer
340           DO ji = 1, jpi
341              DO jj = 1, jpj
342                  jk = 2
343                  DO WHILE ( jk .LE. mbkt(ji,jj) .AND. gdepw_0(ji,jj,jk) < rzisf_tbl(ji,jj) ) ;  jk = jk + 1 ;  END DO
344                  misfkt(ji,jj) = jk-1
345               END DO
346            END DO
347
348         ELSE IF ( nn_isf == 4 ) THEN
349            ! as in nn_isf == 1
350            rhisf_tbl(:,:) = rn_hisf_tbl
351            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
352           
353            ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
354            IF( .NOT.l_isfcpl ) THEN
355               ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror )
356               ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
357               ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) )
358               CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
359               !CALL fld_fill( sf_qisf  , (/ sn_qisf   /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data'       , 'namsbc_isf' )
360            ENDIF
361         END IF
362         ! save initial top boundary layer thickness         
363         rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
364         !
365   END SUBROUTINE sbc_isf_init
366     
367
368
369  INTEGER FUNCTION sbc_isf_alloc()
370      !!----------------------------------------------------------------------
371      !!               ***  FUNCTION sbc_isf_rnf_alloc  ***
372      !!----------------------------------------------------------------------
373      sbc_isf_alloc = 0       ! set to zero if no array to be allocated
374      IF( .NOT. ALLOCATED( qisf ) ) THEN
375         ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj)   , &
376               &    rhisf_tbl(jpi,jpj)    , r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj)  , &
377               &    ttbl(jpi,jpj)         , stbl(jpi,jpj)       , utbl(jpi,jpj)       , &
378               &    vtbl(jpi, jpj)        , risfLeff(jpi,jpj)   , rhisf_tbl_0(jpi,jpj), &
379               &    ralpha(jpi,jpj)       , misfkt(jpi,jpj)     , misfkb(jpi,jpj)     , &
380               &    STAT= sbc_isf_alloc )
381         !
382         IF( lk_mpp                  )   CALL mpp_sum ( sbc_isf_alloc )
383         IF( sbc_isf_alloc /= 0 )   CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.')
384         !
385      ENDIF
386  END FUNCTION
387
388  SUBROUTINE sbc_isf_bg03(kt)
389   !!==========================================================================
390   !!                 *** SUBROUTINE sbcisf_bg03  ***
391   !! add net heat and fresh water flux from ice shelf melting
392   !! into the adjacent ocean using the parameterisation by
393   !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
394   !!     interaction for climate models", Ocean Modelling 5(2003) 157-170.
395   !!  (hereafter BG)
396   !!==========================================================================
397   !!----------------------------------------------------------------------
398   !!   sbc_isf_bg03      : routine called from sbcmod
399   !!----------------------------------------------------------------------
400   !!
401   !! ** Purpose   :   Add heat and fresh water fluxes due to ice shelf melting
402   !! ** Reference :   Beckmann et Goosse, 2003, Ocean Modelling
403   !!
404   !! History :
405   !!      !  06-02  (C. Wang) Original code
406   !!----------------------------------------------------------------------
407
408    INTEGER, INTENT ( in ) :: kt
409
410    INTEGER :: ji, jj, jk, jish  !temporary integer
411    INTEGER :: ijkmin
412    INTEGER :: ii, ij, ik 
413    INTEGER :: inum
414
415    REAL(wp) :: zt_sum      ! sum of the temperature between 200m and 600m
416    REAL(wp) :: zt_ave      ! averaged temperature between 200m and 600m
417    REAL(wp) :: zt_frz      ! freezing point temperature at depth z
418    REAL(wp) :: zpress      ! pressure to compute the freezing point in depth
419   
420    !!----------------------------------------------------------------------
421    IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03')
422     !
423
424    ! This test is false only in the very first time step of a run (JMM ???- Initialy build to skip 1rst year of run )
425    DO ji = 1, jpi
426       DO jj = 1, jpj
427          ik = misfkt(ji,jj)
428          !! Initialize arrays to 0 (each step)
429          zt_sum = 0.e0_wp
430          IF ( ik .GT. 1 ) THEN
431    ! 3. -----------the average temperature between 200m and 600m ---------------------
432             DO jk = misfkt(ji,jj),misfkb(ji,jj)
433             ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!)
434             ! after verif with UNESCO, wrong sign in BG eq. 2
435             ! Calculate freezing temperature
436                zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04 
437                CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, zpress) 
438                zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik)  ! sum temp
439             ENDDO
440             zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
441   
442    ! 4. ------------Net heat flux and fresh water flux due to the ice shelf
443          ! For those corresponding to zonal boundary   
444             qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
445                         & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,ik) 
446             
447             fwfisf(ji,jj) = qisf(ji,jj) / lfusisf          !fresh water flux kg/(m2s)                 
448             fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
449             !add to salinity trend
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  * ( soce / MAX(stbl(ji,jj),zeps))
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! zhtflx is upward heat flux (out of ocean)
587! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
588                     zhtflx = ( zgammat*rau0 - zcfac*zfwflx ) * rcp * (zti(ji,jj) - zfrz(ji,jj) ) 
589! zwflx is upward water flux
590! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
591                     zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj)
592! test convergence and compute gammat
593                     IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE.
594
595                     nit = nit + 1
596                     IF (nit .GE. 51) THEN
597                        WRITE(numout,*) "sbcisf : too many iteration ... ", &
598                            &  zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea
599                        CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
600                     END IF
601! save gammat and compute zhtflx_b
602                     zgammat2d(ji,jj)=zgammat
603                     zgammas2d(ji,jj)=zgammas
604                     zhtflx_b = zhtflx
605
606                  END DO
607! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
608                  qisf(ji,jj) = - zhtflx 
609! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
610                  fwfisf(ji,jj) = zfwflx 
611               ELSE
612                  fwfisf(ji,jj) = 0._wp
613                  qisf(ji,jj)   = 0._wp
614               ENDIF
615               !
616            END DO
617         END DO
618      ENDIF
619      ! lbclnk
620      CALL lbc_lnk(zgammas2d(:,:),'T',1.)
621      CALL lbc_lnk(zgammat2d(:,:),'T',1.)
622      ! output
623      CALL iom_put('isfgammat', zgammat2d)
624      CALL iom_put('isfgammas', zgammas2d)
625      !
626      CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
627      !
628      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav')
629
630   END SUBROUTINE sbc_isf_cav
631
632   SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit )
633      !!----------------------------------------------------------------------
634      !! ** Purpose    : compute the coefficient echange for heat flux
635      !!
636      !! ** Method     : gamma assume constant or depends of u* and stability
637      !!
638      !! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
639      !!                Jenkins et al., 2010, JPO, p2298-2312
640      !!---------------------------------------------------------------------
641      REAL(wp), INTENT(inout) :: gt, gs, zqhisf, zqwisf
642      INTEGER , INTENT(in)    :: ji,jj
643      LOGICAL , INTENT(inout) :: lit
644
645      INTEGER  :: ikt                 ! loop index
646      REAL(wp) :: zut, zvt, zustar           ! U, V at T point and friction velocity
647      REAL(wp) :: zdku, zdkv                 ! U, V shear
648      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number
649      REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
650      REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
651      REAL(wp) :: zhmax                      ! limitation of mol
652      REAL(wp) :: zetastar                   ! stability parameter
653      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence
654      REAL(wp) :: zcoef                      ! temporary coef
655      REAL(wp) :: zdep
656      REAL(wp), PARAMETER :: zxsiN = 0.052   ! dimensionless constant
657      REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number
658      REAL(wp), PARAMETER :: znu   = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1)
659      REAL(wp) ::   rcs      = 1.0e-3_wp        ! conversion: mm/s ==> m/s
660      REAL(wp), DIMENSION(2) :: zts, zab
661      !!---------------------------------------------------------------------
662      !
663      IF( nn_gammablk == 0 ) THEN
664      !! gamma is constant (specified in namelist)
665         gt = rn_gammat0
666         gs = rn_gammas0
667         lit = .FALSE.
668      ELSE IF ( nn_gammablk == 1 ) THEN
669      !! gamma is assume to be proportional to u*
670      !! WARNING in case of Losh 2008 tbl parametrization,
671      !! you have to used the mean value of u in the boundary layer)
672      !! not yet coded
673      !! Jenkins et al., 2010, JPO, p2298-2312
674         ikt = mikt(ji,jj)
675      !! Compute U and V at T points
676   !      zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
677   !      zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
678          zut = utbl(ji,jj)
679          zvt = vtbl(ji,jj)
680
681      !! compute ustar
682         zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
683      !! Compute mean value over the TBL
684
685      !! Compute gammats
686         gt = zustar * rn_gammat0
687         gs = zustar * rn_gammas0
688         lit = .FALSE.
689      ELSE IF ( nn_gammablk == 2 ) THEN
690      !! gamma depends of stability of boundary layer
691      !! WARNING in case of Losh 2008 tbl parametrization,
692      !! you have to used the mean value of u in the boundary layer)
693      !! not yet coded
694      !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
695      !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
696               ikt = mikt(ji,jj)
697
698      !! Compute U and V at T points
699               zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
700               zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
701
702      !! compute ustar
703               zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
704               IF (zustar == 0._wp) THEN           ! only for kt = 1 I think
705                 gt = rn_gammat0
706                 gs = rn_gammas0
707               ELSE
708      !! compute Rc number (as done in zdfric.F90)
709               zcoef = 0.5 / fse3w(ji,jj,ikt)
710               !                                            ! shear of horizontal velocity
711               zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )   &
712                  &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
713               zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )   &
714                  &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
715               !                                            ! richardson number (minimum value set to zero)
716               zRc = rn2(ji,jj,ikt+1) / ( zdku*zdku + zdkv*zdkv + 1.e-20 )
717
718      !! compute Pr and Sc number (can be improved)
719               zPr =   13.8
720               zSc = 2432.0
721
722      !! compute gamma mole
723               zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0
724               zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0
725
726      !! compute bouyancy
727               zts(jp_tem) = ttbl(ji,jj)
728               zts(jp_sal) = stbl(ji,jj)
729               zdep        = fsdepw(ji,jj,ikt)
730               !
731               CALL eos_rab( zts, zdep, zab )
732                  !
733      !! compute length scale
734               zbuofdep = grav * ( zab(jp_tem) * zqhisf - zab(jp_sal) * zqwisf )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
735
736      !! compute Monin Obukov Length
737               ! Maximum boundary layer depth
738               zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) -0.001
739               ! Compute Monin obukhov length scale at the surface and Ekman depth:
740               zmob   = zustar ** 3 / (vkarmn * (zbuofdep + epsln))
741               zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
742
743      !! compute eta* (stability parameter)
744               zetastar = 1 / ( SQRT(1 + MAX(zxsiN * zustar / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0)))
745
746      !! compute the sublayer thickness
747               zhnu = 5 * znu / zustar
748      !! compute gamma turb
749               zgturb = 1/vkarmn * LOG(zustar * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &
750               &      + 1 / ( 2 * zxsiN * zetastar ) - 1/vkarmn
751
752      !! compute gammats
753               gt = zustar / (zgturb + zgmolet)
754               gs = zustar / (zgturb + zgmoles)
755               END IF
756      END IF
757
758   END SUBROUTINE
759
760   SUBROUTINE sbc_isf_tbl( varin, varout, cptin )
761      !!----------------------------------------------------------------------
762      !!                  ***  SUBROUTINE sbc_isf_tbl  ***
763      !!
764      !! ** Purpose : compute mean T/S/U/V in the boundary layer
765      !!
766      !!----------------------------------------------------------------------
767      REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin
768      REAL(wp), DIMENSION(:,:)  , INTENT(out):: varout
769     
770      CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out
771
772      REAL(wp) :: ze3, zhk
773      REAL(wp), DIMENSION(:,:), POINTER :: zikt
774
775      INTEGER :: ji,jj,jk
776      INTEGER :: ikt,ikb
777      INTEGER, DIMENSION(:,:), POINTER :: mkt, mkb
778
779      CALL wrk_alloc( jpi,jpj, mkt, mkb  )
780      CALL wrk_alloc( jpi,jpj, zikt )
781
782      ! get first and last level of tbl
783      mkt(:,:) = misfkt(:,:)
784      mkb(:,:) = misfkb(:,:)
785
786      varout(:,:)=0._wp
787      DO jj = 2,jpj
788         DO ji = 2,jpi
789            IF (ssmask(ji,jj) == 1) THEN
790               ikt = mkt(ji,jj)
791               ikb = mkb(ji,jj)
792
793               ! level fully include in the ice shelf boundary layer
794               DO jk = ikt, ikb - 1
795                  ze3 = fse3t_n(ji,jj,jk)
796                  IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
797                  IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) &
798                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
799                  IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) &
800                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
801               END DO
802
803               ! level partially include in ice shelf boundary layer
804               zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
805               IF (cptin == 'T') &
806                   &  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
807               IF (cptin == 'U') &
808                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk)
809               IF (cptin == 'V') &
810                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk)
811            END IF
812         END DO
813      END DO
814
815      CALL wrk_dealloc( jpi,jpj, mkt, mkb )     
816      CALL wrk_dealloc( jpi,jpj, zikt ) 
817
818      IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.)
819      IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.)
820
821   END SUBROUTINE sbc_isf_tbl
822     
823
824   SUBROUTINE sbc_isf_div( phdivn )
825      !!----------------------------------------------------------------------
826      !!                  ***  SUBROUTINE sbc_isf_div  ***
827      !!       
828      !! ** Purpose :   update the horizontal divergence with the runoff inflow
829      !!
830      !! ** Method  :   
831      !!                CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the
832      !!                          divergence and expressed in m/s
833      !!
834      !! ** Action  :   phdivn   decreased by the runoff inflow
835      !!----------------------------------------------------------------------
836      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
837      !!
838      INTEGER  ::   ji, jj, jk   ! dummy loop indices
839      INTEGER  ::   ikt, ikb 
840      INTEGER  ::   nk_isf
841      REAL(wp)     ::   zhk, z1_hisf_tbl, zhisf_tbl
842      REAL(wp)     ::   zfact     ! local scalar
843      !!----------------------------------------------------------------------
844      !
845      zfact   = 0.5_wp
846      !
847      IF (lk_vvl) THEN     ! need to re compute level distribution of isf fresh water
848         DO jj = 1,jpj
849            DO ji = 1,jpi
850               ikt = misfkt(ji,jj)
851               ikb = misfkt(ji,jj)
852               ! thickness of boundary layer at least the top level thickness
853               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t(ji,jj,ikt))
854
855               ! determine the deepest level influenced by the boundary layer
856               ! test on tmask useless ?????
857               DO jk = ikt, mbkt(ji,jj)
858                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
859               END DO
860               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
861               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
862               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
863
864               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
865               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
866            END DO
867         END DO
868      END IF ! vvl case
869      !
870      DO jj = 1,jpj
871         DO ji = 1,jpi
872               ikt = misfkt(ji,jj)
873               ikb = misfkb(ji,jj)
874               ! level fully include in the ice shelf boundary layer
875               DO jk = ikt, ikb - 1
876                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
877                    &               * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
878               END DO
879               ! level partially include in ice shelf boundary layer
880               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
881                  &             + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
882            !==   ice shelf melting mass distributed over several levels   ==!
883         END DO
884      END DO
885      !
886   END SUBROUTINE sbc_isf_div
887                       
888   FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti )
889      !!----------------------------------------------------------------------
890      !!                 ***  ROUTINE eos_init  ***
891      !!
892      !! ** Purpose :   Compute the in-situ temperature [Celcius]
893      !!
894      !! ** Method  :   
895      !!
896      !! Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
897      !!----------------------------------------------------------------------
898      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptem   ! potential temperature [Celcius]
899      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu]
900      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ppress ! pressure             [dBar]
901      REAL(wp), DIMENSION(:,:), POINTER           ::   pti    ! in-situ temperature [Celcius]
902!      REAL(wp) :: fsatg
903!      REAL(wp) :: pfps, pfpt, pfphp
904      REAL(wp) :: zt, zs, zp, zh, zq, zxk
905      INTEGER  :: ji, jj            ! dummy loop indices
906      !
907      CALL wrk_alloc( jpi,jpj, pti  )
908      !
909      DO jj=1,jpj
910         DO ji=1,jpi
911            zh = ppress(ji,jj)
912! Theta1
913            zt = ptem(ji,jj)
914            zs = psal(ji,jj)
915            zp = 0.0
916            zxk= zh * fsatg( zs, zt, zp )
917            zt = zt + 0.5 * zxk
918            zq = zxk
919! Theta2
920            zp = zp + 0.5 * zh
921            zxk= zh*fsatg( zs, zt, zp )
922            zt = zt + 0.29289322 * ( zxk - zq )
923            zq = 0.58578644 * zxk + 0.121320344 * zq
924! Theta3
925            zxk= zh * fsatg( zs, zt, zp )
926            zt = zt + 1.707106781 * ( zxk - zq )
927            zq = 3.414213562 * zxk - 4.121320344 * zq
928! Theta4
929            zp = zp + 0.5 * zh
930            zxk= zh * fsatg( zs, zt, zp )
931            pti(ji,jj) = zt + ( zxk - 2.0 * zq ) / 6.0
932         END DO
933      END DO
934      !
935      CALL wrk_dealloc( jpi,jpj, pti  )
936      !
937   END FUNCTION tinsitu
938   !
939   FUNCTION fsatg( pfps, pfpt, pfphp )
940      !!----------------------------------------------------------------------
941      !!                 ***  FUNCTION fsatg  ***
942      !!
943      !! ** Purpose    :   Compute the Adiabatic laspse rate [Celcius].[decibar]^-1
944      !!
945      !! ** Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
946      !!
947      !! ** units      :   pressure        pfphp    decibars
948      !!                   temperature     pfpt     deg celsius (ipts-68)
949      !!                   salinity        pfps     (ipss-78)
950      !!                   adiabatic       fsatg    deg. c/decibar
951      !!----------------------------------------------------------------------
952      REAL(wp) :: pfps, pfpt, pfphp 
953      REAL(wp) :: fsatg
954      !
955      fsatg = (((-2.1687e-16*pfpt+1.8676e-14)*pfpt-4.6206e-13)*pfphp         &
956        &    +((2.7759e-12*pfpt-1.1351e-10)*(pfps-35.)+((-5.4481e-14*pfpt    &
957        &    +8.733e-12)*pfpt-6.7795e-10)*pfpt+1.8741e-8))*pfphp             &
958        &    +(-4.2393e-8*pfpt+1.8932e-6)*(pfps-35.)                         &
959        &    +((6.6228e-10*pfpt-6.836e-8)*pfpt+8.5258e-6)*pfpt+3.5803e-5
960      !
961    END FUNCTION fsatg
962    !!======================================================================
963END MODULE sbcisf
Note: See TracBrowser for help on using the repository browser.