Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sbcisf.F90 @ 11200

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 11200

Download in other formats: