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

Context Navigation

source: branches/UKMO/dev_r5518_debug_isf_restart/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 6369

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

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

Download in other formats: