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

Context Navigation

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

Last change on this file since 6755 was 6755, checked in by davestorkey, 8 years ago
UKMO/dev_r5518_GO6_package branch: Merge back changes from UKMO/dev_r5518_GO6_package_MEDUSA_temporary Custom merge into /branches/UKMO/dev_r5518_GO6_package/NEMOGCM: r6749 cf. r6618 of /branches/UKMO/dev_r5518_GO6_package_MEDUSA_temporary/NEMOGCM@6754
File size: 46.5 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: