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

Context Navigation

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 7494

Last change on this file since 7494 was 7494, checked in by timgraham, 7 years ago
Addition of extra diagnostic outputs in CMIP6 data request. NB. Will require update to field_def file from shaconemo
Property svn:keywords set to `Id`
File size: 44.2 KB

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

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

Download in other formats: