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

Context Navigation

source: branches/UKMO/dev_r5518_med_test/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 6200

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

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

Download in other formats: