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

source: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 6006

Last change on this file since 6006 was 6006, checked in by mathiot, 8 years ago

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