Context Navigation

sbcisf.F90 @ 7773

Last change on this file since 7773 was 7773, checked in by mattmartin, 7 years ago

Committing updates after doing the following:

merging the branch dev_r4650_general_vert_coord_obsoper@7763 into this branch
updating it so that the following OBS changes were implemented correctly on top of the simplification changes:
- generalised vertical coordinate for profile obs. This was done so that is now the default option.
- sst bias correction implemented with the new simplified obs code.
- included the biogeochemical obs types int he new simplified obs code.
- included the changes to exclude obs in the boundary for limited area models
- included other changes for the efficiency of the obs operator to remove global arrays.

Property svn:keywords set to Id

File size: 41.8 KB

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

Download in other formats: