Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

PHY_Atm_AT_RUN.f90 @ 2222

Last change on this file since 2222 was 2089, checked in by Laurent Fairhead, 10 years ago

Inclusion de la physique de MAR

Integration of MAR physics

File size: 13.8 KB

Rev	Line
[2089]	1	subroutine PHY_Atm_AT_RUN(FlagSV_KzT,FlagCM)
	2
	3	!------------------------------------------------------------------------------+
	4	! Sat 22-Jun-2013 MAR \|
	5	! MAR PHY_Atm_AT_RUN \|
	6	! subroutine PHY_Atm_AT_RUN performs Turbulent Vertical Diffusion Scheme \|
	7	! \|
	8	! version 3.p.4.1 created by H. Gallee, Wed 13-Mar-2013 \|
	9	! Last Modification by H. Gallee, Sat 22-Jun-2013 \|
	10	! \|
	11	!------------------------------------------------------------------------------+
	12
	13	use Mod_Real
	14	use Mod_PHY____dat
	15	use Mod_PHY____grd
	16	use Mod_PHY_AT_ctr
	17	use Mod_PHY_AT_grd
	18	use Mod_PHY_AT_kkl
	19	use Mod_SISVAT_gpt
	20
	21
	22	IMPLICIT NONE
	23
	24
	25	logical :: FlagSV_KzT
	26	logical :: FlagCM
	27
	28	integer :: i, j, k, ikl
	29
	30
	31
	32
	33	! Assignation of Mod_PHY_AT_kkl
	34	! ================================
	35
	36	DO ikl = 1,kcolp
	37
	38	i = ii__AP(ikl)
	39	j = jj__AP(ikl)
	40
	41	LMO_AT(ikl) = LMO_SV_gpt(ikl)
	42
	43	ENDDO
	44
	45
	46
	47
	48	! Turbulent Kinetic Energy
	49	! ========================
	50
	51	IF (AT_TKE) THEN
	52
	53	! **************
	54	CALL PHY_vdfTKE_RUN
	55	CALL PHY_genTKE_RUN
	56	! **************
	57
	58	ELSE
	59	DO ikl = 1,kcolp
	60	DO k = 1,mzp
	61	TKE_AT(ikl,k)= eps6
	62	eps_AT(ikl,k)= epsn
	63	ENDDO
	64	ENDDO
	65	END IF
	66
	67
	68
	69
	70	! Turbulent Vertical Diffusion
	71	! ============================
	72
	73	NewAAT = .TRUE.
	74	SchmAT = 'u'
	75	! ******
	76	CALL Atm_AT
	77	! ******
	78
	79
	80	NewAAT = .FALSE.
	81	SchmAT = 'v'
	82	! ******
	83	CALL Atm_AT
	84	! ******
	85
	86
	87	IF (.NOT.FlagSV_KzT) THEN
	88
	89	SchmAT = 'T'
	90	! ******
	91	CALL Atm_AT
	92	! ******
	93	ELSE
	94
	95	DO k =1,mzp
	96	DO ikl=1,kcolp
	97	dpktAT(ikl,k) = 0.0000
	98	ENDDO
	99	ENDDO
	100
	101	END IF
	102
	103
	104	SchmAT = 'q'
	105	! ******
	106	CALL Atm_AT
	107	! ******
	108
	109
	110	IF (FlagCM) THEN
	111
	112	SchmAT = 'w'
	113	! ******
	114	CALL Atm_AT
	115	! ******
	116
	117
	118	SchmAT = 'i'
	119	! ******
	120	CALL Atm_AT
	121	! ******
	122
	123
	124	SchmAT = 'C'
	125	! ******
	126	CALL Atm_AT
	127	! ******
	128
	129
	130	SchmAT = 's'
	131	! ******
	132	CALL Atm_AT
	133	! ******
	134
	135
	136	SchmAT = 'r'
	137	! ******
	138	CALL Atm_AT
	139	! ******
	140
	141
	142	END IF
	143
	144
	145	end subroutine PHY_Atm_AT_RUN
	146
	147
	148
	149
	150	subroutine Atm_AT
	151
	152	!------------------------------------------------------------------------------+
	153	! Sat 22-Jun-2013 MAR \|
	154	! MAR Atm_AT \|
	155	! subroutine Atm_AT solves Turbulent Vertical Diffusion Scheme \|
	156	! \|
	157	! version 3.p.4.1 created by H. Gallee, Wed 13-Mar-2013 \|
	158	! Last Modification by H. Gallee, Sat 22-Jun-2013 \|
	159	! \|
	160	!------------------------------------------------------------------------------+
	161
	162	use Mod_Real
	163	use Mod_PHY____dat
	164	use Mod_PHY____grd
	165	use Mod_PHY_AT_grd
	166	use Mod_PHY_AT_kkl
	167	use Mod_PHY_CM_kkl
	168	use Mod_PHY_DY_kkl
	169	use Mod_SISVAT_gpt
	170
	171
	172
	173	IMPLICIT NONE
	174
	175
	176
	177	integer :: i, j, k, ikl
	178	real(kind=real8) :: psa_sq
	179
	180
	181
	182
	183	! ++++++++++++++++
	184	DO ikl = 1,kcolp
	185	! ++++++++++++++++
	186
	187	i = ii__AP(ikl)
	188	j = jj__AP(ikl)
	189
	190
	191
	192	! Tri-Diagonal Matrix: Variable
	193	! =============================
	194
	195	! Wind Speed
	196	! ----------
	197
	198	IF (schmAT.EQ.'a') THEN
	199	DO k=1,mzp
	200	Kz__AT(ikl,k) = 1.0
	201	var_AT(ikl,k) = ua__DY(ikl,k)
	202	ENDDO
	203	k= mzpp
	204	var_AT(ikl,k) = 0.
	205
	206	ELSE IF (schmAT.EQ.'b') THEN
	207	DO k=1,mzp
	208	Kz__AT(ikl,k) = 1.0
	209	var_AT(ikl,k) = va__DY(ikl,k)
	210	ENDDO
	211	k= mzpp
	212	var_AT(ikl,k) = 0.
	213
	214	ELSE IF (schmAT.EQ.'u') THEN
	215	DO k=1,mzp
	216	Kz__AT(ikl,k) = Kzm_AT(ikl,k)
	217	var_AT(ikl,k) = ua__DY(ikl,k)
	218	ENDDO
	219	k= mzpp
	220	var_AT(ikl,k) = 0.
	221
	222	ELSE IF (schmAT.EQ.'v') THEN
	223	DO k=1,mzp
	224	Kz__AT(ikl,k) = Kzm_AT(ikl,k)
	225	var_AT(ikl,k) = va__DY(ikl,k)
	226	ENDDO
	227	k= mzpp
	228	var_AT(ikl,k) = 0.
	229
	230
	231	! Temperature
	232	! -----------
	233
	234	ELSE IF (schmAT.EQ.'T') THEN
	235	DO k=1,mzp
	236	Kz__AT(ikl,k) = Kzh_AT(ikl,k)
	237	var_AT(ikl,k) = pkt_DY(ikl,k)
	238	ENDDO
	239	k= mzpp
	240	var_AT(ikl,k) = pkt_DY(ikl,k)
	241
	242
	243	! Water Species
	244	! -------------
	245
	246	ELSE IF (schmAT.EQ.'q') THEN
	247	DO k=1,mzp
	248	Kz__AT(ikl,k) = Kzh_AT(ikl,k)
	249	var_AT(ikl,k) = qv__DY(ikl,k)
	250	ENDDO
	251	k= mzpp
	252	var_AT(ikl,k) = qv__DY(ikl,mzp) &
	253	& -(Z___DY(ikl,mzp)-Z___DY(ikl,mzpp)) &
	254	& * uqs_SV_gpt(ikl) /Kzh_AT(ikl,mzp )
	255	var_AT(ikl,k) = max(epsq,var_AT(ikl,k ))
	256	var_AT(ikl,k) = min( var_AT(ikl,k ),qvsiCM(ikl,k))
	257	qv__DY(ikl,mzpp) = var_AT(ikl,k )
	258
	259	ELSE IF (schmAT.EQ.'w') THEN
	260	DO k=1,mzp
	261	Kz__AT(ikl,k) = Kzh_AT(ikl,k)
	262	var_AT(ikl,k) = qw__CM(ikl,k)
	263	ENDDO
	264	k= mzpp
	265	var_AT(ikl,k) = qw__CM(ikl,mzp)
	266
	267	ELSE IF (schmAT.EQ.'i') THEN
	268	DO k=1,mzp
	269	Kz__AT(ikl,k) = Kzh_AT(ikl,k)
	270	var_AT(ikl,k) = qi__CM(ikl,k)
	271	ENDDO
	272	k= mzpp
	273	var_AT(ikl,k) = qi__CM(ikl,mzp)
	274
	275	ELSE IF (schmAT.EQ.'C') THEN
	276	DO k=1,mzp
	277	Kz__AT(ikl,k) = Kzh_AT(ikl,k)
	278	var_AT(ikl,k) = CCNiCM(ikl,k)
	279	ENDDO
	280	k= mzpp
	281	var_AT(ikl,k) = CCNiCM(ikl,mzp)
	282
	283	ELSE IF (schmAT.EQ.'s') THEN
	284	DO k=1,mzp
	285	Kz__AT(ikl,k) = Kzh_AT(ikl,k)
	286	var_AT(ikl,k) = qs__CM(ikl,k)
	287	ENDDO
	288	k= mzpp
	289	var_AT(ikl,k) = qs__CM(ikl,k)
	290
	291	ELSE IF (schmAT.EQ.'r') THEN
	292	DO k=1,mzp
	293	Kz__AT(ikl,k) = Kzh_AT(ikl,k)
	294	var_AT(ikl,k) = qr__CM(ikl,k)
	295	ENDDO
	296	k= mzpp
	297	var_AT(ikl,k) = qr__CM(ikl,mzp)
	298
	299	ELSE
	300	write(6,*)
	301	write(6,*) ' schmAT = ',schmAT
	302	write(6,*) ' Inadequate..., NOR for Ekman Spiral, NOR for momentum, NOR for scalars'
	303	write(6,*) ' '
	304
	305	! **************
	306	call MAR________BUG
	307	! **************
	308	STOP
	309	END IF
	310
	311
	312
	313	! Tri-Diagonal Matrix Coefficients
	314	! ================================
	315
	316	IF (NewAAT) THEN
	317	psa_sq = psa_DY(ikl) *psa_DY(ikl)
	318	Ac__AT(ikl,mzp) = Ac0_AT(mzp) &
	319	& * roa_DY(ikl,mzp) *roamDY(ikl,mzp) &
	320	& / psa_sq
	321	DO k=mzp-1,1,-1
	322	Ac__AT(ikl,k ) = Ac0_AT(k ) &
	323	& * roa_DY(ikl,k ) *roamDY(ikl,k ) &
	324	& / psa_sq
	325	Cc__AT(ikl,k+1) = Ac__AT(ikl,k) &
	326	& * roa_DY(ikl,k+1) /roa_DY(ikl,k )
	327	ENDDO
	328	ENDIF
	329
	330	A___AT(ikl,mzp) = Ac__AT(ikl,mzp) * Kz__AT(ikl,mzp)
	331	DO k=mzp-1,1,-1
	332	A___AT(ikl,k ) = Ac__AT(ikl,k ) * Kz__AT(ikl,k )
	333	C___AT(ikl,k+1) = Cc__AT(ikl,k+1) * Kz__AT(ikl,k )
	334	ENDDO
	335	DO k= 2,mzp
	336	B___AT(ikl,k ) = 1.0000 - A___AT(ikl,k ) - C___AT(ikl,k )
	337	ENDDO
	338	B___AT(ikl,1 ) = 1.0000 - A___AT(ikl,1 )
	339	C___AT(ikl,1 ) = 0.0000
	340
	341	D___AT(ikl,1 ) = &
	342	& A___AT(ikl,1 ) * a_b_AT *(var_AT(ikl,1 ) - var_AT(ikl,2 ))
	343
	344	k= mzp
	345	D___AT(ikl,k ) = &
	346	& A___AT(ikl,k ) * (a_b_AT * var_AT(ikl,k ) - var_AT(ikl,k+1) &
	347	& / betaAT) &
	348	& - C___AT(ikl,k ) * a_b_AT *(var_AT(ikl,k-1) - var_AT(ikl,k )) &
	349	& + var_AT(ikl,k )
	350	A___AT(ikl,k ) = 0.0000
	351	DO k= 2,mzp-1
	352	D___AT(ikl,k ) = &
	353	& A___AT(ikl,k ) * a_b_AT *(var_AT(ikl,k ) - var_AT(ikl,k+1)) &
	354	& - C___AT(ikl,k ) * a_b_AT *(var_AT(ikl,k-1) - var_AT(ikl,k )) &
	355	& + var_AT(ikl,k )
	356	ENDDO
	357	D___AT(ikl,1 ) = &
	358	& A___AT(ikl,1 ) * a_b_AT *(var_AT(ikl,1 ) - var_AT(ikl,2 )) &
	359	& + var_AT(ikl,1 )
	360
	361
	362
	363	! Tri-Diagonal Matrix Resolution (Gaussian Elimination / Thomas Algorithm)
	364	! ==============================
	365
	366	! Forward Sweep
	367	! --------------
	368
	369	P___AT(1) = B___AT(ikl,1)
	370	Q___AT(1) = -A___AT(ikl,1) / P___AT(1)
	371	DO k=2,mzp
	372	P___AT(k) = C___AT(ikl,k) * Q___AT(k-1) + B___AT(ikl,k)
	373	Q___AT(k) = -A___AT(ikl,k) / P___AT(k)
	374	ENDDO
	375	X___AT(1) = D___AT(ikl,1) / P___AT(1)
	376	DO k=2,mzp
	377	X___AT(k) = (D___AT(ikl,k) - C___AT(ikl,k) * X___AT(k-1)) / P___AT(k)
	378	END DO
	379
	380
	381	! Backward Sweep
	382	! --------------
	383
	384	DO k= mzp-1,1,-1
	385	X___AT(k) = Q___AT(k) * X___AT(k+1) + X___AT(k)
	386	END DO
	387
	388
	389
	390	! Elimination of too small values
	391	! -------------------------------
	392
	393	IF (SchmAT.EQ.'w'.OR.SchmAT.EQ.'i'.OR.SchmAT.EQ.'C'.OR. &
	394	& SchmAT.EQ.'s'.OR.SchmAT.EQ.'r') THEN
	395	DO k= mzp ,1,-1
	396	IF (X___AT(k).NE.0.0000 .AND. X___AT(k).LT.1.e-18) &
	397	& X___AT(k) = 0.0000
	398	END DO
	399	END IF
	400
	401
	402
	403
	404	! Tendency
	405	! ========
	406
	407	! Wind Speed
	408	! ----------
	409
	410	IF (schmAT.EQ.'a') THEN
	411	DO k=1,mzp
	412	dua_AT(ikl,k) =(X___AT(k) - ua__DY(ikl,k)) / dt__AT
	413	ENDDO
	414	ELSE IF (schmAT.EQ.'b') THEN
	415	DO k=1,mzp
	416	dva_AT(ikl,k) =(X___AT(k) - va__DY(ikl,k)) / dt__AT
	417	ENDDO
	418	ELSE IF (schmAT.EQ.'u') THEN
	419	DO k=1,mzp
	420	dua_AT(ikl,k) =(X___AT(k) - ua__DY(ikl,k)) / dt__AT
	421	ENDDO
	422	ELSE IF (schmAT.EQ.'v') THEN
	423	DO k=1,mzp
	424	dva_AT(ikl,k) =(X___AT(k) - va__DY(ikl,k)) / dt__AT
	425	ENDDO
	426
	427
	428	! Temperature
	429	! -----------
	430
	431	ELSE IF (schmAT.EQ.'T') THEN
	432	DO k=1,mzp
	433	dpktAT(ikl,k) =(X___AT(k) - pkt_DY(ikl,k)) / dt__AT
	434	ENDDO
	435
	436
	437	! Water Species
	438	! -------------
	439
	440	ELSE IF (schmAT.EQ.'q') THEN
	441	DO k=1,mzp
	442	dqv_AT(ikl,k) =(X___AT(k) - qv__DY(ikl,k)) / dt__AT
	443	ENDDO
	444
	445	ELSE IF (schmAT.EQ.'w') THEN
	446	DO k=1,mzp
	447	dqw_AT(ikl,k) =(X___AT(k) - qw__CM(ikl,k)) / dt__AT
	448	ENDDO
	449
	450	ELSE IF (schmAT.EQ.'i') THEN
	451	DO k=1,mzp
	452	dqi_AT(ikl,k) =(X___AT(k) - qi__CM(ikl,k)) / dt__AT
	453	ENDDO
	454
	455	ELSE IF (schmAT.EQ.'C') THEN
	456	DO k=1,mzp
	457	dCi_AT(ikl,k) =(X___AT(k) - CCNiCM(ikl,k)) / dt__AT
	458	ENDDO
	459
	460	ELSE IF (schmAT.EQ.'s') THEN
	461	DO k=1,mzp
	462	dqs_AT(ikl,k) =(X___AT(k) - qs__CM(ikl,k)) / dt__AT
	463	ENDDO
	464
	465	ELSE IF (schmAT.EQ.'r') THEN
	466	DO k=1,mzp
	467	dqr_AT(ikl,k) =(X___AT(k) - qr__CM(ikl,k)) / dt__AT
	468	ENDDO
	469
	470	ELSE
	471	write(6,*)
	472	write(6,*) ' schmAT = ',schmAT
	473	write(6,*) ' Inadequate..., NOR for Ekman Spiral, NOR for momentum, NOR for scalars'
	474	write(6,*) ' '
	475
	476	! **************
	477	call MAR________BUG
	478	! **************
	479	STOP
	480	END IF
	481
	482
	483	! ++++++
	484	END DO
	485	! ++++++
	486
	487
	488
	489
	490	end subroutine Atm_AT

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

Download in other formats:

Original Format