source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 4666

Last change on this file since 4666 was 4666, checked in by mathiot, 6 years ago

#1331 : add ISOMIP config files + ice shelf code

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