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

Last change on this file since 7607 was 7607, checked in by cetlod, 4 years ago

v3.6 stable : add missing features for CMIP6 exercise, see ticket #1834

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