New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcisf.F90 in branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

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

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

ISF cleaning branch: cleaning sbcisf + bug correction in zpshde_isf (ssumask instead of umask(:,:,1))

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