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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

sbcisf.F90 in branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 7033

Last change on this file since 7033 was 7033, checked in by timgraham, 8 years ago
Added extra diagnostics and corrected some of the surface coupling fields.
File size: 48.9 KB

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

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

Download in other formats: