Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sbcisf.F90 @ 11738

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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 11738

Download in other formats: