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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

sbcisf.F90 in branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcisf.F90 @ 10775

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

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

Download in other formats: