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

source: branches/UKMO/r5518_amm15_test/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 7144

Last change on this file since 7144 was 7144, checked in by jcastill, 7 years ago

Remove svn keywords

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