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/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 8046

Last change on this file since 8046 was 8046, checked in by davestorkey, 7 years ago

UKMO/dev_r5518_GO6_package branch: Add extra option for freshwater input from ice shelves in coupled models.

File size: 49.0 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   USE lib_fortran, ONLY: glob_sum
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC   sbc_isf, sbc_isf_div, sbc_isf_alloc  ! routine called in sbcmod and divcur
34
35   ! public in order to be able to output then
36
37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   risf_tsc_b, risf_tsc   
38   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qisf              !: net heat flux from ice shelf
39   REAL(wp), PUBLIC ::   rn_hisf_tbl                 !: thickness of top boundary layer [m]
40   LOGICAL , PUBLIC ::   ln_divisf                   !: flag to correct divergence
41   INTEGER , PUBLIC ::   nn_isfblk                   !:
42   INTEGER , PUBLIC ::   nn_gammablk                 !:
43   LOGICAL , PUBLIC ::   ln_conserve                 !:
44   REAL(wp), PUBLIC ::   rn_gammat0                  !: temperature exchange coeficient
45   REAL(wp), PUBLIC ::   rn_gammas0                  !: salinity    exchange coeficient
46   REAL(wp), PUBLIC ::   rdivisf                     !: flag to test if fwf apply on divergence
47
48   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rzisf_tbl              !:depth of calving front (shallowest point) nn_isf ==2/3
49   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  rhisf_tbl, rhisf_tbl_0 !:thickness of tbl
50   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  r1_hisf_tbl            !:1/thickness of tbl
51   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ralpha                 !:proportion of bottom cell influenced by tbl
52   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  risfLeff               !:effective length (Leff) BG03 nn_isf==2
53   REAL(wp)   , PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  ttbl, stbl, utbl, vtbl !:top boundary layer variable at T point
54   INTEGER,    PUBLIC, ALLOCATABLE, SAVE, DIMENSION (:,:)     ::  misfkt, misfkb         !:Level of ice shelf base
55
56
57   REAL(wp), PUBLIC, SAVE ::   rcpi   = 2000.0_wp     ! phycst ?
58   REAL(wp), PUBLIC, SAVE ::   kappa  = 1.54e-6_wp    ! phycst ?
59   REAL(wp), PUBLIC, SAVE ::   rhoisf = 920.0_wp      ! phycst ?
60   REAL(wp), PUBLIC, SAVE ::   tsurf  = -20.0_wp      ! phycst ?
61   REAL(wp), PUBLIC, SAVE ::   lfusisf= 0.334e6_wp    ! phycst ?
62
63!: Variable used in fldread to read the forcing file (nn_isf == 4 .OR. nn_isf == 3)
64   CHARACTER(len=100), PUBLIC ::   cn_dirisf  = './'    !: Root directory for location of ssr files
65   TYPE(FLD_N)       , PUBLIC ::   sn_qisf, sn_fwfisf     !: information about the runoff file to be read
66   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_qisf, sf_fwfisf
67   TYPE(FLD_N)       , PUBLIC ::   sn_rnfisf              !: information about the runoff file to be read
68   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnfisf           
69   TYPE(FLD_N)       , PUBLIC ::   sn_depmax_isf, sn_depmin_isf, sn_Leff_isf     !: information about the runoff file to be read
70   
71   !! * Substitutions
72#  include "domzgr_substitute.h90"
73   !!----------------------------------------------------------------------
74   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
75   !! $Id$
76   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
77   !!----------------------------------------------------------------------
78
79CONTAINS
80 
81  SUBROUTINE sbc_isf(kt)
82    INTEGER, INTENT(in)          ::   kt         ! ocean time step
83    INTEGER                      ::   ji, jj, jk, ijkmin, inum, ierror
84    INTEGER                      ::   ikt, ikb   ! top and bottom level of the isf boundary layer
85    REAL(wp)                     ::   zgreenland_fwfisf_sum, zantarctica_fwfisf_sum
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
256            IF( lk_oasis) THEN
257            ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
258            IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
259
260              ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern
261              ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets
262              ! to preserve total freshwater conservation in coupled models without an active ice sheet model.
263 
264              ! All related global sums must be done bit reproducibly
265               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
266
267               ! use ABS function because we need to preserve the sign of fwfisf
268               WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  &
269              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) &
270              &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) )
271
272               ! check
273               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
274
275               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
276
277               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum
278
279               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
280
281               ! use ABS function because we need to preserve the sign of fwfisf
282               WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
283              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) &
284              &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) )
285     
286               ! check
287               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
288
289               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
290
291               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
292
293            ENDIF
294            ENDIF
295
296            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
297            stbl(:,:)   = soce
298
299         ELSE IF (nn_isf == 4) THEN
300            ! specified fwf and heat flux forcing beneath the ice shelf
301            CALL fld_read ( kt, nn_fsbc, sf_fwfisf   )
302            !CALL fld_read ( kt, nn_fsbc, sf_qisf   )
303            fwfisf(:,:) = sf_fwfisf(1)%fnow(:,:,1)            ! fwf
304
305            IF( lk_oasis) THEN
306            ! nn_coupled_iceshelf_fluxes uninitialised unless lk_oasis=true
307            IF( nn_coupled_iceshelf_fluxes .gt. 0 ) THEN
308
309              ! Adjust total iceshelf melt rates so that sum of iceberg calving and iceshelf melting in the northern
310              ! and southern hemispheres equals rate of increase of mass of greenland and antarctic ice sheets
311              ! to preserve total freshwater conservation in coupled models without an active ice sheet model.
312
313              ! All related global sums must be done bit reproducibly
314               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
315
316               ! use ABS function because we need to preserve the sign of fwfisf
317               WHERE( greenland_icesheet_mask(:,:) == 1.0 )                                                                  &
318              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( greenland_icesheet_mass_rate_of_change * (1.0-rn_greenland_calving_fraction) &
319              &                           / ( zgreenland_fwfisf_sum + 1.0e-10_wp ) )
320
321               ! check
322               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting climatology (kg/s) : ',zgreenland_fwfisf_sum
323
324               zgreenland_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * greenland_icesheet_mask(:,:) )
325
326               IF(lwp) WRITE(numout, *) 'Greenland iceshelf melting adjusted value (kg/s) : ',zgreenland_fwfisf_sum
327
328               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
329
330               ! use ABS function because we need to preserve the sign of fwfisf
331               WHERE( antarctica_icesheet_mask(:,:) == 1.0 ) &
332              &    fwfisf(:,:) = fwfisf(:,:)  * ABS( antarctica_icesheet_mass_rate_of_change * (1.0-rn_antarctica_calving_fraction) &
333              &                           / ( zantarctica_fwfisf_sum + 1.0e-10_wp ) )
334     
335               ! check
336               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting climatology (kg/s) : ',zantarctica_fwfisf_sum
337
338               zantarctica_fwfisf_sum = glob_sum( fwfisf(:,:) * e1t(:,:) * e2t(:,:) * antarctica_icesheet_mask(:,:) )
339
340               IF(lwp) WRITE(numout, *) 'Antarctica iceshelf melting adjusted value (kg/s) : ',zantarctica_fwfisf_sum
341
342            ENDIF
343            ENDIF
344
345            qisf(:,:)   = fwfisf(:,:) * lfusisf              ! heat        flux
346            !qisf(:,:)   = sf_qisf(1)%fnow(:,:,1)              ! heat flux
347            stbl(:,:)   = soce
348
349         END IF
350         ! compute tsc due to isf
351         ! WARNING water add at temp = 0C, correction term is added, maybe better here but need a 3D variable).
352!         zpress = grav*rau0*fsdept(ji,jj,jk)*1.e-04
353         zt_frz = -1.9 !eos_fzp( tsn(ji,jj,jk,jp_sal), zpress )
354         risf_tsc(:,:,jp_tem) = qisf(:,:) * r1_rau0_rcp - rdivisf * fwfisf(:,:) * zt_frz * r1_rau0 !
355         
356         ! salt effect already take into account in vertical advection
357         risf_tsc(:,:,jp_sal) = (1.0_wp-rdivisf) * fwfisf(:,:) * stbl(:,:) * r1_rau0
358
359         ! output
360         IF( iom_use('qlatisf' ) )   CALL iom_put('qlatisf', qisf)
361         IF( iom_use('fwfisf'  ) )   CALL iom_put('fwfisf' , fwfisf * stbl(:,:) / soce )
362
363         ! if apply only on the trend and not as a volume flux (rdivisf = 0), fwfisf have to be set to 0 now
364         fwfisf(:,:) = rdivisf * fwfisf(:,:)         
365 
366         ! lbclnk
367         CALL lbc_lnk(risf_tsc(:,:,jp_tem),'T',1.)
368         CALL lbc_lnk(risf_tsc(:,:,jp_sal),'T',1.)
369         CALL lbc_lnk(fwfisf(:,:)   ,'T',1.)
370         CALL lbc_lnk(qisf(:,:)     ,'T',1.)
371
372!=============================================================================================================================================
373         IF ( iom_use('fwfisf3d') .OR. iom_use('qlatisf3d') .OR. iom_use('qhcisf3d') .OR. iom_use('qhcisf')) THEN
374            CALL wrk_alloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d )
375            CALL wrk_alloc( jpi,jpj,     zqhcisf2d                        )
376
377            zfwfisf3d(:,:,:) = 0.0_wp                         ! 3d ice shelf melting (kg/m2/s)
378            zqhcisf3d(:,:,:) = 0.0_wp                         ! 3d heat content flux (W/m2)
379            zqlatisf3d(:,:,:)= 0.0_wp                         ! 3d ice shelf melting latent heat flux (W/m2)
380            zqhcisf2d(:,:)   = fwfisf(:,:) * zt_frz * rcp     ! 2d heat content flux (W/m2)
381
382            DO jj = 1,jpj
383               DO ji = 1,jpi
384                  ikt = misfkt(ji,jj)
385                  ikb = misfkb(ji,jj)
386                  DO jk = ikt, ikb - 1
387                     zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
388                     zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
389                     zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * fse3t(ji,jj,jk)
390                  END DO
391                  zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf   (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk)
392                  zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk)
393                  zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf     (ji,jj) * r1_hisf_tbl(ji,jj) * ralpha(ji,jj) * fse3t(ji,jj,jk)
394               END DO
395            END DO
396
397            CALL iom_put('fwfisf3d' , zfwfisf3d (:,:,:))
398            CALL iom_put('qlatisf3d', zqlatisf3d(:,:,:))
399            CALL iom_put('qhcisf3d' , zqhcisf3d (:,:,:))
400            CALL iom_put('qhcisf'   , zqhcisf2d (:,:  ))
401
402            CALL wrk_dealloc( jpi,jpj,jpk, zfwfisf3d, zqhcisf3d, zqlatisf3d )
403            CALL wrk_dealloc( jpi,jpj,     zqhcisf2d                        )
404         END IF
405!=============================================================================================================================================
406
407         IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
408            IF( ln_rstart .AND.    &                     ! Restart: read in restart file
409                 & iom_varid( numror, 'fwf_isf_b', ldstop = .FALSE. ) > 0 ) THEN
410               IF(lwp) WRITE(numout,*) '          nit000-1 isf tracer content forcing fields read in the restart file'
411               CALL iom_get( numror, jpdom_autoglo, 'fwf_isf_b', fwfisf_b(:,:) )   ! before salt content isf_tsc trend
412               CALL iom_get( numror, jpdom_autoglo, 'isf_sc_b', risf_tsc_b(:,:,jp_sal) )   ! before salt content isf_tsc trend
413               CALL iom_get( numror, jpdom_autoglo, 'isf_hc_b', risf_tsc_b(:,:,jp_tem) )   ! before salt content isf_tsc trend
414            ELSE
415               fwfisf_b(:,:)    = fwfisf(:,:)
416               risf_tsc_b(:,:,:)= risf_tsc(:,:,:)
417            END IF
418         ENDIF
419         !
420      END IF
421 
422  END SUBROUTINE sbc_isf
423
424  INTEGER FUNCTION sbc_isf_alloc()
425      !!----------------------------------------------------------------------
426      !!               ***  FUNCTION sbc_isf_rnf_alloc  ***
427      !!----------------------------------------------------------------------
428      sbc_isf_alloc = 0       ! set to zero if no array to be allocated
429      IF( .NOT. ALLOCATED( qisf ) ) THEN
430         ALLOCATE(  risf_tsc(jpi,jpj,jpts), risf_tsc_b(jpi,jpj,jpts), qisf(jpi,jpj)   , &
431               &    rhisf_tbl(jpi,jpj)    , r1_hisf_tbl(jpi,jpj), rzisf_tbl(jpi,jpj)  , &
432               &    ttbl(jpi,jpj)         , stbl(jpi,jpj)       , utbl(jpi,jpj)       , &
433               &    vtbl(jpi, jpj)        , risfLeff(jpi,jpj)   , rhisf_tbl_0(jpi,jpj), &
434               &    ralpha(jpi,jpj)       , misfkt(jpi,jpj)     , misfkb(jpi,jpj)     , &
435               &    STAT= sbc_isf_alloc )
436         !
437         IF( lk_mpp                  )   CALL mpp_sum ( sbc_isf_alloc )
438         IF( sbc_isf_alloc /= 0 )   CALL ctl_warn('sbc_isf_alloc: failed to allocate arrays.')
439         !
440      ENDIF
441  END FUNCTION
442
443  SUBROUTINE sbc_isf_bg03(kt)
444   !!==========================================================================
445   !!                 *** SUBROUTINE sbcisf_bg03  ***
446   !! add net heat and fresh water flux from ice shelf melting
447   !! into the adjacent ocean using the parameterisation by
448   !! Beckmann and Goosse (2003), "A parameterization of ice shelf-ocean
449   !!     interaction for climate models", Ocean Modelling 5(2003) 157-170.
450   !!  (hereafter BG)
451   !!==========================================================================
452   !!----------------------------------------------------------------------
453   !!   sbc_isf_bg03      : routine called from sbcmod
454   !!----------------------------------------------------------------------
455   !!
456   !! ** Purpose   :   Add heat and fresh water fluxes due to ice shelf melting
457   !! ** Reference :   Beckmann et Goosse, 2003, Ocean Modelling
458   !!
459   !! History :
460   !!      !  06-02  (C. Wang) Original code
461   !!----------------------------------------------------------------------
462
463    INTEGER, INTENT ( in ) :: kt
464
465    INTEGER :: ji, jj, jk, jish  !temporary integer
466    INTEGER :: ijkmin
467    INTEGER :: ii, ij, ik 
468    INTEGER :: inum
469
470    REAL(wp) :: zt_sum      ! sum of the temperature between 200m and 600m
471    REAL(wp) :: zt_ave      ! averaged temperature between 200m and 600m
472    REAL(wp) :: zt_frz      ! freezing point temperature at depth z
473    REAL(wp) :: zpress      ! pressure to compute the freezing point in depth
474   
475    !!----------------------------------------------------------------------
476    IF ( nn_timing == 1 ) CALL timing_start('sbc_isf_bg03')
477     !
478
479    ! This test is false only in the very first time step of a run (JMM ???- Initialy build to skip 1rst year of run )
480    DO ji = 1, jpi
481       DO jj = 1, jpj
482          ik = misfkt(ji,jj)
483          !! Initialize arrays to 0 (each step)
484          zt_sum = 0.e0_wp
485          IF ( ik .GT. 1 ) THEN
486    ! 3. -----------the average temperature between 200m and 600m ---------------------
487             DO jk = misfkt(ji,jj),misfkb(ji,jj)
488             ! freezing point temperature  at ice shelf base BG eq. 2 (JMM sign pb ??? +7.64e-4 !!!)
489             ! after verif with UNESCO, wrong sign in BG eq. 2
490             ! Calculate freezing temperature
491                zpress = grav*rau0*fsdept(ji,jj,ik)*1.e-04 
492                CALL eos_fzp(tsb(ji,jj,ik,jp_sal), zt_frz, zpress) 
493                zt_sum = zt_sum + (tsn(ji,jj,ik,jp_tem)-zt_frz) * fse3t(ji,jj,ik) * tmask(ji,jj,ik)  ! sum temp
494             ENDDO
495             zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value
496   
497    ! 4. ------------Net heat flux and fresh water flux due to the ice shelf
498          ! For those corresponding to zonal boundary   
499             qisf(ji,jj) = - rau0 * rcp * rn_gammat0 * risfLeff(ji,jj) * e1t(ji,jj) * zt_ave  &
500                         & / (e1t(ji,jj) * e2t(ji,jj)) * tmask(ji,jj,ik) 
501             
502             fwfisf(ji,jj) = qisf(ji,jj) / lfusisf          !fresh water flux kg/(m2s)                 
503             fwfisf(ji,jj) = fwfisf(ji,jj) * ( soce / stbl(ji,jj) )
504             !add to salinity trend
505          ELSE
506             qisf(ji,jj) = 0._wp ; fwfisf(ji,jj) = 0._wp
507          END IF
508       ENDDO
509    ENDDO
510    !
511    IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_bg03')
512  END SUBROUTINE sbc_isf_bg03
513
514   SUBROUTINE sbc_isf_cav( kt )
515      !!---------------------------------------------------------------------
516      !!                     ***  ROUTINE sbc_isf_cav  ***
517      !!
518      !! ** Purpose :   handle surface boundary condition under ice shelf
519      !!
520      !! ** Method  : -
521      !!
522      !! ** Action  :   utau, vtau : remain unchanged
523      !!                taum, wndm : remain unchanged
524      !!                qns        : update heat flux below ice shelf
525      !!                emp, emps  : update freshwater flux below ice shelf
526      !!---------------------------------------------------------------------
527      INTEGER, INTENT(in)          ::   kt         ! ocean time step
528      !
529      LOGICAL :: ln_isomip = .true.
530      REAL(wp), DIMENSION(:,:), POINTER       ::   zfrz,zpress,zti
531      REAL(wp), DIMENSION(:,:), POINTER       ::   zgammat2d, zgammas2d 
532      !REAL(wp), DIMENSION(:,:), POINTER ::   zqisf, zfwfisf
533      REAL(wp) ::   zlamb1, zlamb2, zlamb3
534      REAL(wp) ::   zeps1,zeps2,zeps3,zeps4,zeps6,zeps7
535      REAL(wp) ::   zaqe,zbqe,zcqe,zaqer,zdis,zsfrz,zcfac
536      REAL(wp) ::   zfwflx, zhtflx, zhtflx_b
537      REAL(wp) ::   zgammat, zgammas
538      REAL(wp) ::   zeps   =  -1.e-20_wp        !==   Local constant initialization   ==!
539      INTEGER  ::   ji, jj     ! dummy loop indices
540      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integers
541      INTEGER  ::   ierror     ! return error code
542      LOGICAL  ::   lit=.TRUE.
543      INTEGER  ::   nit
544      !!---------------------------------------------------------------------
545      !
546      ! coeficient for linearisation of tfreez
547      zlamb1=-0.0575
548      zlamb2=0.0901
549      zlamb3=-7.61e-04
550      IF( nn_timing == 1 )  CALL timing_start('sbc_isf_cav')
551      !
552      CALL wrk_alloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
553
554      zcfac=0.0_wp 
555      IF (ln_conserve)  zcfac=1.0_wp
556      zpress(:,:)=0.0_wp
557      zgammat2d(:,:)=0.0_wp
558      zgammas2d(:,:)=0.0_wp
559      !
560      !
561!CDIR COLLAPSE
562      DO jj = 1, jpj
563         DO ji = 1, jpi
564            ! Crude approximation for pressure (but commonly used)
565            ! 1e-04 to convert from Pa to dBar
566            zpress(ji,jj)=grav*rau0*fsdepw(ji,jj,mikt(ji,jj))*1.e-04
567            !
568         END DO
569      END DO
570
571! Calculate in-situ temperature (ref to surface)
572      zti(:,:)=tinsitu( ttbl, stbl, zpress )
573! Calculate freezing temperature
574      CALL eos_fzp( sss_m(:,:), zfrz(:,:), zpress )
575
576     
577      zhtflx=0._wp ; zfwflx=0._wp
578      IF (nn_isfblk == 1) THEN
579         DO jj = 1, jpj
580            DO ji = 1, jpi
581               IF (mikt(ji,jj) > 1 ) THEN
582                  nit = 1; lit = .TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp
583                  DO WHILE ( lit )
584! compute gamma
585                     CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit)
586! zhtflx is upward heat flux (out of ocean)
587                     zhtflx = zgammat*rcp*rau0*(zti(ji,jj)-zfrz(ji,jj))
588! zwflx is upward water flux
589                     zfwflx = - zhtflx/lfusisf
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. 100) CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
595
596! save gammat and compute zhtflx_b
597                     zgammat2d(ji,jj)=zgammat
598                     zhtflx_b = zhtflx
599                  END DO
600
601                  qisf(ji,jj) = - zhtflx
602! For genuine ISOMIP protocol this should probably be something like
603                  fwfisf(ji,jj) = zfwflx  * ( soce / MAX(stbl(ji,jj),zeps))
604               ELSE
605                  fwfisf(ji,jj) = 0._wp
606                  qisf(ji,jj)   = 0._wp
607               END IF
608            !
609            END DO
610         END DO
611
612      ELSE IF (nn_isfblk == 2 ) THEN
613
614! More complicated 3 equation thermodynamics as in MITgcm
615!CDIR COLLAPSE
616         DO jj = 2, jpj
617            DO ji = 2, jpi
618               IF (mikt(ji,jj) > 1 ) THEN
619                  nit=1; lit=.TRUE.; zgammat=rn_gammat0; zgammas=rn_gammas0; zhtflx_b=0._wp; zhtflx=0._wp
620                  DO WHILE ( lit )
621                     CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, ji, jj, lit)
622
623                     zeps1=rcp*rau0*zgammat
624                     zeps2=lfusisf*rau0*zgammas
625                     zeps3=rhoisf*rcpi*kappa/risfdep(ji,jj)
626                     zeps4=zlamb2+zlamb3*risfdep(ji,jj)
627                     zeps6=zeps4-zti(ji,jj)
628                     zeps7=zeps4-tsurf
629                     zaqe=zlamb1 * (zeps1 + zeps3)
630                     zaqer=0.5/zaqe
631                     zbqe=zeps1*zeps6+zeps3*zeps7-zeps2
632                     zcqe=zeps2*stbl(ji,jj)
633                     zdis=zbqe*zbqe-4.0*zaqe*zcqe               
634! Presumably zdis can never be negative because gammas is very small compared to gammat
635                     zsfrz=(-zbqe-SQRT(zdis))*zaqer
636                     IF (zsfrz .lt. 0.0) zsfrz=(-zbqe+SQRT(zdis))*zaqer
637                     zfrz(ji,jj)=zeps4+zlamb1*zsfrz
638 
639! zfwflx is upward water flux
640                     zfwflx= rau0 * zgammas * ( (zsfrz-stbl(ji,jj)) / zsfrz )
641! zhtflx is upward heat flux (out of ocean)
642! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
643                     zhtflx = ( zgammat*rau0 - zcfac*zfwflx ) * rcp * (zti(ji,jj) - zfrz(ji,jj) ) 
644! zwflx is upward water flux
645! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
646                     zfwflx = ( zgammas*rau0 - zcfac*zfwflx ) * (zsfrz - stbl(ji,jj)) / stbl(ji,jj)
647! test convergence and compute gammat
648                     IF (( zhtflx - zhtflx_b) .LE. 0.01 ) lit = .FALSE.
649
650                     nit = nit + 1
651                     IF (nit .GE. 51) THEN
652                        WRITE(numout,*) "sbcisf : too many iteration ... ", &
653                            &  zhtflx, zhtflx_b, zgammat, zgammas, nn_gammablk, ji, jj, mikt(ji,jj), narea
654                        CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' )
655                     END IF
656! save gammat and compute zhtflx_b
657                     zgammat2d(ji,jj)=zgammat
658                     zgammas2d(ji,jj)=zgammas
659                     zhtflx_b = zhtflx
660
661                  END DO
662! If non conservative we have zcfac=0.0 so zhtflx is as ISOMIP but with different zfrz value
663                  qisf(ji,jj) = - zhtflx 
664! If non conservative we have zcfac=0.0 so what follows is then zfwflx*sss_m/zsfrz
665                  fwfisf(ji,jj) = zfwflx 
666               ELSE
667                  fwfisf(ji,jj) = 0._wp
668                  qisf(ji,jj)   = 0._wp
669               ENDIF
670               !
671            END DO
672         END DO
673      ENDIF
674      ! lbclnk
675      CALL lbc_lnk(zgammas2d(:,:),'T',1.)
676      CALL lbc_lnk(zgammat2d(:,:),'T',1.)
677      ! output
678      CALL iom_put('isfgammat', zgammat2d)
679      CALL iom_put('isfgammas', zgammas2d)
680      !
681      CALL wrk_dealloc( jpi,jpj, zfrz,zpress,zti, zgammat2d, zgammas2d )
682      !
683      IF( nn_timing == 1 )  CALL timing_stop('sbc_isf_cav')
684
685   END SUBROUTINE sbc_isf_cav
686
687   SUBROUTINE sbc_isf_gammats(gt, gs, zqhisf, zqwisf, ji, jj, lit )
688      !!----------------------------------------------------------------------
689      !! ** Purpose    : compute the coefficient echange for heat flux
690      !!
691      !! ** Method     : gamma assume constant or depends of u* and stability
692      !!
693      !! ** References : Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
694      !!                Jenkins et al., 2010, JPO, p2298-2312
695      !!---------------------------------------------------------------------
696      REAL(wp), INTENT(inout) :: gt, gs, zqhisf, zqwisf
697      INTEGER , INTENT(in)    :: ji,jj
698      LOGICAL , INTENT(inout) :: lit
699
700      INTEGER  :: ikt                 ! loop index
701      REAL(wp) :: zut, zvt, zustar           ! U, V at T point and friction velocity
702      REAL(wp) :: zdku, zdkv                 ! U, V shear
703      REAL(wp) :: zPr, zSc, zRc              ! Prandtl, Scmidth and Richardson number
704      REAL(wp) :: zmob, zmols                ! Monin Obukov length, coriolis factor at T point
705      REAL(wp) :: zbuofdep, zhnu             ! Bouyancy length scale, sublayer tickness
706      REAL(wp) :: zhmax                      ! limitation of mol
707      REAL(wp) :: zetastar                   ! stability parameter
708      REAL(wp) :: zgmolet, zgmoles, zgturb   ! contribution of modelecular sublayer and turbulence
709      REAL(wp) :: zcoef                      ! temporary coef
710      REAL(wp) :: zdep
711      REAL(wp), PARAMETER :: zxsiN = 0.052   ! dimensionless constant
712      REAL(wp), PARAMETER :: epsln = 1.0e-20 ! a small positive number
713      REAL(wp), PARAMETER :: znu   = 1.95e-6 ! kinamatic viscosity of sea water (m2.s-1)
714      REAL(wp) ::   rcs      = 1.0e-3_wp        ! conversion: mm/s ==> m/s
715      REAL(wp), DIMENSION(2) :: zts, zab
716      !!---------------------------------------------------------------------
717      !
718      IF( nn_gammablk == 0 ) THEN
719      !! gamma is constant (specified in namelist)
720         gt = rn_gammat0
721         gs = rn_gammas0
722         lit = .FALSE.
723      ELSE IF ( nn_gammablk == 1 ) THEN
724      !! gamma is assume to be proportional to u*
725      !! WARNING in case of Losh 2008 tbl parametrization,
726      !! you have to used the mean value of u in the boundary layer)
727      !! not yet coded
728      !! Jenkins et al., 2010, JPO, p2298-2312
729         ikt = mikt(ji,jj)
730      !! Compute U and V at T points
731   !      zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
732   !      zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
733          zut = utbl(ji,jj)
734          zvt = vtbl(ji,jj)
735
736      !! compute ustar
737         zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
738      !! Compute mean value over the TBL
739
740      !! Compute gammats
741         gt = zustar * rn_gammat0
742         gs = zustar * rn_gammas0
743         lit = .FALSE.
744      ELSE IF ( nn_gammablk == 2 ) THEN
745      !! gamma depends of stability of boundary layer
746      !! WARNING in case of Losh 2008 tbl parametrization,
747      !! you have to used the mean value of u in the boundary layer)
748      !! not yet coded
749      !! Holland and Jenkins, 1999, JPO, p1787-1800, eq 14
750      !! as MOL depends of flux and flux depends of MOL, best will be iteration (TO DO)
751               ikt = mikt(ji,jj)
752
753      !! Compute U and V at T points
754               zut = 0.5 * ( utbl(ji-1,jj  ) + utbl(ji,jj) )
755               zvt = 0.5 * ( vtbl(ji  ,jj-1) + vtbl(ji,jj) )
756
757      !! compute ustar
758               zustar = SQRT( rn_tfri2 * (zut * zut + zvt * zvt) )
759               IF (zustar == 0._wp) THEN           ! only for kt = 1 I think
760                 gt = rn_gammat0
761                 gs = rn_gammas0
762               ELSE
763      !! compute Rc number (as done in zdfric.F90)
764               zcoef = 0.5 / fse3w(ji,jj,ikt)
765               !                                            ! shear of horizontal velocity
766               zdku = zcoef * (  un(ji-1,jj  ,ikt  ) + un(ji,jj,ikt  )   &
767                  &             -un(ji-1,jj  ,ikt+1) - un(ji,jj,ikt+1)  )
768               zdkv = zcoef * (  vn(ji  ,jj-1,ikt  ) + vn(ji,jj,ikt  )   &
769                  &             -vn(ji  ,jj-1,ikt+1) - vn(ji,jj,ikt+1)  )
770               !                                            ! richardson number (minimum value set to zero)
771               zRc = rn2(ji,jj,ikt+1) / ( zdku*zdku + zdkv*zdkv + 1.e-20 )
772
773      !! compute Pr and Sc number (can be improved)
774               zPr =   13.8
775               zSc = 2432.0
776
777      !! compute gamma mole
778               zgmolet = 12.5 * zPr ** (2.0/3.0) - 6.0
779               zgmoles = 12.5 * zSc ** (2.0/3.0) -6.0
780
781      !! compute bouyancy
782               zts(jp_tem) = ttbl(ji,jj)
783               zts(jp_sal) = stbl(ji,jj)
784               zdep        = fsdepw(ji,jj,ikt)
785               !
786               CALL eos_rab( zts, zdep, zab )
787                  !
788      !! compute length scale
789               zbuofdep = grav * ( zab(jp_tem) * zqhisf - zab(jp_sal) * zqwisf )  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
790
791      !! compute Monin Obukov Length
792               ! Maximum boundary layer depth
793               zhmax = fsdept(ji,jj,mbkt(ji,jj)) - fsdepw(ji,jj,mikt(ji,jj)) -0.001
794               ! Compute Monin obukhov length scale at the surface and Ekman depth:
795               zmob   = zustar ** 3 / (vkarmn * (zbuofdep + epsln))
796               zmols  = SIGN(1._wp, zmob) * MIN(ABS(zmob), zhmax) * tmask(ji,jj,ikt)
797
798      !! compute eta* (stability parameter)
799               zetastar = 1 / ( SQRT(1 + MAX(zxsiN * zustar / ( ABS(ff(ji,jj)) * zmols * zRc ), 0.0)))
800
801      !! compute the sublayer thickness
802               zhnu = 5 * znu / zustar
803      !! compute gamma turb
804               zgturb = 1/vkarmn * LOG(zustar * zxsiN * zetastar * zetastar / ( ABS(ff(ji,jj)) * zhnu )) &
805               &      + 1 / ( 2 * zxsiN * zetastar ) - 1/vkarmn
806
807      !! compute gammats
808               gt = zustar / (zgturb + zgmolet)
809               gs = zustar / (zgturb + zgmoles)
810               END IF
811      END IF
812
813   END SUBROUTINE
814
815   SUBROUTINE sbc_isf_tbl( varin, varout, cptin )
816      !!----------------------------------------------------------------------
817      !!                  ***  SUBROUTINE sbc_isf_tbl  ***
818      !!
819      !! ** Purpose : compute mean T/S/U/V in the boundary layer
820      !!
821      !!----------------------------------------------------------------------
822      REAL(wp), DIMENSION(:,:,:), INTENT(in) :: varin
823      REAL(wp), DIMENSION(:,:)  , INTENT(out):: varout
824     
825      CHARACTER(len=1), INTENT(in) :: cptin ! point of variable in/out
826
827      REAL(wp) :: ze3, zhk
828      REAL(wp), DIMENSION(:,:), POINTER :: zikt
829
830      INTEGER :: ji,jj,jk
831      INTEGER :: ikt,ikb
832      INTEGER, DIMENSION(:,:), POINTER :: mkt, mkb
833
834      CALL wrk_alloc( jpi,jpj, mkt, mkb  )
835      CALL wrk_alloc( jpi,jpj, zikt )
836
837      ! get first and last level of tbl
838      mkt(:,:) = misfkt(:,:)
839      mkb(:,:) = misfkb(:,:)
840
841      varout(:,:)=0._wp
842      DO jj = 2,jpj
843         DO ji = 2,jpi
844            IF (ssmask(ji,jj) == 1) THEN
845               ikt = mkt(ji,jj)
846               ikb = mkb(ji,jj)
847
848               ! level fully include in the ice shelf boundary layer
849               DO jk = ikt, ikb - 1
850                  ze3 = fse3t_n(ji,jj,jk)
851                  IF (cptin == 'T' ) varout(ji,jj) = varout(ji,jj) + varin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3
852                  IF (cptin == 'U' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji-1,jj,jk)) &
853                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
854                  IF (cptin == 'V' ) varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,jk) + varin(ji,jj-1,jk)) &
855                     &                                                       * r1_hisf_tbl(ji,jj) * ze3
856               END DO
857
858               ! level partially include in ice shelf boundary layer
859               zhk = SUM( fse3t_n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)
860               IF (cptin == 'T') &
861                   &  varout(ji,jj) = varout(ji,jj) + varin(ji,jj,ikb) * (1._wp - zhk)
862               IF (cptin == 'U') &
863                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji-1,jj,ikb)) * (1._wp - zhk)
864               IF (cptin == 'V') &
865                   &  varout(ji,jj) = varout(ji,jj) + 0.5_wp * (varin(ji,jj,ikb) + varin(ji,jj-1,ikb)) * (1._wp - zhk)
866            END IF
867         END DO
868      END DO
869
870      CALL wrk_dealloc( jpi,jpj, mkt, mkb )     
871      CALL wrk_dealloc( jpi,jpj, zikt ) 
872
873      IF (cptin == 'T') CALL lbc_lnk(varout,'T',1.)
874      IF (cptin == 'U' .OR. cptin == 'V') CALL lbc_lnk(varout,'T',-1.)
875
876   END SUBROUTINE sbc_isf_tbl
877     
878
879   SUBROUTINE sbc_isf_div( phdivn )
880      !!----------------------------------------------------------------------
881      !!                  ***  SUBROUTINE sbc_isf_div  ***
882      !!       
883      !! ** Purpose :   update the horizontal divergence with the runoff inflow
884      !!
885      !! ** Method  :   
886      !!                CAUTION : risf_tsc(:,:,jp_sal) is negative (outflow) increase the
887      !!                          divergence and expressed in m/s
888      !!
889      !! ** Action  :   phdivn   decreased by the runoff inflow
890      !!----------------------------------------------------------------------
891      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
892      !!
893      INTEGER  ::   ji, jj, jk   ! dummy loop indices
894      INTEGER  ::   ikt, ikb 
895      INTEGER  ::   nk_isf
896      REAL(wp)     ::   zhk, z1_hisf_tbl, zhisf_tbl
897      REAL(wp)     ::   zfact     ! local scalar
898      !!----------------------------------------------------------------------
899      !
900      zfact   = 0.5_wp
901      !
902      IF (lk_vvl) THEN     ! need to re compute level distribution of isf fresh water
903         DO jj = 1,jpj
904            DO ji = 1,jpi
905               ikt = misfkt(ji,jj)
906               ikb = misfkt(ji,jj)
907               ! thickness of boundary layer at least the top level thickness
908               rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), fse3t(ji,jj,ikt))
909
910               ! determine the deepest level influenced by the boundary layer
911               ! test on tmask useless ?????
912               DO jk = ikt, mbkt(ji,jj)
913                  IF ( (SUM(fse3t(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk
914               END DO
915               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(fse3t(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness.
916               misfkb(ji,jj) = ikb                                                  ! last wet level of the tbl
917               r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj)
918
919               zhk           = SUM( fse3t(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)  ! proportion of tbl cover by cell from ikt to ikb - 1
920               ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / fse3t(ji,jj,ikb)  ! proportion of bottom cell influenced by boundary layer
921            END DO
922         END DO
923      END IF ! vvl case
924      !
925      DO jj = 1,jpj
926         DO ji = 1,jpi
927               ikt = misfkt(ji,jj)
928               ikb = misfkb(ji,jj)
929               ! level fully include in the ice shelf boundary layer
930               DO jk = ikt, ikb - 1
931                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) + ( fwfisf(ji,jj) + fwfisf_b(ji,jj) ) &
932                    &               * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact
933               END DO
934               ! level partially include in ice shelf boundary layer
935               phdivn(ji,jj,ikb) = phdivn(ji,jj,ikb) + ( fwfisf(ji,jj) &
936                  &             + fwfisf_b(ji,jj) ) * r1_hisf_tbl(ji,jj) * r1_rau0 * zfact * ralpha(ji,jj) 
937            !==   ice shelf melting mass distributed over several levels   ==!
938         END DO
939      END DO
940      !
941   END SUBROUTINE sbc_isf_div
942                       
943   FUNCTION tinsitu( ptem, psal, ppress ) RESULT( pti )
944      !!----------------------------------------------------------------------
945      !!                 ***  ROUTINE eos_init  ***
946      !!
947      !! ** Purpose :   Compute the in-situ temperature [Celcius]
948      !!
949      !! ** Method  :   
950      !!
951      !! Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
952      !!----------------------------------------------------------------------
953      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptem   ! potential temperature [Celcius]
954      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   psal   ! salinity             [psu]
955      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ppress ! pressure             [dBar]
956      REAL(wp), DIMENSION(:,:), POINTER           ::   pti    ! in-situ temperature [Celcius]
957!      REAL(wp) :: fsatg
958!      REAL(wp) :: pfps, pfpt, pfphp
959      REAL(wp) :: zt, zs, zp, zh, zq, zxk
960      INTEGER  :: ji, jj            ! dummy loop indices
961      !
962      CALL wrk_alloc( jpi,jpj, pti  )
963      !
964      DO jj=1,jpj
965         DO ji=1,jpi
966            zh = ppress(ji,jj)
967! Theta1
968            zt = ptem(ji,jj)
969            zs = psal(ji,jj)
970            zp = 0.0
971            zxk= zh * fsatg( zs, zt, zp )
972            zt = zt + 0.5 * zxk
973            zq = zxk
974! Theta2
975            zp = zp + 0.5 * zh
976            zxk= zh*fsatg( zs, zt, zp )
977            zt = zt + 0.29289322 * ( zxk - zq )
978            zq = 0.58578644 * zxk + 0.121320344 * zq
979! Theta3
980            zxk= zh * fsatg( zs, zt, zp )
981            zt = zt + 1.707106781 * ( zxk - zq )
982            zq = 3.414213562 * zxk - 4.121320344 * zq
983! Theta4
984            zp = zp + 0.5 * zh
985            zxk= zh * fsatg( zs, zt, zp )
986            pti(ji,jj) = zt + ( zxk - 2.0 * zq ) / 6.0
987         END DO
988      END DO
989      !
990      CALL wrk_dealloc( jpi,jpj, pti  )
991      !
992   END FUNCTION tinsitu
993   !
994   FUNCTION fsatg( pfps, pfpt, pfphp )
995      !!----------------------------------------------------------------------
996      !!                 ***  FUNCTION fsatg  ***
997      !!
998      !! ** Purpose    :   Compute the Adiabatic laspse rate [Celcius].[decibar]^-1
999      !!
1000      !! ** Reference  :   Bryden,h.,1973,deep-sea res.,20,401-408
1001      !!
1002      !! ** units      :   pressure        pfphp    decibars
1003      !!                   temperature     pfpt     deg celsius (ipts-68)
1004      !!                   salinity        pfps     (ipss-78)
1005      !!                   adiabatic       fsatg    deg. c/decibar
1006      !!----------------------------------------------------------------------
1007      REAL(wp) :: pfps, pfpt, pfphp 
1008      REAL(wp) :: fsatg
1009      !
1010      fsatg = (((-2.1687e-16*pfpt+1.8676e-14)*pfpt-4.6206e-13)*pfphp         &
1011        &    +((2.7759e-12*pfpt-1.1351e-10)*(pfps-35.)+((-5.4481e-14*pfpt    &
1012        &    +8.733e-12)*pfpt-6.7795e-10)*pfpt+1.8741e-8))*pfphp             &
1013        &    +(-4.2393e-8*pfpt+1.8932e-6)*(pfps-35.)                         &
1014        &    +((6.6228e-10*pfpt-6.836e-8)*pfpt+8.5258e-6)*pfpt+3.5803e-5
1015      !
1016    END FUNCTION fsatg
1017    !!======================================================================
1018END MODULE sbcisf
Note: See TracBrowser for help on using the repository browser.