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

source: branches/NERC/dev_r5589_is_oce_cpl/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 5802

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

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