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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_isfprint/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 13072

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

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

Download in other formats: