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

source: branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 5834

Last change on this file since 5834 was 5834, checked in by mathiot, 9 years ago

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