source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

File size: 41.8 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
58   REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ?
59   REAL(wp), PUBLIC, SAVE ::   kappa  = 1.54e-6_wp    ! phycst ?
60   REAL(wp), PUBLIC, SAVE ::   rhoisf = 920.0_wp      ! phycst ?
61   REAL(wp), PUBLIC, SAVE ::   tsurf  = -20.0_wp      ! phycst ?
62   REAL(wp), PUBLIC, SAVE ::   lfusisf= 0.334e6_wp    ! phycst ?
63
64!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
65   CHARACTER(len=100), PUBLIC ::   cn_dirisf  = './'    !: Root directory for location of ssr files
66   TYPE(FLD_N)       , PUBLIC ::   sn_qisf, sn_fwfisf     !: information about the runoff file to be read
67   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qisf, sf_fwfisf
68   TYPE(FLD_N)       , PUBLIC ::   sn_rnfisf              !: information about the runoff file to be read
69   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnfisf           
70   TYPE(FLD_N)       , PUBLIC ::   sn_depmax_isf, sn_depmin_isf, sn_Leff_isf     !: information about the runoff file to be read
71   
72   !! * Substitutions
73#  include "domzgr_substitute.h90"
74   !!----------------------------------------------------------------------
75   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
76   !! $Id$
77   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
78   !!----------------------------------------------------------------------
79
80CONTAINS
81 
82  SUBROUTINE sbc_isf(kt)
83    INTEGER, INTENT(in)          ::   kt         ! ocean time step
84    INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
85    INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
86    REAL(wp)                     ::   rmin
87    REAL(wp)                     ::   zhk
88    REAL(wp)                     ::   zt_frz, zpress
89    CHARACTER(len=256)           ::   cfisf , cvarzisf, cvarhisf   ! name for isf file
90    CHARACTER(LEN=256)           :: cnameis                     ! name of iceshelf file
91    CHARACTER (LEN=32)           :: cvarLeff                    ! variable name for efficient Length scale
92    INTEGER           ::   ios           ! Local integer output status for namelist read
93      !
94      !!---------------------------------------------------------------------
95      NAMELIST/namsbc_isf/ nn_isfblk, rn_hisf_tbl, ln_divisf, ln_conserve, rn_gammat0, rn_gammas0, nn_gammablk, &
96                         & sn_fwfisf, sn_qisf, sn_rnfisf, sn_depmax_isf, sn_depmin_isf, sn_Leff_isf
97      !
98      !
99      !                                         ! ====================== !
100      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
101         !                                      ! ====================== !
102         REWIND( numnam_ref )              ! Namelist namsbc_rnf in reference namelist : Runoffs
103         READ  ( numnam_ref, namsbc_isf, IOSTAT = ios, ERR = 901)
104901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in reference namelist', lwp )
105
106         REWIND( numnam_cfg )              ! Namelist namsbc_rnf in configuration namelist : Runoffs
107         READ  ( numnam_cfg, namsbc_isf, IOSTAT = ios, ERR = 902 )
108902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_isf in configuration namelist', lwp )
109         IF(lwm) WRITE ( numond, namsbc_isf )
110
111
112         IF ( lwp ) WRITE(numout,*)
113         IF ( lwp ) WRITE(numout,*) 'sbc_isf: heat flux of the ice shelf'
114         IF ( lwp ) WRITE(numout,*) '~~~~~~~~~'
115         IF ( lwp ) WRITE(numout,*) 'sbcisf :' 
116         IF ( lwp ) WRITE(numout,*) '~~~~~~~~'
117         IF ( lwp ) WRITE(numout,*) '        nn_isf      = ', nn_isf
118         IF ( lwp ) WRITE(numout,*) '        nn_isfblk   = ', nn_isfblk
119         IF ( lwp ) WRITE(numout,*) '        rn_hisf_tbl = ', rn_hisf_tbl
120         IF ( lwp ) WRITE(numout,*) '        ln_divisf   = ', ln_divisf 
121         IF ( lwp ) WRITE(numout,*) '        nn_gammablk = ', nn_gammablk 
122         IF ( lwp ) WRITE(numout,*) '        rn_tfri2    = ', rn_tfri2 
123         IF (ln_divisf) THEN       ! keep it in the namelist ??? used true anyway as for runoff ? (PM)
124            rdivisf = 1._wp
125         ELSE
126            rdivisf = 0._wp
127         END IF
128         !
129         ! Allocate public variable
130         IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' )
131         !
132         ! initialisation
133         qisf(:,:)        = 0._wp  ; fwfisf(:,:) = 0._wp
134         risf_tsc(:,:,:)  = 0._wp
135         !
136         ! define isf tbl tickness, top and bottom indice
137         IF      (nn_isf == 1) THEN
138            rhisf_tbl(:,:) = rn_hisf_tbl
139            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
140         ELSE IF ((nn_isf == 3) .OR. (nn_isf == 2)) THEN
141            ALLOCATE( sf_rnfisf(1), STAT=ierror )
142            ALLOCATE( sf_rnfisf(1)%fnow(jpi,jpj,1), sf_rnfisf(1)%fdta(jpi,jpj,1,2) )
143            CALL fld_fill( sf_rnfisf, (/ sn_rnfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
144
145            !: read effective lenght (BG03)
146            IF (nn_isf == 2) THEN
147               ! Read Data and save some integral values
148               CALL iom_open( sn_Leff_isf%clname, inum )
149               cvarLeff  = 'soLeff'               !: variable name for Efficient Length scale
150               CALL iom_get( inum, jpdom_data, cvarLeff, risfLeff , 1)
151               CALL iom_close(inum)
152               !
153               risfLeff = risfLeff*1000           !: convertion in m
154            END IF
155
156           ! read depth of the top and bottom of the isf top boundary layer (in this case, isf front depth and grounding line depth)
157            CALL iom_open( sn_depmax_isf%clname, inum )
158            cvarhisf = TRIM(sn_depmax_isf%clvar)
159            CALL iom_get( inum, jpdom_data, cvarhisf, rhisf_tbl, 1) !: depth of deepest point of the ice shelf base
160            CALL iom_close(inum)
161            !
162            CALL iom_open( sn_depmin_isf%clname, inum )
163            cvarzisf = TRIM(sn_depmin_isf%clvar)
164            CALL iom_get( inum, jpdom_data, cvarzisf, rzisf_tbl, 1) !: depth of shallowest point of the ice shelves base
165            CALL iom_close(inum)
166            !
167            rhisf_tbl(:,:) = rhisf_tbl(:,:) - rzisf_tbl(:,:)        !: tickness isf boundary layer
168
169           !! compute first level of the top boundary layer
170           DO ji = 1, jpi
171              DO jj = 1, jpj
172                  jk = 2
173                  DO WHILE ( jk .LE. mbkt(ji,jj) .AND. gdepw_0(ji,jj,jk) < rzisf_tbl(ji,jj) ) ;  jk = jk + 1 ;  END DO
174                  misfkt(ji,jj) = jk-1
175               END DO
176            END DO
177
178         ELSE IF ( nn_isf == 4 ) THEN
179            ! as in nn_isf == 1
180            rhisf_tbl(:,:) = rn_hisf_tbl
181            misfkt(:,:)    = mikt(:,:)         ! same indice for bg03 et cav => used in isfdiv
182           
183            ! load variable used in fldread (use for temporal interpolation of isf fwf forcing)
184            ALLOCATE( sf_fwfisf(1), sf_qisf(1), STAT=ierror )
185            ALLOCATE( sf_fwfisf(1)%fnow(jpi,jpj,1), sf_fwfisf(1)%fdta(jpi,jpj,1,2) )
186            ALLOCATE( sf_qisf(1)%fnow(jpi,jpj,1), sf_qisf(1)%fdta(jpi,jpj,1,2) )
187            CALL fld_fill( sf_fwfisf, (/ sn_fwfisf /), cn_dirisf, 'sbc_isf_init', 'read fresh water flux isf data', 'namsbc_isf' )
188            !CALL fld_fill( sf_qisf  , (/ sn_qisf   /), cn_dirisf, 'sbc_isf_init', 'read heat flux isf data'       , 'namsbc_isf' )
189         END IF
190         
191         ! save initial top boundary layer thickness         
192         rhisf_tbl_0(:,:) = rhisf_tbl(:,:)
193
194      END IF
195
196      !                                            ! ---------------------------------------- !
197      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
198         !                                         ! ---------------------------------------- !
199         fwfisf_b  (:,:  ) = fwfisf  (:,:  )               ! Swap the ocean forcing fields except at nit000
200         risf_tsc_b(:,:,:) = risf_tsc(:,:,:)               ! where before fields are set at the end of the routine
201         !
202      ENDIF
203
204      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN
205
206         ! compute bottom level of isf tbl and thickness of tbl below the ice shelf
207         DO jj = 1,jpj
208            DO ji = 1,jpi
209               ikt = misfkt(ji,jj)
210               ikb = misfkt(ji,jj)
211               ! thickness of boundary layer at least the top level thickness
212               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t_n(ji,jj,ikt))
213
214               ! determine the deepest level influenced by the boundary layer
215               DO jk = ikt, mbkt(ji,jj)
216                  IF ( (SUM(fse3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
217               END DO
218               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
219               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
220               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
221
222               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
223               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
224            END DO
225         END DO
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) CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
475
476! save gammat and compute zhtflx_b
477                     zgammat2d(ji,jj)=zgammat
478                     zhtflx_b = zhtflx
479                  END DO
480
481                  qisf(ji,jj) = - zhtflx
482! For genuine ISOMIP protocol this should probably be something like
483                  fwfisf(ji,jj) = zfwflx  * ( soce / MAX(stbl(ji,jj),zeps))
484               ELSE
485                  fwfisf(ji,jj) = 0._wp
486                  qisf(ji,jj)   = 0._wp
487               END IF
488            !
489            END DO
490         END DO
491
492      ELSE IF (nn_isfblk == 2 ) THEN
493
494! More complicated 3 equation thermodynamics as in MITgcm
495!CDIR COLLAPSE
496         DO jj = 2, jpj
497            DO ji = 2, jpi
498               IF (mikt(ji,jj) > 1 ) THEN
499                  nit=1; lit=.TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp; zhtflx=0._wp
500                  DO WHILE ( lit )
501                     CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit)
502
503                     zeps1=rcp*rau0*zgammat
504                     zeps2=lfusisf*rau0*zgammas
505                     zeps3=rhoisf*rcpi*kappa/risfdep(ji,jj)
506                     zeps4=zlamb2+zlamb3*risfdep(ji,jj)
507                     zeps6=zeps4-zti(ji,jj)
508                     zeps7=zeps4-tsurf
509                     zaqe=zlamb1 * (zeps1 + zeps3)
510                     zaqer=0.5/zaqe
511                     zbqe=zeps1*zeps6+zeps3*zeps7-zeps2
512                     zcqe=zeps2*stbl(ji,jj)
513                     zdis=zbqe*zbqe-4.0*zaqe*zcqe               
514! Presumably zdis can never be negative because gammas is very small compared to gammat
515                     zsfrz=(-zbqe-SQRT(zdis))*zaqer
516                     IF (zsfrz .lt. 0.0) zsfrz=(-zbqe+SQRT(zdis))*zaqer
517                     zfrz(ji,jj)=zeps4+zlamb1*zsfrz
518 
519! zfwflx is upward water flux
520                     zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / zsfrz )
521! zhtflx is upward heat flux (out of ocean)
522! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
523                     zhtflx = ( zgammat*rau0 - zcfac*zfwflx ) * rcp * (zti(ji,jj) - zfrz(ji,jj) ) 
524! zwflx is upward water flux
525! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
526                     zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj)
527! test convergence and compute gammat
528                     IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE.
529
530                     nit = nit + 1
531                     IF (nit .GE. 51) THEN
532                        WRITE(numout,*) "sbcisf : too many iteration ... ", &
533                            &  zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea
534                        CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
535                     END IF
536! save gammat and compute zhtflx_b
537                     zgammat2d(ji,jj)=zgammat
538                     zgammas2d(ji,jj)=zgammas
539                     zhtflx_b = zhtflx
540
541                  END DO
542! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
543                  qisf(ji,jj) = - zhtflx 
544! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
545                  fwfisf(ji,jj) = zfwflx 
546               ELSE
547                  fwfisf(ji,jj) = 0._wp
548                  qisf(ji,jj)   = 0._wp
549               ENDIF
550               !
551            END DO
552         END DO
553      ENDIF
554      ! lbclnk
555      CALL lbc_lnk(zgammas2d(:,:),'T',1.)
556      CALL lbc_lnk(zgammat2d(:,:),'T',1.)
557      ! output
558      CALL iom_put('isfgammat', zgammat2d)
559      CALL iom_put('isfgammas', zgammas2d)
560      !
561      CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
562      !
563      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav')
564
565   END SUBROUTINE sbc_isf_cav
566
567   SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit )
568      !!----------------------------------------------------------------------
569      !! ** Purpose    : compute the coefficient echange for heat flux
570      !!
571      !! ** Method     : gamma assume constant or depends of u* and stability
572      !!
573      !! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
574      !!                Jenkins et al., 2010, JPO, p2298-2312
575      !!---------------------------------------------------------------------
576      REAL(wp), INTENT(inout) :: gt, gs, zqhisf, zqwisf
577      INTEGER , INTENT(in)    :: ji,jj
578      LOGICAL , INTENT(inout) :: lit
579
580      INTEGER  :: ikt                 ! loop index
581      REAL(wp) :: zut, zvt, zustar           ! U, V at T point and friction velocity
582      REAL(wp) :: zdku, zdkv                 ! U, V shear
583      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number
584      REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
585      REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
586      REAL(wp) :: zhmax                      ! limitation of mol
587      REAL(wp) :: zetastar                   ! stability parameter
588      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence
589      REAL(wp) :: zcoef                      ! temporary coef
590      REAL(wp) :: zdep
591      REAL(wp), PARAMETER :: zxsiN = 0.052   ! dimensionless constant
592      REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number
593      REAL(wp), PARAMETER :: znu   = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1)
594      REAL(wp) ::   rcs      = 1.0e-3_wp        ! conversion: mm/s ==> m/s
595      REAL(wp), DIMENSION(2) :: zts, zab
596      !!---------------------------------------------------------------------
597      !
598      IF( nn_gammablk == 0 ) THEN
599      !! gamma is constant (specified in namelist)
600         gt = rn_gammat0
601         gs = rn_gammas0
602         lit = .FALSE.
603      ELSE IF ( nn_gammablk == 1 ) THEN
604      !! gamma is assume to be proportional to u*
605      !! WARNING in case of Losh 2008 tbl parametrization,
606      !! you have to used the mean value of u in the boundary layer)
607      !! not yet coded
608      !! Jenkins et al., 2010, JPO, p2298-2312
609         ikt = mikt(ji,jj)
610      !! Compute U and V at T points
611   !      zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
612   !      zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
613          zut = utbl(ji,jj)
614          zvt = vtbl(ji,jj)
615
616      !! compute ustar
617         zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
618      !! Compute mean value over the TBL
619
620      !! Compute gammats
621         gt = zustar * rn_gammat0
622         gs = zustar * rn_gammas0
623         lit = .FALSE.
624      ELSE IF ( nn_gammablk == 2 ) THEN
625      !! gamma depends of stability of boundary layer
626      !! WARNING in case of Losh 2008 tbl parametrization,
627      !! you have to used the mean value of u in the boundary layer)
628      !! not yet coded
629      !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
630      !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
631               ikt = mikt(ji,jj)
632
633      !! Compute U and V at T points
634               zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
635               zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
636
637      !! compute ustar
638               zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
639               IF (zustar == 0._wp) THEN           ! only for kt = 1 I think
640                 gt = rn_gammat0
641                 gs = rn_gammas0
642               ELSE
643      !! compute Rc number (as done in zdfric.F90)
644               zcoef = 0.5 / fse3w(ji,jj,ikt)
645               !                                            ! shear of horizontal velocity
646               zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )   &
647                  &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
648               zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )   &
649                  &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
650               !                                            ! richardson number (minimum value set to zero)
651               zRc = rn2(ji,jj,ikt+1) / ( zdku*zdku + zdkv*zdkv + 1.e-20 )
652
653      !! compute Pr and Sc number (can be improved)
654               zPr =   13.8
655               zSc = 2432.0
656
657      !! compute gamma mole
658               zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0
659               zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0
660
661      !! compute bouyancy
662               zts(jp_tem) = ttbl(ji,jj)
663               zts(jp_sal) = stbl(ji,jj)
664               zdep        = fsdepw(ji,jj,ikt)
665               !
666               CALL eos_rab( zts, zdep, zab )
667                  !
668      !! compute length scale
669               zbuofdep = grav * ( zab(jp_tem) * zqhisf - zab(jp_sal) * zqwisf )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
670
671      !! compute Monin Obukov Length
672               ! Maximum boundary layer depth
673               zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) -0.001
674               ! Compute Monin obukhov length scale at the surface and Ekman depth:
675               zmob   = zustar ** 3 / (vkarmn * (zbuofdep + epsln))
676               zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
677
678      !! compute eta* (stability parameter)
679               zetastar = 1 / ( SQRT(1 + MAX(zxsiN * zustar / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0)))
680
681      !! compute the sublayer thickness
682               zhnu = 5 * znu / zustar
683      !! compute gamma turb
684               zgturb = 1/vkarmn * LOG(zustar * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &
685               &      + 1 / ( 2 * zxsiN * zetastar ) - 1/vkarmn
686
687      !! compute gammats
688               gt = zustar / (zgturb + zgmolet)
689               gs = zustar / (zgturb + zgmoles)
690               END IF
691      END IF
692
693   END SUBROUTINE
694
695   SUBROUTINE sbc_isf_tbl( varin, varout, cptin )
696      !!----------------------------------------------------------------------
697      !!                  ***  SUBROUTINE sbc_isf_tbl  ***
698      !!
699      !! ** Purpose : compute mean T/S/U/V in the boundary layer
700      !!
701      !!----------------------------------------------------------------------
702      REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin
703      REAL(wp), DIMENSION(:,:)  , INTENT(out):: varout
704     
705      CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out
706
707      REAL(wp) :: ze3, zhk
708      REAL(wp), DIMENSION(:,:), POINTER :: zikt
709
710      INTEGER :: ji,jj,jk
711      INTEGER :: ikt,ikb
712      INTEGER, DIMENSION(:,:), POINTER :: mkt, mkb
713
714      CALL wrk_alloc( jpi,jpj, mkt, mkb  )
715      CALL wrk_alloc( jpi,jpj, zikt )
716
717      ! get first and last level of tbl
718      mkt(:,:) = misfkt(:,:)
719      mkb(:,:) = misfkb(:,:)
720
721      varout(:,:)=0._wp
722      DO jj = 2,jpj
723         DO ji = 2,jpi
724            IF (ssmask(ji,jj) == 1) THEN
725               ikt = mkt(ji,jj)
726               ikb = mkb(ji,jj)
727
728               ! level fully include in the ice shelf boundary layer
729               DO jk = ikt, ikb - 1
730                  ze3 = fse3t_n(ji,jj,jk)
731                  IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
732                  IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) &
733                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
734                  IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) &
735                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
736               END DO
737
738               ! level partially include in ice shelf boundary layer
739               zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
740               IF (cptin == 'T') &
741                   &  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
742               IF (cptin == 'U') &
743                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk)
744               IF (cptin == 'V') &
745                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk)
746            END IF
747         END DO
748      END DO
749
750      CALL wrk_dealloc( jpi,jpj, mkt, mkb )     
751      CALL wrk_dealloc( jpi,jpj, zikt ) 
752
753      IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.)
754      IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.)
755
756   END SUBROUTINE sbc_isf_tbl
757     
758
759   SUBROUTINE sbc_isf_div( phdivn )
760      !!----------------------------------------------------------------------
761      !!                  ***  SUBROUTINE sbc_isf_div  ***
762      !!       
763      !! ** Purpose :   update the horizontal divergence with the runoff inflow
764      !!
765      !! ** Method  :   
766      !!                CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the
767      !!                          divergence and expressed in m/s
768      !!
769      !! ** Action  :   phdivn   decreased by the runoff inflow
770      !!----------------------------------------------------------------------
771      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
772      !!
773      INTEGER  ::   ji, jj, jk   ! dummy loop indices
774      INTEGER  ::   ikt, ikb 
775      INTEGER  ::   nk_isf
776      REAL(wp)     ::   zhk, z1_hisf_tbl, zhisf_tbl
777      REAL(wp)     ::   zfact     ! local scalar
778      !!----------------------------------------------------------------------
779      !
780      zfact   = 0.5_wp
781      !
782      IF (lk_vvl) THEN     ! need to re compute level distribution of isf fresh water
783         DO jj = 1,jpj
784            DO ji = 1,jpi
785               ikt = misfkt(ji,jj)
786               ikb = misfkt(ji,jj)
787               ! thickness of boundary layer at least the top level thickness
788               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t(ji,jj,ikt))
789
790               ! determine the deepest level influenced by the boundary layer
791               ! test on tmask useless ?????
792               DO jk = ikt, mbkt(ji,jj)
793                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
794               END DO
795               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
796               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
797               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
798
799               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
800               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
801            END DO
802         END DO
803      END IF ! vvl case
804      !
805      DO jj = 1,jpj
806         DO ji = 1,jpi
807               ikt = misfkt(ji,jj)
808               ikb = misfkb(ji,jj)
809               ! level fully include in the ice shelf boundary layer
810               DO jk = ikt, ikb - 1
811                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
812                    &               * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
813               END DO
814               ! level partially include in ice shelf boundary layer
815               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
816                  &             + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
817            !==   ice shelf melting mass distributed over several levels   ==!
818         END DO
819      END DO
820      !
821   END SUBROUTINE sbc_isf_div
822                       
823   FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti )
824      !!----------------------------------------------------------------------
825      !!                 ***  ROUTINE eos_init  ***
826      !!
827      !! ** Purpose :   Compute the in-situ temperature [Celcius]
828      !!
829      !! ** Method  :   
830      !!
831      !! Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
832      !!----------------------------------------------------------------------
833      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptem   ! potential temperature [Celcius]
834      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu]
835      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ppress ! pressure             [dBar]
836      REAL(wp), DIMENSION(:,:), POINTER           ::   pti    ! in-situ temperature [Celcius]
837!      REAL(wp) :: fsatg
838!      REAL(wp) :: pfps, pfpt, pfphp
839      REAL(wp) :: zt, zs, zp, zh, zq, zxk
840      INTEGER  :: ji, jj            ! dummy loop indices
841      !
842      CALL wrk_alloc( jpi,jpj, pti  )
843      !
844      DO jj=1,jpj
845         DO ji=1,jpi
846            zh = ppress(ji,jj)
847! Theta1
848            zt = ptem(ji,jj)
849            zs = psal(ji,jj)
850            zp = 0.0
851            zxk= zh * fsatg( zs, zt, zp )
852            zt = zt + 0.5 * zxk
853            zq = zxk
854! Theta2
855            zp = zp + 0.5 * zh
856            zxk= zh*fsatg( zs, zt, zp )
857            zt = zt + 0.29289322 * ( zxk - zq )
858            zq = 0.58578644 * zxk + 0.121320344 * zq
859! Theta3
860            zxk= zh * fsatg( zs, zt, zp )
861            zt = zt + 1.707106781 * ( zxk - zq )
862            zq = 3.414213562 * zxk - 4.121320344 * zq
863! Theta4
864            zp = zp + 0.5 * zh
865            zxk= zh * fsatg( zs, zt, zp )
866            pti(ji,jj) = zt + ( zxk - 2.0 * zq ) / 6.0
867         END DO
868      END DO
869      !
870      CALL wrk_dealloc( jpi,jpj, pti  )
871      !
872   END FUNCTION tinsitu
873   !
874   FUNCTION fsatg( pfps, pfpt, pfphp )
875      !!----------------------------------------------------------------------
876      !!                 ***  FUNCTION fsatg  ***
877      !!
878      !! ** Purpose    :   Compute the Adiabatic laspse rate [Celcius].[decibar]^-1
879      !!
880      !! ** Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
881      !!
882      !! ** units      :   pressure        pfphp    decibars
883      !!                   temperature     pfpt     deg celsius (ipts-68)
884      !!                   salinity        pfps     (ipss-78)
885      !!                   adiabatic       fsatg    deg. c/decibar
886      !!----------------------------------------------------------------------
887      REAL(wp) :: pfps, pfpt, pfphp 
888      REAL(wp) :: fsatg
889      !
890      fsatg = (((-2.1687e-16*pfpt+1.8676e-14)*pfpt-4.6206e-13)*pfphp         &
891        &    +((2.7759e-12*pfpt-1.1351e-10)*(pfps-35.)+((-5.4481e-14*pfpt    &
892        &    +8.733e-12)*pfpt-6.7795e-10)*pfpt+1.8741e-8))*pfphp             &
893        &    +(-4.2393e-8*pfpt+1.8932e-6)*(pfps-35.)                         &
894        &    +((6.6228e-10*pfpt-6.836e-8)*pfpt+8.5258e-6)*pfpt+3.5803e-5
895      !
896    END FUNCTION fsatg
897    !!======================================================================
898END MODULE sbcisf
Note: See TracBrowser for help on using the repository browser.