Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

PHY_Atm_CP_RUN.f90 @ 4794

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

Inclusion de la physique de MAR

Integration of MAR physics

File size: 18.6 KB

Rev	Line
[2089]	1	subroutine PHY_Atm_CP_RUN(mzc,kcolc)
	2
	3	!------------------------------------------------------------------------------+
	4	! Mon 17-Jun-2013 MAR \|
	5	! subroutine PHY_Atm_CP_RUN interfaces \|
	6	! Bechtold et al. (2001) Convective Parameterization with MAR \|
	7	! \|
	8	! version 3.p.4.1 created by H. Gallee, Mon 8-Apr-2013 \|
	9	! Last Modification by H. Gallee, Mon 17-Jun-2013 \|
	10	! \|
	11	!------------------------------------------------------------------------------+
	12	! \|
	13	! INPUT \|
	14	! ^^^^^ it_EXP : Experiment Iteration Counter \|
	15	! dt__CP : Mass Flux Scheme: Time Step \|
	16	! dxHOST : grid spacing of HOST MODEL [m] \|
	17	! Ta__DY(kcolp,mzpp) : air temperature [K] \|
	18	! qs__CM(kcolp,mzpp) : air snow Particl. concentr. [kg/kg] \|
	19	! qr__CM(kcolp,mzp ) : air rain drops concentr. [kg/kg] \|
	20	! \|
	21	! INPUT/OUTPUT \|
	22	! ^^^^^^^^^^^^ qv__DY(kcolp,mzpp) : Specific Humidity [kg/kg] \|
	23	! qw__CM(kcolp,mzp ) : air cloud droplets concentr. [kg/kg] \|
	24	! qi__CM(kcolp,mzp ) : air cloud crystals concentr. [kg/kg] \|
	25	! snowCP(kcolp ) : Snow (convective) [m] \|
	26	! snowCM(kcolp ) : Snow (convective + stratiform) [m] \|
	27	! rainCP(kcolp ) : Rain (convective) [m] \|
	28	! rainCM(kcolp ) : Rain (convective + stratiform) [m] \|
	29	! \|
	30	! OUTPUT dpktCP(kcolp,mzp ) : Reduc. Pot.Temperat.Tendency [K/X/s] \|
	31	! ^^^^^^ dqv_CP(kcolp,mzp ) : Specific Humidity Tendency [kg/kg/s] \|
	32	! dqw_CP(kcolp,mzp ) : cloud dropl.Concent.Tendency [kg/kg/s] \|
	33	! dqi_CP(kcolp,mzp ) : cloud cryst.Concent.Tendency [kg/kg/s] \|
	34	! \|
	35	! dss_CP(kcolp ) : Snow (convective) Tendency [m/s] \|
	36	! drr_CP(kcolp ) : Rain (convective) Tendency [m/s] \|
	37	! \|
	38	! pkt_DY(kcolp,mzp ) : Reduced Potential Temperature [K/X] \|
	39	! \|
	40	! CAPECP(kcolp ) : Convective Avail.Potent.Energy \|
	41	! \|
	42	! REFER. : MesoNH MASS FLUX CONVECTION Routine \|
	43	! ^^^^^^^^ (Bechtold et al., 2001, QJRMS 127, pp 869-886) \|
	44	! \|
	45	! # OPTION : #EW Energy and Water Conservation \|
	46	! # ^^^^^^^^ \|
	47	! \|
	48	! MODIF. HGallee: 18-11-2004: Adaptation to CVAmnh.f90.laurent \|
	49	! ^^^^^^ (Argument kensbl of CONVECTION removed) \|
	50	! \|
	51	!------------------------------------------------------------------------------+
	52
	53	use Mod_Real
	54	use Mod_PHY____dat
	55	use Mod_PHY____grd
	56	use Mod_PHY_CP_ctr
	57	use Mod_PHY_CP_dat
	58	use Mod_PHY_CP_grd
	59	use Mod_PHY_CP_kkl
	60	use Mod_PHY_DY_kkl
	61	use Mod_PHY_AT_kkl
	62	use Mod_PHY_CM_kkl
	63
	64	use Mod_Atm_CP_RUN
	65
	66	IMPLICIT NONE
	67
	68
	69
	70	! Global Variables
	71	! ================
	72
	73	integer :: mzc ! Nb of levels
	74	integer :: kcolc ! Nb of columns
	75
	76
	77
	78
	79	! Local Variables
	80	! ================
	81
	82	real(kind=real8) :: bANA !
	83	real(kind=real8) :: zANA !
	84
	85	integer :: i ! x-Axis Index
	86	integer :: j ! y-Axis Index
	87	integer :: k ! Level Index (from top to bottom)
	88	integer :: klc ! Level Index (from bottom to top )
	89	integer :: ikl ! Column Index
	90
	91	REAL :: Pdxdy0(kcolc)
	92	REAL :: P_pa_0(kcolc,mzc)
	93	REAL :: P_za_0(kcolc,mzc)
	94	REAL :: P_Ta_0(kcolc,mzc)
	95	REAL :: P_Qa_0(kcolc,mzc)
	96	REAL :: P_Qw_0(kcolc,mzc)
	97	REAL :: P_Qi_0(kcolc,mzc)
	98	REAL :: P_Ua_0(kcolc,mzc)
	99	REAL :: P_Va_0(kcolc,mzc)
	100	REAL :: P_Wa_0(kcolc,mzc)
	101
	102	integer :: Kstep1(kcolc) ! convective counter
	103	integer :: K_CbT1(kcolc) ! cloud top level
	104	integer :: K_CbB1(kcolc) ! cloud base level
	105
	106	integer, parameter :: KTCCH0=1 !
	107	REAL :: P_CH_0(kcolc,mzc,KTCCH0)
	108	REAL :: PdCH_1(kcolc,mzc,KTCCH0)
	109
	110	REAL :: PdTa_1(kcolc,mzc)
	111	REAL :: PdQa_1(kcolc,mzc)
	112	REAL :: PdQw_1(kcolc,mzc)
	113	REAL :: PdQi_1(kcolc,mzc)
	114	REAL :: Pdrr_1(kcolc)
	115	REAL :: Pdss_1(kcolc)
	116	REAL :: PuMF_1(kcolc,mzc) ! Upward Mass Flux
	117	REAL :: PdMF_1(kcolc,mzc) ! Downward Mass Flux
	118	REAL :: Pfrr_1(kcolc,mzc) ! Liquid Precipitation Flux
	119	REAL :: Pfss_1(kcolc,mzc) ! Solid Precipitation Flux
	120	REAL :: Pcape1(kcolc) ! CAPE [J/kg]
	121
	122
	123	! Diagnostic Variables
	124	! --------------------
	125
	126	! #EW integer :: irmx ,jrmx ,iter_0 !
	127	! #EW real(kind=real8) :: rr_max,temp_r,energ0 !
	128	! #EW real(kind=real8) :: water0,waterb !
	129
	130
	131
	132	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	133	! !
	134	! ALLOCATION
	135	! ----------
	136
	137	IF (it_RUN.EQ.1 .OR. FlagDALLOC) THEN !
	138	allocate ( wa_ANA(kcolc,mzc) ) ! ANAbatic Wind Speed [m/s]
	139	END IF
	140	! !
	141	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	142
	143
	144
	145
	146	! Update Convective Mass Flux
	147	! ===========================
	148
	149	IF ( mod(iitCV0,jjtCV0).EQ.0 ) THEN
	150
	151	! Martin control
	152	PRINT*,'Dans PHY_Atm_CP_RUN'
	153	PRINT*,'size(ii__AP)',size(ii__AP)
	154	PRINT*,'size(jj__AP)',size(jj__AP)
	155	PRINT*,'size(wa_ANA)=',size(wa_ANA)
	156	PRINT*,'ii__AP(1)=',ii__AP(1)
	157	PRINT*,'jj__AP(1)=',jj__AP(1)
	158	PRINT*,'ii__AP(kcolp)=',ii__AP(kcolp)
	159	PRINT*,'jj__AP(kcolp)=',jj__AP(kcolp)
	160	PRINT*,'mzp=',mzp
	161	! Martin control
	162
	163	DO ikl = 1,kcolp
	164
	165	! PRINT*,'ikl=',ikl
	166
	167	i = ii__AP(ikl)
	168	j = jj__AP(ikl)
	169
	170	! PRINT*,'ii__AP(',ikl,')=',ii__AP(ikl)
	171	! PRINT*,'jj__AP(',ikl,')=',jj__AP(ikl)
	172
	173
	174	! Contribution from Subgrid Mountain Breeze
	175	! -----------------------------------------
	176
	177	DO k=1,mzp
	178	wa_ANA(ikl,k ) = 0.0
	179	END DO
	180	IF(Lo_ANA) THEN
	181	DO k=1,mzp
	182	bANA = min(Z___DY(ikl,k ), zi__AT(ikl))
	183	zANA = hANA(ikl) + bANA * 2.0
	184	IF (Z___DY(ikl,k ) .LE. zANA .AND. &
	185	& Ta__DY(ikl,mzpp) .GT. Ta__DY(ikl,mzp)) THEN
	186	wa_ANA(ikl,k ) = rANA * bANA * 0.5 ! Half Integrated Horizontal Divergence
	187
	188	END IF
	189	END DO
	190	END IF
	191
	192
	193
	194	! Mass Flux convective Scheme: Set Up Vertical Profiles
	195	! -----------------------------------------------------
	196
	197	Pdxdy0(ikl) = dxHOST * dxHOST ! grid area [m2]
	198	DO klc= 1,mzp
	199	k = mzpp-klc
	200	P_pa_0(ikl,klc) = (psa_DY( ikl)sigma(k) + pt__DY) 1.e3 ! pressure in layer [Pa]
	201	P_za_0(ikl,klc) = Z___DY(ikl,k) ! height of model layer [m]
	202	P_Ta_0(ikl,klc) = Ta__DY(ikl,k) ! grid scale T at time t [K]
	203	P_Qa_0(ikl,klc) = qv__DY(ikl,k) ! grid scale water vapor at time t [kg/kg]
	204	P_Qw_0(ikl,klc) = qw__CM(ikl,k) / (1.0-qw__CM(ikl,k)) ! grid scale Cloud drops at time t [kg/kg]
	205	P_Qi_0(ikl,klc) = qi__CM(ikl,k) / (1.0-qi__CM(ikl,k)) ! grid scale Cloud ice at time t [kg/kg]
	206	P_Ua_0(ikl,klc) = ua__DY(ikl,k) ! grid scale hor. wind u at time t [m/s]
	207	P_Va_0(ikl,klc) = va__DY(ikl,k) ! grid scale hor. wind v at time t [m/s]
	208	P_Wa_0(ikl,klc) = wa__DY(ikl,k) + wa_ANA(ikl,k) ! grid scale vertic.wind at time t [m/s]
	209	END DO
	210
	211	END DO ! ikl = 1,kcolp
	212
	213
	214	! Mass Flux convective Scheme: Bechtold et al. (2001) Convective Parameterization
	215	! -------------------------------------------------------------------------------
	216
	217	! ***************
	218	call CONVECTION( &
	219	& kcolc , mzc , kidia0, kfdia0, kbdia0, ktdia0, &
	220	& pdtCV , Odeep , Oshal , Orset0, Odown0, kIce_0, &
	221	& OsetA0, PTdcv , PTscv , &
	222	& kensbl, &
	223	& P_pa_0, P_za_0, Pdxdy0, &
	224	& P_Ta_0, P_Qa_0, P_Qw_0, P_Qi_0, P_Ua_0, P_Va_0, P_Wa_0, &
	225	& Kstep1, PdTa_1, PdQa_1, PdQw_1, PdQi_1, &
	226	& Pdrr_1, Pdss_1, &
	227	& PuMF_1, PdMF_1, Pfrr_1, Pfss_1, Pcape1, K_CbT1, K_CbB1, &
	228	& OCvTC0, KTCCH0, P_CH_0, PdCH_1)
	229	! ***************
	230
	231
	232
	233	! Mass Flux convective Scheme: products
	234	! -------------------------------------
	235
	236	DO ikl = 1,kcolp
	237	CAPECP(ikl) = Pcape1(ikl)
	238	timeCP(ikl) = Kstep1(ikl) *dt__CP
	239	drr_CP(ikl) = Pdrr_1(ikl)
	240	dss_CP(ikl) = Pdss_1(ikl)
	241
	242	DO klc= 1,mzp
	243	k = mzpp - klc
	244	dpktCP(ikl,k) = PdTa_1(ikl,klc) /ExnrDY(ikl,k)
	245	dqv_CP(ikl,k) = PdQa_1(ikl,klc)
	246	dqw_CP(ikl,k) = PdQw_1(ikl,klc)
	247	dqi_CP(ikl,k) = PdQi_1(ikl,klc)
	248
	249	END DO
	250
	251	END DO
	252
	253
	254	END IF ! mod(iitCV0,jjtCV0).EQ.0 (UPDATE) THEN
	255
	256
	257
	258
	259	! Vertical Integrated Energy and Water Content
	260	! ============================================
	261
	262	! #EW irmx = i_x0
	263	! #EW jrmx = j_y0
	264	! #EW rr_max = 0.0
	265
	266	DO ikl = 1,kcolp
	267	i = ii__AP(ikl)
	268	j = jj__AP(ikl)
	269
	270	! #EW enr0EW(ikl) = 0.0
	271	! #EW wat0EW(ikl) = 0.0
	272
	273	! #EW DO k=1,mzp
	274	! #EW temp_r = pkt_DY(ikl,k)*ExnrDY(ikl,k)
	275	! #EW enr0EW(ikl) = enr0EW(ikl) &
	276	! #EW& +(temp_r &
	277	! #EW& -(qw__CM(ikl,k)+qr__CM(ikl,k)) * Lv_CPd &
	278	! #EW& -(qi__CM(ikl,k)+qs__CM(ikl,k)) * Ls_CPd)*dsigmi(k)
	279	! #EW wat0EW(ikl) = wat0EW(ikl) &
	280	! #EW& +(qv__DY(ikl,k) &
	281	! #EW& + qw__CM(ikl,k)+qr__CM(ikl,k) &
	282	! #EW& + qi__CM(ikl,k)+qs__CM(ikl,k) )*dsigmi(k)
	283	! #EW END DO
	284
	285	! #EW enr0EW(ikl) = enr0EW(ikl) * psa_DY( ikl) * Grav_I
	286	! #EW wat0EW(ikl) = wat0EW(ikl) * psa_DY( ikl) * Grav_I
	287	! .. wat0EW [m] contains implicit factor 1.d3 [kPa-->Pa] /ro_Wat
	288
	289	! #EW energ0 = energ0 - enr0EW(ikl)
	290	! #EW water0 = water0 - wat0EW(ikl)
	291
	292
	293
	294
	295	! Update of Mass Flux convective Tendencies
	296	! =========================================
	297
	298	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	299
	300	IF (timeCP(ikl) .GT. 0.) THEN
	301
	302
	303	! Temporal tendencies on pkt_DY, qv__DY and rainCM
	304	! ------------------------------------------------
	305
	306	DO k=1,mzp
	307	dqv_CP(ikl,k) = min(dqv_CP(ikl,k),(qv__DY(ikl,k)-epsq)/dt__CP)
	308	dqw_CP(ikl,k) = min(dqw_CP(ikl,k),(qw__CM(ikl,k)-epsq)/dt__CP)
	309	dqi_CP(ikl,k) = min(dqi_CP(ikl,k),(qi__CM(ikl,k)-epsq)/dt__CP)
	310	ENDDO
	311
	312	rainCM(ikl) = rainCM(ikl) + drr_CP(ikl) *dt__CP
	313	rainCP(ikl) = rainCP(ikl) + drr_CP(ikl) *dt__CP
	314
	315	snowCM(ikl) = snowCM(ikl) + dss_CP(ikl) *dt__CP
	316	snowCP(ikl) = snowCP(ikl) + dss_CP(ikl) *dt__CP
	317
	318	ELSE
	319
	320	DO k=1,mzp
	321	dpktCP(ikl,k) = 0.00
	322	dqv_CP(ikl,k) = 0.00
	323	dqw_CP(ikl,k) = 0.00
	324	dqi_CP(ikl,k) = 0.00
	325	ENDDO
	326
	327
	328	ENDIF ! {timeCP(ikl) .gt. 0}
	329
	330	timeCP(ikl) = max(timeCP(ikl) - dt__CP,zer0)
	331	! .. ^^^^ remaining time before the end of convection
	332	! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	333
	334
	335
	336
	337	! Vertical Integrated Energy and Water Content
	338	! ============================================
	339
	340	! #EW enr1EW(ikl) = 0.0
	341	! #EW wat1EW(ikl) = 0.0
	342	! #EW watfEW(ikl) =-drr_CP(ikl)
	343
	344	! #EW DO k=1,mz
	345	! #EW temp_r = pkt_DY(ikl,k)*ExnrDY(ikl,k)
	346	! #EW enr1EW(ikl) = enr1EW(ikl) &
	347	! #EW& +(temp_r &
	348	! #EW& -(qw__CM(ikl,k)+qr__CM(ikl,k)) *Lv_CPd &
	349	! #EW& -(qi__CM(ikl,k)+qs__CM(ikl,k)) Ls_CPd)dsigmi(k)
	350	! #EW wat1EW(ikl) = wat1EW(ikl) &
	351	! #EW& +(qv__DY(ikl,k) &
	352	! #EW& + qw__CM(ikl,k)+qr__CM(ikl,k) &
	353	! #EW& + qi__CM(ikl,k)+qs__CM(ikl,k) )*dsigmi(k)
	354	! #EW END DO
	355
	356	! #EW enr1EW(ikl) = enr1EW(ikl) * psa_DY( ikl) *Grav_I &
	357	! #EW& - drr_CP(ikl) *Lv_CPd
	358	! #EW wat1EW(ikl) = wat1EW(ikl) * psa_DY( ikl) *Grav_I
	359	! .. wat1EW [m] contains implicit factor 1.d3 [kPa-->Pa] /rhoWat
	360
	361	! #EW energ0 = energ0 + enr1EW(ikl)
	362	! #EW water0 = water0 + wat1EW(ikl)
	363	! #EW iter_0 = iter_0 + 1
	364
	365
	366
	367
	368	! Vertical Integrated Energy and Water Content: OUTPUT
	369	! ====================================================
	370
	371	! #EW IF (drr_CP(ikl).gt.rr_max) THEN
	372	! #EW rr_max = drr_CP(ikl)
	373	! #EW irmx = i
	374	! #EW jrmx = j
	375	! #EW END IF
	376
	377	END DO
	378
	379	! #EW waterb = wat1EW(irmx,jrmx)-wat0EW(irmx,jrmx)-watfEW(irmx,jrmx)
	380	! #EW write(6,606) it_EXP,enr0EW(irmx,jrmx),1.d3*wat0EW(irmx,jrmx), &
	381	! #EW& irmx,jrmx,enr1EW(irmx,jrmx),1.d3*wat1EW(irmx,jrmx), &
	382	! #EW& 1.d3*watfEW(irmx,jrmx), &
	383	! #EW& 1.d3*waterb , &
	384	! #EW& energ0/iter_0 , water0/iter_0
	385	606 format(i9,' Before CVAj: E0 =',f12.6,' W0 = ',f9.6, &
	386	& /,i5,i4,' After CVAj: E1 =',f12.6,' W1 = ',f9.6, &
	387	& ' W Flux =',f9.6, &
	388	& ' Div(W) =',e9.3, &
	389	& /,9x,' Mean dE =',f12.9,' dW = ',e9.3)
	390
	391
	392
	393
	394	! Incrementation Step Nb
	395	! ======================
	396
	397	iitCV0 = iitCV0 + 1
	398
	399
	400	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	401	! !
	402	! DE-ALLOCATION
	403	! =============
	404
	405	IF (FlagDALLOC) THEN !
	406	deallocate ( wa_ANA ) ! ANAbatic Wind Speed [m/s]
	407	END IF
	408
	409	! Bug corrigé par Martin:
	410	! IF (it_RUN.EQ.1 .OR. FlagDALLOC) THEN !
	411	! deallocate ( wa_ANA ) ! ANAbatic Wind Speed [m/s]
	412	! END IF
	413	! !
	414	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	415
	416
	417
	418	return
	419	end subroutine PHY_Atm_CP_RUN

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

Download in other formats:

Original Format