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.

ocesbc.F90 in branches/dev_002_LIM/NEMO/OPA_SRC/SBC – NEMO

Context Navigation

source: branches/dev_002_LIM/NEMO/OPA_SRC/SBC/ocesbc.F90 @ 827

Last change on this file since 827 was 827, checked in by ctlod, 16 years ago
dev_002_LIM : adapt ocean modules to both LIM 2.0 & 3.0 using key_lim2 & key_lim3
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 36.3 KB

Rev	Line
[3]	1	MODULE ocesbc
	2	!!======================================================================
	3	!! * MODULE ocesbc *
	4	!! Ocean surface boundary conditions
	5	!!======================================================================
	6
	7	!!----------------------------------------------------------------------
[84]	8	!! oce_sbc : ???
	9	!! oce_sbc_dmp : ???
[3]	10	!!----------------------------------------------------------------------
	11	!! * Modules used
	12	USE oce ! dynamics and tracers variables
	13	USE dom_oce ! ocean space domain variables
	14	USE cpl_oce ! coupled ocean-atmosphere variables
[258]	15	USE ice_oce ! sea-ice variable
	16	USE blk_oce ! bulk variables
	17	USE flx_oce ! sea-ice/ocean forcings variables
[3]	18	USE phycst ! Define parameters for the routines
[258]	19	USE taumod ! surface stress forcing
	20	USE flxmod ! thermohaline fluxes
	21	USE flxrnf ! runoffs forcing
[3]	22	USE tradmp ! damping salinity trend
[258]	23	USE dtatem ! ocean temperature data
	24	USE dtasal ! ocean salinity data
	25	USE ocfzpt ! surface ocean freezing point
	26	USE lbclnk ! ocean lateral boundary condition
	27	USE lib_mpp ! distribued memory computing library
[3]	28	USE in_out_manager ! I/O manager
[258]	29	USE prtctl ! Print control
[3]	30
	31	IMPLICIT NONE
	32	PRIVATE
	33
	34	!! * Accessibility
	35	PUBLIC oce_sbc ! routine called by step
	36
	37	!! * Shared module variables
[84]	38	REAL(wp), PUBLIC :: & !:
[19]	39	aplus, aminus, & !:
	40	empold = 0.e0 !: current year freshwater budget correction
[84]	41	REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !:
[19]	42	qt , & !: total surface heat flux (w/m2)
	43	qsr , & !: solar radiation (w/m2)
[632]	44	qla , & !: Latent heat flux (W/m2)
	45	qlw , & !: Long Wave Heat Flux (W/m2)
	46	qsb , & !: Sensible Heat Flux (W/m2)
[19]	47	emp , & !: evaporation minus precipitation (kg/m2/s = mm/s)
	48	emps, & !: evaporation - precipitation (free surface)
	49	qrp , & !: heat flux damping (w/m2)
	50	erp !: evaporation damping (kg/m2/s = mm/s)
[359]	51	#if ! defined key_dynspg_rl
[19]	52	REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !:
	53	dmp !: internal dampind term
[3]	54	#endif
	55
	56	# include "domzgr_substitute.h90"
	57	# include "vectopt_loop_substitute.h90"
[247]	58
[3]	59	!!----------------------------------------------------------------------
[247]	60	!! OPA 9.0 , LOCEAN-IPSL (2005)
	61	!! $Header$
	62	!! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
[3]	63	!!----------------------------------------------------------------------
	64	CONTAINS
	65
[826]	66	#if defined key_lim3 \|\| defined key_lim2
[3]	67	!!----------------------------------------------------------------------
[827]	68	!! 'key_lim3' : LIM2 sea-ice model
	69	!! 'key_lim2' : LIM3 sea-ice model
[3]	70	!!----------------------------------------------------------------------
	71	# if defined key_coupled
	72	!!----------------------------------------------------------------------
	73	!! 'key_coupled' : Coupled Ocean/Atmosphere
	74	!!----------------------------------------------------------------------
	75
	76	SUBROUTINE oce_sbc( kt )
	77	!!---------------------------------------------------------------------
	78	!! * ROUTINE oce_sbc *
	79	!!
	80	!! ** Purpose : Ocean surface boundaries conditions with
	81	!! Louvain la Neuve Sea Ice Model in coupled mode
	82	!!
	83	!! History :
	84	!! 1.0 ! 00-10 (O. Marti) Original code
	85	!! 2.0 ! 02-12 (G. Madec) F90: Free form and module
[359]	86	!! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization
[3]	87	!!----------------------------------------------------------------------
[153]	88	!! * Arguments
	89	INTEGER, INTENT( in ) :: kt ! ocean time step
	90
[3]	91	!! * Local declarations
	92	INTEGER :: ji, jj ! dummy loop indices
[440]	93	REAL(wp) :: ztx, ztaux, zty, ztauy
[3]	94	REAL(wp) :: ztdta, ztgel, zqrp
	95	!!----------------------------------------------------------------------
[632]	96
[3]	97	! 1. initialization to zero at kt = nit000
	98	! ---------------------------------------
[632]	99
[3]	100	IF( kt == nit000 ) THEN
	101	qsr (:,:) = 0.e0
	102	freeze(:,:) = 0.e0
	103	qt (:,:) = 0.e0
	104	qrp (:,:) = 0.e0
	105	emp (:,:) = 0.e0
	106	emps (:,:) = 0.e0
	107	erp (:,:) = 0.e0
[359]	108	#if ! defined key_dynspg_rl
[3]	109	dmp (:,:) = 0.e0
	110	#endif
	111	ENDIF
	112
	113	IF( MOD( kt-1, nfice ) == 0 ) THEN
	114
	115	CALL oce_sbc_dmp ! Computation of internal and evaporation damping terms
	116
	117	! Surface heat flux (W/m2)
	118	! -----------------------
	119
	120	! restoring heat flux
	121	DO jj = 1, jpj
	122	DO ji = 1, jpi
	123	ztgel = fzptn(ji,jj)
	124	#if defined key_dtasst
	125	ztdta = MAX( sst(ji,jj), ztgel )
	126	#else
	127	ztdta = MAX( t_dta(ji,jj,1), ztgel )
	128	#endif
