source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 9406

Last change on this file since 9406 was 9406, checked in by mathiot, 2 years ago

fix for ticket #2066 in nemo 3.6

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