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

source: branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 5682

Last change on this file since 5682 was 5682, checked in by mattmartin, 9 years ago

OBS simplification changes committed to branch after running SETTE tests to make sure we get the same results as the trunk for ORCA2_LIM_OBS.

  • Property svn:keywords set to Id
File size: 42.1 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   REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ?
58   REAL(wp), PUBLIC, SAVE ::   kappa  = 1.54e-6_wp    ! phycst ?
59   REAL(wp), PUBLIC, SAVE ::   rhoisf = 920.0_wp      ! phycst ?
60   REAL(wp), PUBLIC, SAVE ::   tsurf  = -20.0_wp      ! phycst ?
61   REAL(wp), PUBLIC, SAVE ::   lfusisf= 0.334e6_wp    ! phycst ?
62
63!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
64   CHARACTER(len=100), PUBLIC ::   cn_dirisf  = './'    !: Root directory for location of ssr files
65   TYPE(FLD_N)       , PUBLIC ::   sn_qisf, sn_fwfisf     !: information about the runoff file to be read
66   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qisf, sf_fwfisf
67   TYPE(FLD_N)       , PUBLIC ::   sn_rnfisf              !: information about the runoff file to be read
68   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnfisf           
69   TYPE(FLD_N)       , PUBLIC ::   sn_depmax_isf, sn_depmin_isf, sn_Leff_isf     !: information about the runoff file to be read
70   
71   !! * Substitutions
72#  include "domzgr_substitute.h90"
73   !!----------------------------------------------------------------------
74   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
75   !! $Id$
76   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
77   !!----------------------------------------------------------------------
78
79CONTAINS
80 
81  SUBROUTINE sbc_isf(kt)
82    INTEGER, INTENT(in)          ::   kt         ! ocean time step
83    INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
84    INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
85    REAL(wp)                     ::   rmin
86    REAL(wp)                     ::   zhk
87    REAL(wp)                     ::   zt_frz, zpress
88    CHARACTER(len=256)           ::   cfisf , cvarzisf, cvarhisf   ! name for isf file
89    CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file
90    CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale
91    INTEGER           ::   ios           ! Local integer output status for namelist read
92      !
93      !!---------------------------------------------------------------------
94      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, ln_divisf, ln_conserve, rn_gammat0, rn_gammas0, nn_gammablk, &
95                         & sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
96      !
97      !
98      !                                         ! ====================== !
99      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
100         !                                      ! ====================== !
101         REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs
102         READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
103901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp )
104
105         REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
106         READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
107902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
108         IF(lwm) WRITE ( numond, namsbc_isf )
109
110
111         IF ( lwp ) WRITE(numout,*)
112         IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf'
113         IF ( lwp ) WRITE(numout,*) '~~~~~~~~~'
114         IF ( lwp ) WRITE(numout,*) 'sbcisf :' 
115         IF ( lwp ) WRITE(numout,*) '~~~~~~~~'
116         IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf
117         IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk
118         IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl
119         IF ( lwp ) WRITE(numout,*) '        ln_divisf   = ', ln_divisf 
120         IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk 
121         IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2 
122         IF (ln_divisf) THEN       ! keep it in the namelist ??? used true anyway as for runoff ? (PM)
123            rdivisf = 1._wp
124         ELSE
125            rdivisf = 0._wp
126         END IF
127         !
128         ! Allocate public variable
129         IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
130         !
131         ! initialisation
132         qisf(:,:)        = 0._wp  ; fwfisf(:,:) = 0._wp
133         risf_tsc(:,:,:)  = 0._wp
134         !
135         ! define isf tbl tickness, top and bottom indice
136         IF      (nn_isf == 1) THEN
137            rhisf_tbl(:,:) = rn_hisf_tbl
138            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
139         ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN
140            ALLOCATE( sf_rnfisf(1), STAT=ierror )
141            ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
142            CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
143
144            !: read effective lenght (BG03)
145            IF (nn_isf == 2) THEN
146               ! Read Data and save some integral values
147               CALL iom_open( sn_Leff_isf%clname, inum )
148               cvarLeff  = 'soLeff'               !: variable name for Efficient Length scale
149               CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
150               CALL iom_close(inum)
151               !
152               risfLeff = risfLeff*1000           !: convertion in m
153            END IF
154
155           ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
156            CALL iom_open( sn_depmax_isf%clname, inum )
157            cvarhisf = TRIM(sn_depmax_isf%clvar)
158            CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
159            CALL iom_close(inum)
160            !
161            CALL iom_open( sn_depmin_isf%clname, inum )
162            cvarzisf = TRIM(sn_depmin_isf%clvar)
163            CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
164            CALL iom_close(inum)
165            !
166            rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer
167
168           !! compute first level of the top boundary layer
169           DO ji = 1, jpi
170              DO jj = 1, jpj
171                  jk = 2
172                  DO WHILE ( jk .LE. mbkt(ji,jj) .AND. fsdepw(ji,jj,jk) < rzisf_tbl(ji,jj) ) ;  jk = jk + 1 ;  END DO
173                  misfkt(ji,jj) = jk-1
174               END DO
175            END DO
176
177         ELSE IF ( nn_isf == 4 ) THEN
178            ! as in nn_isf == 1
179            rhisf_tbl(:,:) = rn_hisf_tbl
180            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
181           
182            ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
183            ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror )
184            ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
185            ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) )
186            CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
187            !CALL fld_fill( sf_qisf  , (/ sn_qisf   /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data'       , 'namsbc_isf' )
188         END IF
189         
190         rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
191
192         ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
193         DO jj = 1,jpj
194            DO ji = 1,jpi
195               ikt = misfkt(ji,jj)
196               ikb = misfkt(ji,jj)
197               ! thickness of boundary layer at least the top level thickness
198               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt))
199
200               ! determine the deepest level influenced by the boundary layer
201               ! test on tmask useless ?????
202               DO jk = ikt, mbkt(ji,jj)
203                  IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
204               END DO
205               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
206               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
207               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
208
209               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
210               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
211            END DO
212         END DO
213         
214      END IF
215
216      !                                            ! ---------------------------------------- !
217      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
218         !                                         ! ---------------------------------------- !
219         fwfisf_b  (:,:  ) = fwfisf  (:,:  )               ! Swap the ocean forcing fields except at nit000
220         risf_tsc_b(:,:,:) = risf_tsc(:,:,:)               ! where before fields are set at the end of the routine
221         !
222      ENDIF
223
224      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN
225
226
227         ! compute salf and heat flux
228         IF (nn_isf == 1) THEN
229            ! realistic ice shelf formulation
230            ! compute T/S/U/V for the top boundary layer
231            CALL sbc_isf_tbl(tsn(:,:,:,jp_tem),ttbl(:,:),'T')
232            CALL sbc_isf_tbl(tsn(:,:,:,jp_sal),stbl(:,:),'T')
233            CALL sbc_isf_tbl(un(:,:,:),utbl(:,:),'U')
234            CALL sbc_isf_tbl(vn(:,:,:),vtbl(:,:),'V')
235            ! iom print
236            CALL iom_put('ttbl',ttbl(:,:))
237            CALL iom_put('stbl',stbl(:,:))
238            CALL iom_put('utbl',utbl(:,:))
239            CALL iom_put('vtbl',vtbl(:,:))
240            ! compute fwf and heat flux
241            CALL sbc_isf_cav (kt)
242
243         ELSE IF (nn_isf == 2) THEN
244            ! Beckmann and Goosse parametrisation
245            stbl(:,:)   = soce
246            CALL sbc_isf_bg03(kt)
247
248         ELSE IF (nn_isf == 3) THEN
249            ! specified runoff in depth (Mathiot et al., XXXX in preparation)
250            CALL fld_read ( kt, nn_fsbc, sf_rnfisf   )
251            fwfisf(:,:) = - sf_rnfisf(1)%fnow(:,:,1)         ! fresh water flux from the isf (fwfisf <0 mean melting)
252            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
253            stbl(:,:)   = soce
254
255         ELSE IF (nn_isf == 4) THEN
256            ! specified fwf and heat flux forcing beneath the ice shelf
257            CALL fld_read ( kt, nn_fsbc, sf_fwfisf   )
258            !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
259            fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf
260            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
261            !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux
262            stbl(:,:)   = soce
263
264         END IF
265         ! compute tsc due to isf
266         ! WARNING water add at temp = 0C, correction term is added, maybe better here but need a 3D variable).
267!         zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04
268         zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress )
269         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - rdivisf * fwfisf(:,:) * zt_frz * r1_rau0 !
270         
271         ! salt effect already take into account in vertical advection
272         risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0
273
274         ! output
275         IF( iom_use('qisf'  ) )   CALL iom_put('qisf'  , qisf)
276         IF( iom_use('fwfisf') )   CALL iom_put('fwfisf', fwfisf * stbl(:,:) / soce )
277
278         ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now
279         fwfisf(:,:) = rdivisf * fwfisf(:,:)         
280 
281         ! lbclnk
282         CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.)
283         CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.)
284         CALL lbc_lnk(fwfisf(:,:)   ,'T',1.)
285         CALL lbc_lnk(qisf(:,:)     ,'T',1.)
286
287         IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
288            IF( ln_rstart .AND.    &                     ! Restart: read in restart file
289                 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN
290               IF(lwp) WRITE(numout,*) '          nit000-1 isf tracer content forcing fields read in the restart file'
291               CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:) )   ! before salt content isf_tsc trend
292               CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b', risf_tsc_b(:,:,jp_sal) )   ! before salt content isf_tsc trend
293               CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b', risf_tsc_b(:,:,jp_tem) )   ! before salt content isf_tsc trend
294            ELSE
295               fwfisf_b(:,:)    = fwfisf(:,:)
296               risf_tsc_b(:,:,:)= risf_tsc(:,:,:)
297            END IF
298         ENDIF
299         !
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                CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, 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      CALL eos_fzp( sss_m(:,:), zfrz(:,:), 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, zgammat2d, zgammas2d )
565      !
566      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav')
567
568   END SUBROUTINE sbc_isf_cav
569
570   SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit )
571      !!----------------------------------------------------------------------
572      !! ** Purpose    : compute the coefficient echange for heat flux
573      !!
574      !! ** Method     : gamma assume constant or depends of u* and stability
575      !!
576      !! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
577      !!                Jenkins et al., 2010, JPO, p2298-2312
578      !!---------------------------------------------------------------------
579      REAL(wp), INTENT(inout) :: gt, gs, zqhisf, zqwisf
580      INTEGER , INTENT(in)    :: ji,jj
581      LOGICAL , INTENT(inout) :: lit
582
583      INTEGER  :: ikt                 ! loop index
584      REAL(wp) :: zut, zvt, zustar           ! U, V at T point and friction velocity
585      REAL(wp) :: zdku, zdkv                 ! U, V shear
586      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number
587      REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
588      REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
589      REAL(wp) :: zhmax                      ! limitation of mol
590      REAL(wp) :: zetastar                   ! stability parameter
591      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence
592      REAL(wp) :: zcoef                      ! temporary coef
593      REAL(wp) :: zdep
594      REAL(wp), PARAMETER :: zxsiN = 0.052   ! dimensionless constant
595      REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number
596      REAL(wp), PARAMETER :: znu   = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1)
597      REAL(wp) ::   rcs      = 1.0e-3_wp        ! conversion: mm/s ==> m/s
598      REAL(wp), DIMENSION(2) :: zts, zab
599      !!---------------------------------------------------------------------
600      !
601      IF( nn_gammablk == 0 ) THEN
602      !! gamma is constant (specified in namelist)
603         gt = rn_gammat0
604         gs = rn_gammas0
605         lit = .FALSE.
606      ELSE IF ( nn_gammablk == 1 ) THEN
607      !! gamma is assume to be proportional to u*
608      !! WARNING in case of Losh 2008 tbl parametrization,
609      !! you have to used the mean value of u in the boundary layer)
610      !! not yet coded
611      !! Jenkins et al., 2010, JPO, p2298-2312
612         ikt = mikt(ji,jj)
613      !! Compute U and V at T points
614   !      zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
615   !      zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
616          zut = utbl(ji,jj)
617          zvt = vtbl(ji,jj)
618
619      !! compute ustar
620         zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
621      !! Compute mean value over the TBL
622
623      !! Compute gammats
624         gt = zustar * rn_gammat0
625         gs = zustar * rn_gammas0
626         lit = .FALSE.
627      ELSE IF ( nn_gammablk == 2 ) THEN
628      !! gamma depends of stability of boundary layer
629      !! WARNING in case of Losh 2008 tbl parametrization,
630      !! you have to used the mean value of u in the boundary layer)
631      !! not yet coded
632      !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
633      !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
634               ikt = mikt(ji,jj)
635
636      !! Compute U and V at T points
637               zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
638               zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
639
640      !! compute ustar
641               zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
642               IF (zustar == 0._wp) THEN           ! only for kt = 1 I think
643                 gt = rn_gammat0
644                 gs = rn_gammas0
645               ELSE
646      !! compute Rc number (as done in zdfric.F90)
647               zcoef = 0.5 / fse3w(ji,jj,ikt)
648               !                                            ! shear of horizontal velocity
649               zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )   &
650                  &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
651               zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )   &
652                  &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
653               !                                            ! richardson number (minimum value set to zero)
654               zRc = rn2(ji,jj,ikt+1) / ( zdku*zdku + zdkv*zdkv + 1.e-20 )
655
656      !! compute Pr and Sc number (can be improved)
657               zPr =   13.8
658               zSc = 2432.0
659
660      !! compute gamma mole
661               zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0
662               zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0
663
664      !! compute bouyancy
665               zts(jp_tem) = ttbl(ji,jj)
666               zts(jp_sal) = stbl(ji,jj)
667               zdep        = fsdepw(ji,jj,ikt)
668               !
669               CALL eos_rab( zts, zdep, zab )
670                  !
671      !! compute length scale
672               zbuofdep = grav * ( zab(jp_tem) * zqhisf - zab(jp_sal) * zqwisf )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
673
674      !! compute Monin Obukov Length
675               ! Maximum boundary layer depth
676               zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) -0.001
677               ! Compute Monin obukhov length scale at the surface and Ekman depth:
678               zmob   = zustar ** 3 / (vkarmn * (zbuofdep + epsln))
679               zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
680
681      !! compute eta* (stability parameter)
682               zetastar = 1 / ( SQRT(1 + MAX(zxsiN * zustar / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0)))
683
684      !! compute the sublayer thickness
685               zhnu = 5 * znu / zustar
686      !! compute gamma turb
687               zgturb = 1/vkarmn * LOG(zustar * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &
688               &      + 1 / ( 2 * zxsiN * zetastar ) - 1/vkarmn
689
690      !! compute gammats
691               gt = zustar / (zgturb + zgmolet)
692               gs = zustar / (zgturb + zgmoles)
693               END IF
694      END IF
695
696   END SUBROUTINE
697
698   SUBROUTINE sbc_isf_tbl( varin, varout, cptin )
699      !!----------------------------------------------------------------------
700      !!                  ***  SUBROUTINE sbc_isf_tbl  ***
701      !!
702      !! ** Purpose : compute mean T/S/U/V in the boundary layer
703      !!
704      !!----------------------------------------------------------------------
705      REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin
706      REAL(wp), DIMENSION(:,:)  , INTENT(out):: varout
707     
708      CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out
709
710      REAL(wp) :: ze3, zhk
711      REAL(wp), DIMENSION(:,:), POINTER :: zikt
712
713      INTEGER :: ji,jj,jk
714      INTEGER :: ikt,ikb
715      INTEGER, DIMENSION(:,:), POINTER :: mkt, mkb
716
717      CALL wrk_alloc( jpi,jpj, mkt, mkb  )
718      CALL wrk_alloc( jpi,jpj, zikt )
719
720      ! get first and last level of tbl
721      mkt(:,:) = misfkt(:,:)
722      mkb(:,:) = misfkb(:,:)
723
724      varout(:,:)=0._wp
725      DO jj = 2,jpj
726         DO ji = 2,jpi
727            IF (ssmask(ji,jj) == 1) THEN
728               ikt = mkt(ji,jj)
729               ikb = mkb(ji,jj)
730
731               ! level fully include in the ice shelf boundary layer
732               DO jk = ikt, ikb - 1
733                  ze3 = fse3t_n(ji,jj,jk)
734                  IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
735                  IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) &
736                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
737                  IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) &
738                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
739               END DO
740
741               ! level partially include in ice shelf boundary layer
742               zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
743               IF (cptin == 'T') &
744                   &  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
745               IF (cptin == 'U') &
746                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk)
747               IF (cptin == 'V') &
748                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk)
749            END IF
750         END DO
751      END DO
752
753      CALL wrk_dealloc( jpi,jpj, mkt, mkb )     
754      CALL wrk_dealloc( jpi,jpj, zikt ) 
755
756      IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.)
757      IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.)
758
759   END SUBROUTINE sbc_isf_tbl
760     
761
762   SUBROUTINE sbc_isf_div( phdivn )
763      !!----------------------------------------------------------------------
764      !!                  ***  SUBROUTINE sbc_isf_div  ***
765      !!       
766      !! ** Purpose :   update the horizontal divergence with the runoff inflow
767      !!
768      !! ** Method  :   
769      !!                CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the
770      !!                          divergence and expressed in m/s
771      !!
772      !! ** Action  :   phdivn   decreased by the runoff inflow
773      !!----------------------------------------------------------------------
774      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
775      !!
776      INTEGER  ::   ji, jj, jk   ! dummy loop indices
777      INTEGER  ::   ikt, ikb 
778      INTEGER  ::   nk_isf
779      REAL(wp)     ::   zhk, z1_hisf_tbl, zhisf_tbl
780      REAL(wp)     ::   zfact     ! local scalar
781      !!----------------------------------------------------------------------
782      !
783      zfact   = 0.5_wp
784      !
785      IF (lk_vvl) THEN     ! need to re compute level distribution of isf fresh water
786         DO jj = 1,jpj
787            DO ji = 1,jpi
788               ikt = misfkt(ji,jj)
789               ikb = misfkt(ji,jj)
790               ! thickness of boundary layer at least the top level thickness
791               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t(ji,jj,ikt))
792
793               ! determine the deepest level influenced by the boundary layer
794               ! test on tmask useless ?????
795               DO jk = ikt, mbkt(ji,jj)
796!                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
797               END DO
798               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
799               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
800               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
801
802               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
803               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
804            END DO
805         END DO
806      END IF ! vvl case
807      !
808      DO jj = 1,jpj
809         DO ji = 1,jpi
810               ikt = misfkt(ji,jj)
811               ikb = misfkb(ji,jj)
812               ! level fully include in the ice shelf boundary layer
813               DO jk = ikt, ikb - 1
814                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
815                    &               * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
816               END DO
817               ! level partially include in ice shelf boundary layer
818               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
819                  &             + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
820            !==   ice shelf melting mass distributed over several levels   ==!
821         END DO
822      END DO
823      !
824   END SUBROUTINE sbc_isf_div
825                       
826   FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti )
827      !!----------------------------------------------------------------------
828      !!                 ***  ROUTINE eos_init  ***
829      !!
830      !! ** Purpose :   Compute the in-situ temperature [Celcius]
831      !!
832      !! ** Method  :   
833      !!
834      !! Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
835      !!----------------------------------------------------------------------
836      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptem   ! potential temperature [Celcius]
837      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu]
838      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ppress ! pressure             [dBar]
839      REAL(wp), DIMENSION(:,:), POINTER           ::   pti    ! in-situ temperature [Celcius]
840!      REAL(wp) :: fsatg
841!      REAL(wp) :: pfps, pfpt, pfphp
842      REAL(wp) :: zt, zs, zp, zh, zq, zxk
843      INTEGER  :: ji, jj            ! dummy loop indices
844      !
845      CALL wrk_alloc( jpi,jpj, pti  )
846      !
847      DO jj=1,jpj
848         DO ji=1,jpi
849            zh = ppress(ji,jj)
850! Theta1
851            zt = ptem(ji,jj)
852            zs = psal(ji,jj)
853            zp = 0.0
854            zxk= zh * fsatg( zs, zt, zp )
855            zt = zt + 0.5 * zxk
856            zq = zxk
857! Theta2
858            zp = zp + 0.5 * zh
859            zxk= zh*fsatg( zs, zt, zp )
860            zt = zt + 0.29289322 * ( zxk - zq )
861            zq = 0.58578644 * zxk + 0.121320344 * zq
862! Theta3
863            zxk= zh * fsatg( zs, zt, zp )
864            zt = zt + 1.707106781 * ( zxk - zq )
865            zq = 3.414213562 * zxk - 4.121320344 * zq
866! Theta4
867            zp = zp + 0.5 * zh
868            zxk= zh * fsatg( zs, zt, zp )
869            pti(ji,jj) = zt + ( zxk - 2.0 * zq ) / 6.0
870         END DO
871      END DO
872      !
873      CALL wrk_dealloc( jpi,jpj, pti  )
874      !
875   END FUNCTION tinsitu
876   !
877   FUNCTION fsatg( pfps, pfpt, pfphp )
878      !!----------------------------------------------------------------------
879      !!                 ***  FUNCTION fsatg  ***
880      !!
881      !! ** Purpose    :   Compute the Adiabatic laspse rate [Celcius].[decibar]^-1
882      !!
883      !! ** Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
884      !!
885      !! ** units      :   pressure        pfphp    decibars
886      !!                   temperature     pfpt     deg celsius (ipts-68)
887      !!                   salinity        pfps     (ipss-78)
888      !!                   adiabatic       fsatg    deg. c/decibar
889      !!----------------------------------------------------------------------
890      REAL(wp) :: pfps, pfpt, pfphp 
891      REAL(wp) :: fsatg
892      !
893      fsatg = (((-2.1687e-16*pfpt+1.8676e-14)*pfpt-4.6206e-13)*pfphp         &
894        &    +((2.7759e-12*pfpt-1.1351e-10)*(pfps-35.)+((-5.4481e-14*pfpt    &
895        &    +8.733e-12)*pfpt-6.7795e-10)*pfpt+1.8741e-8))*pfphp             &
896        &    +(-4.2393e-8*pfpt+1.8932e-6)*(pfps-35.)                         &
897        &    +((6.6228e-10*pfpt-6.836e-8)*pfpt+8.5258e-6)*pfpt+3.5803e-5
898      !
899    END FUNCTION fsatg
900    !!======================================================================
901END MODULE sbcisf
Note: See TracBrowser for help on using the repository browser.