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/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/2015/dev_r5204_CNRS_PISCES_dcy/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 5222

Last change on this file since 5222 was 5120, checked in by acc, 9 years ago

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