source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 5541

Last change on this file since 5541 was 5541, checked in by jchanut, 5 years ago

change eos_fzp function into a subroutine (needed for AGRIF), #1530

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