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 @ 6488

Last change on this file since 6488 was 6488, checked in by davestorkey, 8 years ago

Commit changes from dev_r5518_coupling_GSI7_GSI8_landice and its ancestor branch dev_r5518_CICE_coupling_GSI7_GSI8.
Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r6023 cf. r5668 of /branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice/NEMOGCM@6487

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