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

Last change on this file since 7494 was 7494, checked in by timgraham, 4 years ago

Addition of extra diagnostic outputs in CMIP6 data request. NB. Will require update to field_def file from shaconemo

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