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

Context Navigation

source: branches/UKMO/dev_r5518_coupling_GSI7_GSI8_landice_bitcomp/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 6675

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

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

Download in other formats: