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

Context Navigation

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 9366

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

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

Download in other formats: