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

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

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

Last change on this file since 7784 was 7784, checked in by cetlod, 7 years ago

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