[440]	129	zqrp = dqdt0 * ( tb(ji,jj,1) - ztdta )
[3]	130
	131	qrp(ji,jj) = (1.0-freeze(ji,jj) ) * zqrp
	132	END DO
	133	END DO
	134
	135	! non solar heat flux + solar flux + restoring
[296]	136	qt (:,:) = fnsolar(:,:) + fsolar(:,:) + qrp(:,:)
[3]	137
	138	! solar flux
	139	qsr(:,:) = fsolar(:,:)
	140
[359]	141	#if ! defined key_dynspg_rl
[3]	142	! total concentration/dilution effect (use on SSS)
	143	emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)
	144
	145	! total volume flux (use on sea-surface height)
	146	emp (:,:) = fmass(:,:) - dmp(:,:) + runoff(:,:) + erp(:,:)
	147	#else
	148	! Rigid-lid (emp=emps=E-P-R+Erp)
	149	! freshwater flux
	150	emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:)
	151	emp (:,:) = emps(:,:)
	152	#endif
	153
	154	DO jj = 1, jpjm1
	155	DO ji = 1, fs_jpim1 ! vertor opt.
	156	ztx = 0.5 * ( freeze(ji+1,jj) + freeze(ji+1,jj+1) )
[827]	157	zty = 0.5 * ( freeze(ji,jj+1) + freeze(ji+1,jj+1) )
	158	#if defined key_lim3
	159	ztaux = ftaux(ji,jj)
	160	ztauy = ftauy(ji,jj)
	161	#elif defined key_lim2
[3]	162	ztaux = 0.5 * ( ftaux (ji+1,jj) + ftaux (ji+1,jj+1) )
[827]	163	ztauy = 0.5 * ( ftauy (ji,jj+1) + ftauy (ji+1,jj+1) )
	164	#endif
[3]	165	taux(ji,jj) = (1.0-ztx) * taux(ji,jj) + ztx * ztaux
	166	tauy(ji,jj) = (1.0-zty) * tauy(ji,jj) + zty * ztauy
	167	END DO
	168	END DO
	169	CALL lbc_lnk( taux, 'U', -1. )
	170	CALL lbc_lnk( tauy, 'V', -1. )
	171
	172	! Re-initialization of fluxes
[84]	173	sst_io(:,:) = 0.e0
	174	sss_io(:,:) = 0.e0
	175	u_io (:,:) = 0.e0
	176	v_io (:,:) = 0.e0
	177	gtaux (:,:) = 0.e0
	178	gtauy (:,:) = 0.e0
[3]	179
	180	ENDIF
	181
	182	END SUBROUTINE oce_sbc
	183
[473]	184	# elif defined key_flx_bulk_monthly \|\| defined key_flx_bulk_daily \|\| defined key_flx_core
[3]	185	!!----------------------------------------------------------------------
[826]	186	!! 'key_lim3' with LIM sea-ice model
[3]	187	!!----------------------------------------------------------------------
	188
	189	SUBROUTINE oce_sbc( kt )
	190	!!---------------------------------------------------------------------
	191	!! * ROUTINE oce_sbc *
	192	!!
	193	!! ** Purpose : - Ocean surface boundary conditions with LIM sea-ice
	194	!! model in forced mode using bulk formulea
	195	!!
	196	!! History :
	197	!! 1.0 ! 99-11 (M. Imbard) Original code
	198	!! ! 01-03 (D. Ludicone, E. Durand, G. Madec) free surf.
	199	!! 2.0 ! 02-09 (G. Madec, C. Ethe) F90: Free form and module
[359]	200	!! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization
[3]	201	!!----------------------------------------------------------------------
	202	!! * arguments
	203	INTEGER, INTENT( in ) :: kt ! ocean time step
	204
	205	!! * Local declarations
[19]	206	INTEGER :: ji, jj ! dummy loop indices
	207	REAL(wp) :: ztx, ztaux, zty, ztauy
[3]	208	!!----------------------------------------------------------------------
	209
	210	! 1. initialization to zero at kt = nit000
	211	! ---------------------------------------
[632]	212
[3]	213	IF( kt == nit000 ) THEN
	214	qsr (:,:) = 0.e0
	215	qt (:,:) = 0.e0
[632]	216	qla (:,:) = 0.e0
	217	qlw (:,:) = 0.e0
	218	qsb (:,:) = 0.e0
[3]	219	qrp (:,:) = 0.e0
	220	emp (:,:) = 0.e0
	221	emps (:,:) = 0.e0
	222	erp (:,:) = 0.e0
[359]	223	#if ! defined key_dynspg_rl
[3]	224	dmp (:,:) = 0.e0
	225	#endif
	226	ENDIF
[632]	227	#if defined key_flx_core
	228	qla(:,:) = (1 - freeze(:,:))qla_oce + freeze(:,:)qla_ice(:,:)
	229	qsb(:,:) = (1 - freeze(:,:))qsb_oce + freeze(:,:)qsb_ice(:,:)
	230	qlw(:,:) = (1 - freeze(:,:))qlw_oce + freeze(:,:)qlw_ice(:,:)
	231	#endif
[3]	232	IF( MOD( kt-1, nfice ) == 0 ) THEN
	233
	234	CALL oce_sbc_dmp ! Computation of internal and evaporation damping terms
	235
	236	! Surface Ocean fluxes
	237	! ====================
	238
	239	! Surface heat flux (W/m2)
	240	! -----------------
	241
[296]	242	qt (:,:) = fnsolar(:,:) + fsolar(:,:) ! non solar heat flux + solar flux
[3]	243	qsr (:,:) = fsolar(:,:) ! solar flux
	244
[359]	245	#if ! defined key_dynspg_rl
[3]	246	! total concentration/dilution effect (use on SSS)
	247	emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:) + empold
	248
	249	! total volume flux (use on sea-surface height)
	250	emp (:,:) = fmass(:,:) - dmp(:,:) + runoff(:,:) + erp(:,:) + empold
	251	#else
	252	! Rigid-lid (emp=emps=E-P-R+Erp)
	253	emps(:,:) = fmass(:,:) + fsalt(:,:) + runoff(:,:) + erp(:,:) ! freshwater flux
	254	emp (:,:) = emps(:,:)
	255
	256	#endif
	257
	258	! Surface stress
	259	! --------------
	260
	261	! update the stress beloww sea-ice area
	262	DO jj = 1, jpjm1
	263	DO ji = 1, fs_jpim1 ! vertor opt.
	264	ztx = MAX( freezn(ji,jj), freezn(ji,jj+1) ) ! ice/ocean indicator at U- and V-points
	265	zty = MAX( freezn(ji,jj), freezn(ji+1,jj) )
[827]	266	#if defined key_lim3
	267	ztaux = ftaux(ji,jj)
	268	ztauy = ftauy(ji,jj)
	269	#elif defined key_lim2
	270	ztaux = 0.5 * ( ftaux (ji+1,jj) + ftaux (ji+1,jj+1) )
	271	ztauy = 0.5 * ( ftauy (ji,jj+1) + ftauy (ji+1,jj+1) )
	272	#endif
[3]	273	taux(ji,jj) = (1.-ztx) * taux(ji,jj) + ztx * ztaux ! stress at the ocean surface
	274	tauy(ji,jj) = (1.-zty) * tauy(ji,jj) + zty * ztauy
	275	END DO
	276	END DO
	277
	278	! boundary condition on the stress (taux,tauy)
	279	CALL lbc_lnk( taux, 'U', -1. )
	280	CALL lbc_lnk( tauy, 'V', -1. )
	281
	282	! Re-initialization of fluxes
[84]	283	sst_io(:,:) = 0.e0
	284	sss_io(:,:) = 0.e0
	285	u_io (:,:) = 0.e0
	286	v_io (:,:) = 0.e0
[3]	287
	288	ENDIF
	289
	290	END SUBROUTINE oce_sbc
	291
	292	# else
	293	!!----------------------------------------------------------------------
	294	!! Error option LIM sea-ice model requires bulk formulea
	295	!!----------------------------------------------------------------------
	296	This line forced a compilation error
	297	# endif
	298
	299	#else
	300	!!----------------------------------------------------------------------
	301	!! Default option NO LIM sea-ice model
	302	!!----------------------------------------------------------------------
	303	# if defined key_coupled
	304	!!----------------------------------------------------------------------
	305	!! 'key_coupled' : Coupled Ocean/Atmosphere
	306	!!----------------------------------------------------------------------
	307
	308	SUBROUTINE oce_sbc( kt )
	309	!!---------------------------------------------------------------------
	310	!! * ROUTINE oce_sbc *
	311	!!
	312	!! ** Purpose : Ocean surface boundaries conditions in
	313	!! coupled ocean/atmosphere case without sea-ice
	314	!!
	315	!! History :
	316	!! 1.0 ! 00-10 (O. Marti) Original code
	317	!! 2.0 ! 02-12 (G. Madec) F90: Free form and module
	318	!!----------------------------------------------------------------------
	319	!! * Modules used
[258]	320	USE cpl_oce ! coupled ocean-atmosphere variables
[3]	321
	322	!! * Arguments
	323	INTEGER, INTENT( in ) :: kt ! ocean time step index
	324
	325	!! * Local declarations
	326	INTEGER :: ji, jj, jf ! dummy loop indices
[440]	327	REAL(wp) :: ztgel, & ! temporary scalars
[3]	328	zice, zhemis, zqrp, zqri, & ! " "
	329	zq, zqi, zerp, ze, zei, zro ! " "
	330	!!----------------------------------------------------------------------
	331
	332	! Compute fluxes
	333	! --------------
	334
	335	DO jj = 1, jpj
	336	DO ji = 1, jpi
	337
	338	ztgel = fzptn(ji,jj) ! local freezing temperature
	339
	340	! opa model ice freeze()
	341
	342	zice = tmask(ji,jj,1)
	343	IF( tn(ji,jj,1) >= ztgel ) zice = 0.
	344	freeze(ji,jj) = zice
	345
	346	! hemisphere indicator (=1 north, =-1 south)
[632]	347
[3]	348	zhemis = float(isign(1, mjg(jj)-(jpjglo/2+1)))
[632]	349
[3]	350	! a) net downward radiative flux qsr()
	351	! - AGCM qsrc if no ice
	352	! - zero under ice (zice=1)
	353
	354	qsr(ji,jj) = (1.-zice)qsrc(ji,jj)tmask(ji,jj,1)
	355
	356	! b) heat flux damping term qrp()
	357	! - no damping if no ice (zice=0)
	358	! - gamma*min(0,t-tgel) if ice (zice=1)
	359
	360	zqrp = 0.
[440]	361	zqri = dqdt0*MIN( 0., tb(ji,jj,1)-ztgel )
[3]	362	qrp(ji,jj) = ( ( 1. - zice ) * zqrp + zice * zqri ) * tmask(ji,jj,1)
	363
	364
	365	! c) net downward heat flux q() = q0 + qrp()
	366	! for q0
	367	! - AGCM qc if no ice (zice=0)
	368	! - -2 watt/m2 (arctic) or -4 watt/m2 (antarctic) if ice (zice=1)
	369	zq = qc(ji,jj)
	370	zqi = -3. + zhemis
	371	qt(ji,jj) = ( (1.-zice) * zq + zice * zqi ) * tmask(ji,jj,1) + qrp(ji,jj)
[632]	372
[3]	373	! d) water flux damping term erp()
	374	! - no damping
	375	zerp = 0.
	376	erp(ji,jj) = zerp
[632]	377
[3]	378	! e) net upward water flux e() = eo + runoff() + erp()
	379	! for e0
	380	! - AGCM if no ice (zice=0)
	381	! - 1.mm/day if climatological and opa ice (zice=1)
	382	ze = ec(ji,jj)
	383	zei = 1./rday
	384	zro = runoff(ji,jj)
	385	emp(ji,jj) = ( ( 1. - zice ) * ze + zice * zei + zro ) * tmask(ji,jj,1) + erp(ji,jj)
[632]	386
[160]	387	! f) net upward water flux for the salinity surface
	388	! boundary condition
	389	emps(:,:) = emp(:,:)
	390
[3]	391	END DO
	392	END DO
	393
	394	END SUBROUTINE oce_sbc
[473]	395	# elif defined key_flx_bulk_monthly \|\| defined key_flx_bulk_daily \|\| defined key_flx_forced_daily \|\| defined key_flx_core
[3]	396	!!-------------------------------------------------------------------------
[632]	397	!! 'key_flx_bulk_monthly' or 'key_flx_bulk_daily' or core bulk formulea
[3]	398	!! 'key_flx_forced_daily' or no bulk case
	399	!!-------------------------------------------------------------------------
	400
	401	SUBROUTINE oce_sbc( kt )
	402	!!---------------------------------------------------------------------
	403	!! * ROUTINE oce_sbc *
	404	!!
	405	!! ** Purpose : Ocean surface boundary conditions in forced mode
	406	!! using either flux or bulk formulation.
	407	!!
	408	!! History :
	409	!! 1.0 ! 99-11 (M. Imbard) Original code
	410	!! ! 01-03 (D. Ludicone, E. Durand, G. Madec) free surf.
	411	!! 2.0 ! 02-09 (G. Madec, C. Ethe) F90: Free form and module
[359]	412	!! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization
[3]	413	!!----------------------------------------------------------------------
	414	!! * Modules used
[258]	415	USE daymod ! calendar
[3]	416	#if ! defined key_dtasst
[258]	417	USE dtasst, ONLY : rclice ! sea surface temperature data
[3]	418	#endif
[632]	419	#if defined key_flx_bulk_monthly \|\| defined key_flx_bulk_daily \|\| defined key_flx_core
[258]	420	USE blk_oce ! bulk variables
[3]	421	#endif
	422	#if defined key_flx_forced_daily
[258]	423	USE flx_oce ! sea-ice/ocean forcings variables
[3]	424	#endif
	425
	426	!! * arguments
	427	INTEGER, INTENT( in ) :: kt ! ocean time step
	428
	429	!! * local declarations
	430	INTEGER :: ji, jj ! dummy loop arguments
	431	INTEGER :: i15, ifreq !
	432	REAL(wp) :: zxy
	433	REAL(wp) :: zsice, zqri, zqrp, ztdta, zqrj
[440]	434	REAL(wp) :: zq, zqi, zhemis
[3]	435	REAL(wp), DIMENSION(jpi,jpj) :: zeri, zerps, ziclim
	436	REAL(wp), DIMENSION(jpi,jpj) :: zqt, zqsr, zemp
	437	!!----------------------------------------------------------------------
	438
[473]	439
	440	#if defined key_flx_core
	441	CALL ctl_stop( 'flxcore and no ice model not tested yet' )
	442	#endif
	443
[3]	444	! 1. initialization to zero at kt = nit000
	445	! ---------------------------------------
[632]	446
[3]	447	IF( kt == nit000 ) THEN
	448	qsr (:,:) = 0.e0
	449	freeze (:,:) = 0.e0
	450	qt (:,:) = 0.e0
[632]	451	qla (:,:) = 0.e0
	452	qlw (:,:) = 0.e0
	453	qsb (:,:) = 0.e0
[3]	454	qrp (:,:) = 0.e0
	455	emp (:,:) = 0.e0
	456	emps (:,:) = 0.e0
	457	erp (:,:) = 0.e0
[359]	458	#if ! defined key_dynspg_rl
[3]	459	dmp (:,:) = 0.e0
	460	#endif
	461	ENDIF
	462
[632]	463	#if defined key_flx_bulk_monthly \|\| defined key_flx_bulk_daily \|\| defined key_flx_core
[3]	464	ifreq = nfbulk
	465	zqt (:,:) = qsr_oce(:,:) + qnsr_oce(:,:)
	466	zqsr(:,:) = qsr_oce(:,:)
	467	zemp(:,:) = evap(:,:) - tprecip(:,:)
	468	#endif
	469
	470	#if defined key_flx_forced_daily
	471	ifreq = 1
	472	zqt (:,:) = p_qt (:,:)
	473	zqsr(:,:) = p_qsr(:,:)
	474	zemp(:,:) = p_emp(:,:)
	475	#endif
[632]	476	#if defined key_flx_core
	477	qla(:,:) = (1 - freeze(:,:))qla_oce + freeze(:,:)qla_ice(:,:)
	478	qsb(:,:) = (1 - freeze(:,:))qsb_oce + freeze(:,:)qsb_ice(:,:)
	479	qlw(:,:) = (1 - freeze(:,:))qlw_oce + freeze(:,:)qlw_ice(:,:)
	480	#endif
[3]	481	IF( MOD( kt-1, ifreq) == 0 ) THEN
	482	! Computation of internal and evaporation damping terms
	483	CALL oce_sbc_dmp
	484
	485	zsice = - 0.04 / 0.8 ! ratio of isohaline compressibility over isotherme compressibility
	486	! ( d rho / dt ) / ( d rho / ds ) ( s = 34, t = -1.8 )
	487	! Flux computation
	488	DO jj = 1, jpj
	489	DO ji = 1, jpi
	490	! climatological ice
	491	ziclim(ji,jj) = FLOAT( NINT( rclice(ji,jj,1) ) )
	492
	493	! avoid surfreezing point
	494	tn(ji,jj,1) = MAX( tn(ji,jj,1), fzptn(ji,jj) )
	495
	496	! hemisphere indicator (=1 north, =-1 south)
	497	zhemis = FLOAT( isign(1, mjg(jj) - (jpjdta/2+1) ) )
	498
	499	! restoring temperature (ztdta >= to local freezing temperature)
	500	#if defined key_dtasst
	501	ztdta = MAX( sst(ji,jj), fzptn(ji,jj) )
	502	#else
	503	ztdta = MAX( t_dta(ji,jj,1), fzptn(ji,jj) )
	504	#endif
	505
	506	! a) net downward radiative flux qsr()
[227]	507	qsr(ji,jj) = (1.-ziclim(ji,jj)) * zqsr(ji,jj) * tmask(ji,jj,1)
[3]	508
	509	! b) heat flux damping term qrp()
	510	! - gamma*(t-tlevitus) if no climatological ice (ziclim=0)
	511	! - gamma*(t-(tgel-1.)) if climatological ice and no opa ice (ziclim=1 zicopa=0)
	512	! - gamma*min(0,t-tgel) if climatological and opa ice (ziclim=1 zicopa=1)
	513
[440]	514	zqri = dqdt0 * ( tb(ji,jj,1) - ( fzptn(ji,jj) - 1.) )
	515	zqrj = dqdt0 * MIN( 0., tb(ji,jj,1) - fzptn(ji,jj) )
[3]	516
[227]	517	qrp(ji,jj) = ( ziclim(ji,jj) * ( (1 - freeze(ji,jj)) * zqri &
	518	& + freeze(ji,jj) * zqrj ) ) * tmask(ji,jj,1)
	519
[632]	520	#if ! defined key_flx_bulk_monthly && ! defined key_flx_bulk_daily && ! defined key_flx_core
[440]	521	zqrp = dqdt0 * ( tb(ji,jj,1) - ztdta )
[227]	522	qrp(ji,jj) = qrp(ji,jj) + (1. - ziclim(ji,jj)) * zqrp
	523	# endif
	524
[3]	525	! c) net downward heat flux q() = q0 + qrp()
	526	! for q0
	527	! - ECMWF fluxes if no climatological ice (ziclim=0)
	528	! - qrp if climatological ice and no opa ice (ziclim=1 zicopa=0)
	529	! - -2 watt/m2 (arctic) or -4 watt/m2 (antarctic) if climatological and opa ice
	530	! (ziclim=1 zicopa=1)
	531	zq = zqt(ji,jj)
	532	zqi = -3. + zhemis
	533	qt (ji,jj) = ( (1.-ziclim(ji,jj)) * zq &
	534	+ziclim(ji,jj) * freeze(ji,jj) * zqi ) &
	535	* tmask(ji,jj,1) &
	536	+ qrp(ji,jj)
	537
	538	END DO
	539	END DO
	540
[359]	541	#if ! defined key_dynspg_rl
[3]	542	! Free-surface
	543
	544	! Water flux for zero buoyancy flux if no opa ice and ice clim
	545	zeri(:,:) = -zsice * qrp(:,:) * ro0cpr * rauw / 34.0
	546	zerps(:,:) = ziclim(:,:) * ( (1-freeze(:,:)) * zeri(:,:) )
	547
	548	! Contribution to sea level:
	549	! net upward water flux emp() = e-p + runoff() + erp() + dmp + empold
	550	emp (:,:) = zemp(:,:) & ! e-p data
	551	& + runoff(:,:) & ! runoff data
	552	& + erp(:,:) & ! restoring term to SSS data
	553	& + dmp(:,:) & ! freshwater flux associated with internal damping
	554	& + empold ! domain averaged annual mean correction
	555
	556	! Contribution to salinity:
	557	! net upward water flux emps() = e-p + runoff() + erp() + zerps + empold
	558	emps(:,:) = zemp(:,:) &
	559	& + runoff(:,:) &
	560	& + erp(:,:) &
	561	& + zerps(:,:) &
	562	& + empold
	563	#else
	564	! Rigid-lid (emp=emps=E-P-R+Erp)
	565	! freshwater flux
	566	zeri(:,:) = -zsice * qrp(:,:) * ro0cpr * rauw / 34.0
	567	zerps(:,:) = ziclim(:,:) * ( (1-freeze(:,:)) * zeri(:,:) )
	568	emps (:,:) = zemp(:,:) &
	569	& + runoff(:,:) &
	570	& + erp(:,:) &
	571	& + zerps(:,:)
	572	emp (:,:) = emps(:,:)
	573	#endif
	574
	575	! Boundary condition on emp for free surface option
	576	! -------------------------------------------------
	577	CALL lbc_lnk( emp, 'T', 1. )
[632]	578
[3]	579	ENDIF
	580
	581	END SUBROUTINE oce_sbc
	582
	583	# else
	584	!!----------------------------------------------------------------------
	585	!! Default option : Analytical forcing
	586	!!----------------------------------------------------------------------
	587
	588	SUBROUTINE oce_sbc( kt )
	589	!!---------------------------------------------------------------------
	590	!! * ROUTINE oce_sbc *
	591	!!
	592	!! ** Purpose : provide the thermohaline fluxes (heat and freshwater)
[93]	593	!! to the ocean at each time step.
[3]	594	!!
	595	!! ** Method : Constant surface fluxes (read in namelist (namflx))
	596	!!
[296]	597	!! ** Action : - qt, qsr, emp, emps, qrp, erp
[3]	598	!!
	599	!! History :
	600	!! ! 91-03 () Original code
	601	!! 8.5 ! 02-09 (G. Madec) F90: Free form and module
[93]	602	!! 9.0 ! 04-05 (A. Koch-Larrouy) Add Gyre configuration
[3]	603	!!----------------------------------------------------------------------
	604	!! * Modules used
	605	USE flxrnf ! ocean runoffs
[93]	606	USE daymod, ONLY : nyear ! calendar
[440]	607	USE dtasss ! sea surface salinity data
[3]	608
	609	!! * arguments
	610	INTEGER, INTENT( in ) :: kt ! ocean time step
	611
	612	!! * local declarations
	613	REAL(wp) :: & !!! surface fluxes namelist (namflx)
	614	q0 = 0.e0, & ! net heat flux
	615	qsr0 = 0.e0, & ! solar heat flux
	616	emp0 = 0.e0 ! net freshwater flux
[440]	617	REAL(wp) :: zsrp, &
	618	zemp_S, zemp_N, zemp_sais, &
	619	zTstar, zcos_sais1, zconv, &
	620	zcos_sais2
[93]	621	REAL(wp) :: &
	622	zsumemp, & ! tampon used for the emp sum
	623	zsurf, & ! tampon used for the domain sum
	624	ztime, & ! time in hour
[434]	625	ztimemax1, ztimemin1, & ! 21th june, and 21th december if date0 = 1st january
	626	ztimemax2, ztimemin2 ! 21th august, and 21th february if date0 = 1st january
[93]	627	REAL(wp), DIMENSION(jpi,jpj) :: t_star
[434]	628	INTEGER :: ji, jj ! dummy loop indices
	629
[93]	630	INTEGER :: &
	631	zyear0, & ! initial year
	632	zmonth0, & ! initial month
	633	zday0, & ! initial day
[434]	634	zday_year0 ! initial day since january 1st
[3]	635
	636	NAMELIST/namflx/ q0, qsr0, emp0
	637	!!---------------------------------------------------------------------
	638
[93]	639	!same temperature, E-P as in HAZELEGER 2000
[3]	640
[93]	641	IF( cp_cfg == 'gyre' ) THEN
[3]	642
[434]	643	zyear0 = ndate0 / 10000 ! initial year
	644	zmonth0 = ( ndate0 - zyear0 * 10000 ) / 100 ! initial month
	645	zday0 = ndate0 - zyear0 * 10000 - zmonth0 * 100 ! initial day betwen 1 and 30
	646	zday_year0= (zmonth0-1)*30.+zday0 ! initial day betwen 1 and 360
[3]	647
[434]	648	! current day (in hours) since january the 1st of the current year
	649	ztime = FLOAT( kt ) * rdt / (rmmss * rhhmm) & ! total incrementation (in hours)
	650	& - (nyear - 1) * rjjhh * raajj ! minus years since beginning of experiment (in hours)
	651
	652	! 21th june at 24h in hours
	653	ztimemax1 = ((5.30.)+21.) 24.
	654	! 21th december in hours
[167]	655	! rjjhh * raajj / 4 = 1 seasonal cycle in hours
[434]	656	ztimemin1 = ztimemax1 + rjjhh * raajj / 2
	657	! 21th july at 24h in hours
	658	ztimemax2 = ((6.30.)+21.) 24.
	659	! 21th january day in hours
	660	! rjjhh * raajj / 4 = 1 seasonal cycle in hours
	661	ztimemin2 = ztimemax2 - rjjhh * raajj / 2
	662
[93]	663	! amplitudes
	664	zemp_S = 0.7 ! intensity of COS in the South
	665	zemp_N = 0.8 ! intensity of COS in the North
	666	zemp_sais= 0.1
[167]	667	zTstar = 28.3 ! intemsity from 28.3 a -5 deg
[434]	668
[167]	669	! 1/2 period between 21th June and 21th December
[434]	670	zcos_sais1 = COS( (ztime - ztimemax1) / (ztimemin1 - ztimemax1) * rpi )
	671
	672	! 1/2 period between 21th July and 21th January
	673	zcos_sais2 = COS( (ztime - ztimemax2) / (ztimemax2 - ztimemin2) * rpi )
	674
[93]	675	zconv = 3.16e-5 ! convert 1m/yr->3.16e-5mm/s
	676	DO jj = 1, jpj
[167]	677	DO ji = 1, jpi
	678	! domain from 15 deg to 50 deg between 27 and 28 degC at 15N, -3
	679	! and 13 degC at 50N 53.5 + or - 11 = 1/4 period :
	680	! 64.5 in summer, 42.5 in winter
[434]	681	t_star (ji,jj) = zTstar * ( 1 + 1. / 50. * zcos_sais2 ) &
[167]	682	& * COS( rpi * (gphit(ji,jj) - 5.) &
[434]	683	& / (53.5 * ( 1 + 11 / 53.5 * zcos_sais2 ) * 2.) )
[440]	684	qt (ji,jj) = dqdt0 * ( tb(ji,jj,1) - t_star(ji,jj) )
[167]	685	IF( gphit(ji,jj) >= 14.845 .AND. 37.2 >= gphit(ji,jj)) THEN
	686	! zero at 37.8 deg, max at 24.6 deg
	687	emp (ji,jj) = zemp_S * zconv &
	688	& * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (24.6 - 37.2) ) &
[434]	689	& * ( 1 - zemp_sais / zemp_S * zcos_sais1)
[167]	690	emps (ji,jj) = emp (ji,jj)
	691	ELSE
	692	! zero at 37.8 deg, max at 46.8 deg
	693	emp (ji,jj) = - zemp_N * zconv &
	694	& * SIN( rpi / 2 * (gphit(ji,jj) - 37.2) / (46.8 - 37.2) ) &
[434]	695	& * ( 1 - zemp_sais / zemp_N * zcos_sais1 )
	696	emps (ji,jj) = emp (ji,jj)
[167]	697	ENDIF
	698	! 23.5 deg : tropics
[434]	699	qsr (ji,jj) = 230 * COS( 3.1415 * ( gphit(ji,jj) - 23.5 * zcos_sais1 ) / ( 0.9 * 180 ) )
[167]	700	END DO
[93]	701	END DO
[531]	702
	703	! Compute the emp flux such as its integration on the whole domain at each time is zero
	704	IF( nbench /= 1 .AND. nbit_cmp /= 1 ) THEN
	705	zsumemp = 0.e0 ; zsurf = 0.e0
	706	DO jj = 1, jpj
	707	DO ji = 1, jpi
	708	zsumemp = zsumemp + emp(ji,jj) * tmask(ji,jj,1) * tmask_i(ji,jj)
	709	zsurf = zsurf + tmask(ji,jj,1) * tmask_i(ji,jj)
	710	END DO
[434]	711	END DO
[3]	712
[531]	713	IF( lk_mpp ) CALL mpp_sum( zsumemp ) ! sum over the global domain
	714	IF( lk_mpp ) CALL mpp_sum( zsurf ) ! sum over the global domain
[93]	715
	716	! Default GYRE configuration
	717	zsumemp = zsumemp / zsurf
[531]	718	ELSE
	719	! Benchmark GYRE configuration (to allow the bit to bit comparison between Mpp/Mono case)
	720	zsumemp = 0.e0 ; zsurf = 0.e0
[93]	721	ENDIF
[434]	722
[440]	723	IF( lk_dtasss ) THEN ! Sea surface salinity damping
	724	! deds0 from (mm/day) => (mm/s)
	725	zsrp = deds0 / rday
	726	erp(:,:) = zsrp * ( sss(:,:) - sb(:,:,1) ) / ( sn(:,:,1) + 1.e-20 )
	727	ELSE ! NO Sea surface salinity damping
	728	zsrp = 0.e0
	729	erp(:,:) = 0.e0
	730	ENDIF
	731
	732	! freshwater fluxes terms
[93]	733	DO jj = 1, jpj
[434]	734	DO ji = 1, jpi
	735	emp(ji, jj) = emp(ji,jj) - zsumemp * tmask(ji,jj,1)
[440]	736	emps(ji, jj)= emp(ji,jj) + erp(ji,jj)
[434]	737	END DO
[93]	738	END DO
	739
[434]	740	IF( kt == nit000 .AND. lwp ) THEN
	741	WRITE(numout,*)' ocesbc : analytical formulation for gyre'
	742	WRITE(numout,*)' ~~~~~~~ '
	743	WRITE(numout,*)' nyear = ', nyear
	744	WRITE(numout,*)' ztime = ',ztime
	745	WRITE(numout,*)' ztimemax1 = ',ztimemax1
	746	WRITE(numout,*)' ztimemin1 = ',ztimemin1
	747	WRITE(numout,*)' ztimemax2 = ',ztimemax2
	748	WRITE(numout,*)' ztimemin2 = ',ztimemin2
	749	WRITE(numout,*)' zyear0 = ', zyear0
	750	WRITE(numout,*)' zmonth0 = ', zmonth0
	751	WRITE(numout,*)' zday0 = ', zday0
	752	WRITE(numout,*)' zday_year0 = ',zday_year0
	753	WRITE(numout,*)' raajj = ', raajj
	754	WRITE(numout,*)' zemp_S = ',zemp_S
	755	WRITE(numout,*)' zemp_N = ',zemp_N
	756	WRITE(numout,*)' zemp_sais = ',zemp_sais
	757	WRITE(numout,*)' zTstar = ',zTstar
	758	WRITE(numout,*)' zsumemp = ',zsumemp
	759	WRITE(numout,*)' zsurf = ',zsurf
[440]	760	WRITE(numout,*)' dqdt0 = ',dqdt0
	761	WRITE(numout,*)' deds0 = ',deds0
[434]	762	WRITE(numout,*)' zconv = ',zconv
	763	ENDIF
	764
[93]	765	ELSE
	766
[440]	767	! Constant surface fluxes
	768
[93]	769	IF( kt == nit000 ) THEN
[440]	770	IF(lwp) THEN
	771	WRITE(numout,*)' '
	772	WRITE(numout,*)' ocesbc : Constant surface fluxes read in namelist'
	773	WRITE(numout,*)' ~~~~~~~ '
	774	WRITE(numout,*)' Namelist namflx: set the constant flux values'
	775	WRITE(numout,*)' net heat flux q0 = ', q0 , ' W/m2'
	776	WRITE(numout,*)' solar heat flux qsr0 = ', qsr0, ' W/m2'
	777	WRITE(numout,*)' net heat flux emp0 = ', emp0, ' W/m2'
	778	ENDIF
[93]	779
	780	qt (:,:) = q0
	781	qsr (:,:) = qsr0
	782	emp (:,:) = emp0
	783	emps (:,:) = emp0
	784	qrp (:,:) = 0.e0
	785	erp (:,:) = 0.e0
[632]	786
[93]	787	runoff(:,:) = 0.e0
	788	ENDIF
[440]	789
[19]	790	ENDIF
[3]	791
[19]	792	END SUBROUTINE oce_sbc
[3]	793
	794	# endif
	795	#endif
	796
[19]	797	#if defined key_dtasal
	798	!!----------------------------------------------------------------------
	799	!! 'key_dtasal' salinity data
	800	!!----------------------------------------------------------------------
[3]	801	SUBROUTINE oce_sbc_dmp
	802	!!---------------------------------------------------------------------
[19]	803	!! * ROUTINE oce_sbc_dmp *
[3]	804	!!
	805	!! ** Purpose : Computation of internal and evaporation damping terms
	806	!! for ocean surface boundary conditions
	807	!!
	808	!! History :
[19]	809	!! 9.0 ! 04-01 (G. Madec, C. Ethe) Original code
[359]	810	!! 9.0 ! 05-11 (V. Garnier) Surface pressure gradient organization
[3]	811	!!----------------------------------------------------------------------
	812	!! * Local declarations
	813	INTEGER :: ji, jj ! dummy loop indices
	814	REAL(wp), DIMENSION(jpi,jpj) :: zsss, zfreeze
[440]	815	REAL(wp) :: zerp, zsrp
[258]	816	CHARACTER (len=71) :: charout
[359]	817	#if ! defined key_dynspg_rl
[19]	818	REAL(wp) :: zwei
	819	REAL(wp) :: zerpplus(jpi,jpj), zerpminus(jpi,jpj)
	820	REAL(wp) :: zplus, zminus, zadefi
[3]	821	# if defined key_tradmp
	822	INTEGER jk
	823	REAL(wp), DIMENSION(jpi,jpj) :: zstrdmp
	824	# endif
	825	#endif
	826	!!----------------------------------------------------------------------
	827
[826]	828	#if defined key_lim3 \|\| defined key_lim2
[3]	829	! sea ice indicator (1 or 0)
	830	DO jj = 1, jpj
	831	DO ji = 1, jpi
[84]	832	freezn(ji,jj) = MAX(0., SIGN(1., freeze(ji,jj)-rsmall) )
[3]	833	END DO
	834	END DO
	835	zsss (:,:) = sss_io(:,:)
	836	zfreeze(:,:) = freezn(:,:)
	837	#else
[227]	838	zsss (:,:) = sb (:,:,1)
[3]	839	zfreeze(:,:) = freeze(:,:)
	840	#endif
	841
	842	! Initialisation
	843	! --------------
	844	! Restoring coefficients on SST and SSS
[440]	845	zsrp = dqdt0 * ro0cpr * rauw ! (Kg/m2/s)
[3]	846
[359]	847	#if ! defined key_dynspg_rl
[3]	848	! Free-surface
	849
	850	! Internal damping
	851	# if defined key_tradmp
	852	! Vertical mean of dampind trend (computed in tradmp module)
	853	zstrdmp(:,:) = 0.e0
	854	DO jk = 1, jpk
	855	zstrdmp(:,:) = zstrdmp(:,:) + strdmp(:,:,jk) * fse3t(:,:,jk)
	856	END DO
	857	! volume flux associated to internal damping to climatology
[84]	858	dmp(:,:) = zstrdmp(:,:) * rauw / ( zsss(:,:) + 1.e-20 )
[3]	859	# else
	860	dmp(:,:) = 0.e0 ! No internal damping
	861	# endif
	862
	863	! evaporation damping term ( Surface restoring )
	864	zerpplus (:,:) = 0.e0
	865	zerpminus(:,:) = 0.e0
	866	zplus = 15. / rday
	867	zminus = -15. / rday
	868
	869	DO jj = 1, jpj
	870	DO ji = 1, jpi
	871	zerp = ( 1. - 2.upsrnfh(ji,jj) ) zsrp &
	872	& * ( zsss(ji,jj) - s_dta(ji,jj,1) ) &
[84]	873	& / ( zsss(ji,jj) + 1.e-20 )
[3]	874
	875	zerp = MIN( zerp, zplus )
	876	zerp = MAX( zerp, zminus )
	877	erp(ji,jj) = zerp
	878	zerpplus (ji,jj) = MAX( erp(ji,jj), 0.e0 )
	879	zerpminus(ji,jj) = MIN( erp(ji,jj), 0.e0 )
	880	END DO
	881	END DO
	882
	883	aplus = 0.e0
	884	aminus = 0.e0
[632]	885
[531]	886	IF( nbit_cmp == 1) THEN
[632]	887
[531]	888	IF(ln_ctl) THEN
	889	WRITE(charout,FMT="('oce_sbc_dmp : a+ = ',D23.16, ' a- = ',D23.16)") aplus, aminus
	890	CALL prt_ctl_info(charout)
	891	ENDIF
	892	erp(:,:) = 0.e0
	893
	894	ELSE
	895
	896	DO jj = 1, jpj
	897	DO ji = 1, jpi
	898	zwei = e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj)
	899	aplus = aplus + zerpplus (ji,jj) * zwei
	900	aminus = aminus - zerpminus(ji,jj) * zwei
	901	END DO
[3]	902	END DO
[531]	903	IF( lk_mpp ) CALL mpp_sum( aplus ) ! sums over the global domain
	904	IF( lk_mpp ) CALL mpp_sum( aminus )
[3]	905
[531]	906	IF(ln_ctl) THEN
	907	WRITE(charout,FMT="('oce_sbc_dmp : a+ = ',D23.16, ' a- = ',D23.16)") aplus, aminus
	908	CALL prt_ctl_info(charout)
	909	ENDIF
[258]	910
[531]	911	zadefi = MIN( aplus, aminus )
	912	IF( zadefi == 0.e0 ) THEN
	913	erp(:,:) = 0.e0
	914	ELSE
	915	erp(:,:) = zadefi * ( zerpplus(:,:) / aplus + zerpminus(:,:) / aminus )
	916	ENDIF
	917
	918	END IF
[3]	919	#else
	920	! Rigid-lid (emp=emps=E-P-R+Erp)
[632]	921
[3]	922	erp(:,:) = ( 1. - zfreeze(:,:) ) * zsrp & ! surface restoring term
	923	& * ( zsss(:,:) - s_dta(:,:,1) ) &
[84]	924	& / ( zsss(:,:) + 1.e-20 )
[3]	925	#endif
	926
	927	END SUBROUTINE oce_sbc_dmp
	928
[19]	929	#else
	930	!!----------------------------------------------------------------------
	931	!! Dummy routine NO salinity data
	932	!!----------------------------------------------------------------------
	933	SUBROUTINE oce_sbc_dmp ! Dummy routine
	934	WRITE(,) 'oce_sbc_dmp: you should not have seen that print! error?'
	935	END SUBROUTINE oce_sbc_dmp
	936	#endif
	937
[3]	938	!!======================================================================
	939	END MODULE ocesbc

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

Download in other formats: