Context Navigation

sbcisf.F90 @ 9383

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

Download in other formats: