Context Navigation

traclmdz_mod.F90 @ 5445

Last change on this file since 5445 was 2408, checked in by Laurent Fairhead, 9 years ago
Merged trunk changes r2298:2396 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 23.4 KB

Rev	Line
[1191]	1	!$Id $
	2	!
	3	MODULE traclmdz_mod
[1750]	4
[1191]	5	!
	6	! In this module all tracers specific to LMDZ are treated. This module is used
	7	! only if running without any other chemestry model as INCA or REPROBUS.
	8	!
[1403]	9	IMPLICIT NONE
	10
[1191]	11	REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: masktr ! Masque reservoir de sol traceur
	12	!$OMP THREADPRIVATE(masktr) ! Masque de l'echange avec la surface (1 = reservoir) ou (possible >= 1 )
	13	REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: fshtr ! Flux surfacique dans le reservoir de sol
	14	!$OMP THREADPRIVATE(fshtr)
	15	REAL,DIMENSION(:),ALLOCATABLE,SAVE :: hsoltr ! Epaisseur equivalente du reservoir de sol
	16	!$OMP THREADPRIVATE(hsoltr)
	17	!
	18	!Radioelements:
	19	!--------------
	20	!
	21	REAL,DIMENSION(:),ALLOCATABLE,SAVE :: tautr ! Constante de decroissance radioactive
	22	!$OMP THREADPRIVATE(tautr)
	23	REAL,DIMENSION(:),ALLOCATABLE,SAVE :: vdeptr ! Vitesse de depot sec dans la couche Brownienne
	24	!$OMP THREADPRIVATE(vdeptr)
	25	REAL,DIMENSION(:),ALLOCATABLE,SAVE :: scavtr ! Coefficient de lessivage
	26	!$OMP THREADPRIVATE(scavtr)
	27	REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: srcbe ! Production du beryllium7 dans l atmosphere (U/s/kgA)
	28	!$OMP THREADPRIVATE(srcbe)
	29
	30	LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: radio ! radio(it) = true => decroisssance radioactive
	31	!$OMP THREADPRIVATE(radio)
	32
	33	REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: trs ! Conc. radon ds le sol
	34	!$OMP THREADPRIVATE(trs)
	35
[1409]	36	INTEGER,SAVE :: id_aga ! Identification number for tracer : Age of stratospheric air
	37	!$OMP THREADPRIVATE(id_aga)
	38	INTEGER,SAVE :: lev_1p5km ! Approximative vertical layer number at 1.5km above surface, used for calculation of the age of air. The result shouldn't be that sensible to the exactness of this value as long as it is in the lower troposphere.
	39	!$OMP THREADPRIVATE(lev_1p5km)
	40
	41	INTEGER,SAVE :: id_rn, id_pb ! Identification number for tracer : radon (Rn222), lead (Pb210)
	42	!$OMP THREADPRIVATE(id_rn, id_pb)
	43
[1191]	44	INTEGER,SAVE :: id_be ! Activation et position du traceur Be7 [ id_be=0 -> desactive ]
	45	!$OMP THREADPRIVATE(id_be)
	46
[1403]	47	INTEGER,SAVE :: id_pcsat, id_pcocsat, id_pcq ! traceurs pseudo-vapeur CL qsat, qsat_oc, q
	48	!$OMP THREADPRIVATE(id_pcsat, id_pcocsat, id_pcq)
[1421]	49	INTEGER,SAVE :: id_pcs0, id_pcos0, id_pcq0 ! traceurs pseudo-vapeur CL qsat, qsat_oc, q
	50	! ! qui ne sont pas transportes par la convection
[1403]	51	!$OMP THREADPRIVATE(id_pcs0, id_pcos0, id_pcq0)
	52
	53	INTEGER, SAVE:: id_o3
[1421]	54	!$OMP THREADPRIVATE(id_o3)
	55	! index of ozone tracer with Cariolle parameterization
	56	! 0 means no ozone tracer
[1403]	57
[1409]	58	LOGICAL,SAVE :: rnpb=.FALSE. ! Presence du couple Rn222, Pb210
[1191]	59	!$OMP THREADPRIVATE(rnpb)
	60
	61
	62	CONTAINS
	63
	64	SUBROUTINE traclmdz_from_restart(trs_in)
	65	! This subroutine initialize the module saved variable trs with values from restart file (startphy.nc).
	66	! This subroutine is called from phyetat0 after the field trs_in has been read.
	67
	68	USE dimphy
[2408]	69	USE infotrac_phy
[1191]	70
	71	! Input argument
	72	REAL,DIMENSION(klon,nbtr), INTENT(IN) :: trs_in
	73
	74	! Local variables
	75	INTEGER :: ierr
	76
	77	! Allocate restart variables trs
	78	ALLOCATE( trs(klon,nbtr), stat=ierr)
[2408]	79	IF (ierr /= 0) CALL abort_physic('traclmdz_from_restart', 'pb in allocation 1',1)
[1191]	80
	81	! Initialize trs with values read from restart file
	82	trs(:,:) = trs_in(:,:)
	83
	84	END SUBROUTINE traclmdz_from_restart
	85
	86
[1665]	87	SUBROUTINE traclmdz_init(pctsrf, xlat, xlon, ftsol, tr_seri, t_seri, pplay, sh, pdtphys, aerosol, lessivage)
[1191]	88	! This subroutine allocates and initialize module variables and control variables.
[1421]	89	! Initialization of the tracers should be done here only for those not found in the restart file.
[1191]	90	USE dimphy
[2408]	91	USE infotrac_phy
[1403]	92	USE regr_pr_comb_coefoz_m, ONLY: alloc_coefoz
	93	USE press_coefoz_m, ONLY: press_coefoz
[1227]	94	USE carbon_cycle_mod, ONLY : carbon_cycle_init, carbon_cycle_tr, carbon_cycle_cpl
[1421]	95	USE mod_grid_phy_lmdz
	96	USE mod_phys_lmdz_para
[1795]	97	USE indice_sol_mod
[2408]	98	USE print_control_mod, ONLY: lunout
[1191]	99
	100	! Input variables
[1227]	101	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
[1665]	102	REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes en degres pour chaque point
	103	REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes en degres pour chaque point
[1227]	104	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol ! Temperature du sol (surf)(Kelvin)
[1403]	105	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri! Concentration Traceur [U/KgA]
	106	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
	107	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
	108	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
[1454]	109	REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
[1227]	110
[1191]	111	! Output variables
	112	LOGICAL,DIMENSION(nbtr), INTENT(OUT) :: aerosol
	113	LOGICAL,INTENT(OUT) :: lessivage
	114
	115	! Local variables
[1403]	116	INTEGER :: ierr, it, iiq, i, k
[1421]	117	REAL, DIMENSION(klon_glo,klev) :: varglo ! variable temporaire sur la grille global
	118	REAL, DIMENSION(klev) :: mintmp, maxtmp
	119	LOGICAL :: zero
[1750]	120	! RomP >>> profil initial Be7
	121	integer ilesfil
	122	parameter (ilesfil=1)
	123	integer irr,kradio
	124	real beryllium(klon,klev)
	125	! profil initial Pb210
	126	integer ilesfil2
	127	parameter (ilesfil2=1)
	128	integer irr2,kradio2
	129	real plomb(klon,klev)
	130	!! RomP <<<
[1191]	131	! --------------------------------------------
	132	! Allocation
	133	! --------------------------------------------
	134	ALLOCATE( scavtr(nbtr), stat=ierr)
[2408]	135	IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 9',1)
[1191]	136	scavtr(:)=1.
	137
	138	ALLOCATE( radio(nbtr), stat=ierr)
[2408]	139	IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 11',1)
[1191]	140	radio(:) = .false. ! Par defaut pas decroissance radioactive
	141
	142	ALLOCATE( masktr(klon,nbtr), stat=ierr)
[2408]	143	IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 2',1)
[1191]	144
	145	ALLOCATE( fshtr(klon,nbtr), stat=ierr)
[2408]	146	IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 3',1)
[1191]	147
	148	ALLOCATE( hsoltr(nbtr), stat=ierr)
[2408]	149	IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 4',1)
[1191]	150
	151	ALLOCATE( tautr(nbtr), stat=ierr)
[2408]	152	IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 5',1)
[1191]	153	tautr(:) = 0.
	154
	155	ALLOCATE( vdeptr(nbtr), stat=ierr)
[2408]	156	IF (ierr /= 0) CALL abort_physic('traclmdz_init', 'pb in allocation 6',1)
[1191]	157	vdeptr(:) = 0.
	158
	159
	160	lessivage = .TRUE.
[1750]	161	!!jyg(20130206) : le choix d activation du lessivage est fait dans phytrac avec iflag_lscav
	162	!! call getin('lessivage',lessivage)
	163	!! if(lessivage) then
	164	!! print*,'lessivage lsc ON'
	165	!! else
	166	!! print*,'lessivage lsc OFF'
	167	!! endif
[1191]	168	aerosol(:) = .FALSE. ! Tous les traceurs sont des gaz par defaut
	169
	170	!
[1403]	171	! Recherche des traceurs connus : Be7, O3, CO2,...
[1191]	172	! --------------------------------------------
[1409]	173	id_rn=0; id_pb=0; id_aga=0; id_be=0; id_o3=0
[1421]	174	id_pcsat=0; id_pcocsat=0; id_pcq=0; id_pcs0=0; id_pcos0=0; id_pcq0=0
[1191]	175	DO it=1,nbtr
[2298]	176	!! iiq=niadv(it+2) ! jyg
	177	iiq=niadv(it+nqo) ! jyg
	178	IF ( tname(iiq) == "RN" ) THEN
[1409]	179	id_rn=it ! radon
	180	ELSE IF ( tname(iiq) == "PB") THEN
	181	id_pb=it ! plomb
[1750]	182	! RomP >>> profil initial de PB210
	183	open (ilesfil2,file='prof.pb210',status='old',iostat=irr2)
	184	IF (irr2 == 0) THEN
	185	read(ilesfil2,*) kradio2
	186	print*,'number of levels for pb210 profile ',kradio2
	187	do k=kradio2,1,-1
	188	read (ilesfil2,*) plomb(:,k)
	189	enddo
	190	close(ilesfil2)
	191	do k=1,klev
	192	do i=1,klon
	193	tr_seri(i,k,id_pb)=plomb(i,k)
	194	!! print*, 'tr_seri',i,k,tr_seri(i,k,id_pb)
	195	enddo
	196	enddo
	197	ELSE
	198	print *, 'Prof.pb210 does not exist: use restart values'
	199	ENDIF
	200	! RomP <<<
[1409]	201	ELSE IF ( tname(iiq) == "Aga" .OR. tname(iiq) == "AGA" ) THEN
	202	! Age of stratospheric air
	203	id_aga=it
	204	radio(id_aga) = .FALSE.
	205	aerosol(id_aga) = .FALSE.
	206	pbl_flg(id_aga) = 0
	207
	208	! Find the first model layer above 1.5km from the surface
	209	IF (klev>=30) THEN
	210	lev_1p5km=6 ! NB! This value is for klev=39
	211	ELSE IF (klev>=10) THEN
	212	lev_1p5km=5 ! NB! This value is for klev=19
	213	ELSE
	214	lev_1p5km=klev/2
	215	END IF
	216	ELSE IF ( tname(iiq) == "BE" .OR. tname(iiq) == "Be" .OR. &
[1191]	217	tname(iiq) == "BE7" .OR. tname(iiq) == "Be7" ) THEN
[1212]	218	! Recherche du Beryllium 7
[1191]	219	id_be=it
	220	ALLOCATE( srcbe(klon,klev) )
	221	radio(id_be) = .TRUE.
	222	aerosol(id_be) = .TRUE. ! le Be est un aerosol
[1750]	223	!jyg le 13/03/2013 ; ajout de pplay en argument de init_be
	224	!!! CALL init_be(pctsrf,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
	225	CALL init_be(pctsrf,pplay,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
[1421]	226	WRITE(lunout,*) 'Initialisation srcBe: OK'
[1750]	227	! RomP >>> profil initial de Be7
	228	open (ilesfil,file='prof.be7',status='old',iostat=irr)
	229	IF (irr == 0) THEN
	230	read(ilesfil,*) kradio
	231	print*,'number of levels for Be7 profile ',kradio
	232	do k=kradio,1,-1
	233	read (ilesfil,*) beryllium(:,k)
	234	enddo
	235	close(ilesfil)
	236	do k=1,klev
	237	do i=1,klon
	238	tr_seri(i,k,id_be)=beryllium(i,k)
	239	!! print*, 'tr_seri',i,k,tr_seri(i,k,id_be)
	240	enddo
	241	enddo
	242	ELSE
	243	print *, 'Prof.Be7 does not exist: use restart values'
	244	ENDIF
	245	! RomP <<<
[1403]	246	ELSE IF (tname(iiq)=="O3" .OR. tname(iiq)=="o3") THEN
	247	! Recherche de l'ozone : parametrization de la chimie par Cariolle
	248	id_o3=it
	249	CALL alloc_coefoz ! allocate ozone coefficients
	250	CALL press_coefoz ! read input pressure levels
[1421]	251	ELSE IF ( tname(iiq) == "pcsat" .OR. tname(iiq) == "Pcsat" ) THEN
	252	id_pcsat=it
	253	ELSE IF ( tname(iiq) == "pcocsat" .OR. tname(iiq) == "Pcocsat" ) THEN
	254	id_pcocsat=it
	255	ELSE IF ( tname(iiq) == "pcq" .OR. tname(iiq) == "Pcq" ) THEN
	256	id_pcq=it
	257	ELSE IF ( tname(iiq) == "pcs0" .OR. tname(iiq) == "Pcs0" ) THEN
	258	id_pcs0=it
	259	conv_flg(it)=0 ! No transport by convection for this tracer
	260	ELSE IF ( tname(iiq) == "pcos0" .OR. tname(iiq) == "Pcos0" ) THEN
	261	id_pcos0=it
	262	conv_flg(it)=0 ! No transport by convection for this tracer
	263	ELSE IF ( tname(iiq) == "pcq0" .OR. tname(iiq) == "Pcq0" ) THEN
	264	id_pcq0=it
	265	conv_flg(it)=0 ! No transport by convection for this tracer
	266	ELSE
	267	WRITE(lunout,*) 'This is an unknown tracer in LMDZ : ', trim(tname(iiq))
[1403]	268	END IF
	269	END DO
	270
[1191]	271	!
	272	! Valeurs specifiques pour les traceurs Rn222 et Pb210
	273	! ----------------------------------------------
[1409]	274	IF ( id_rn/=0 .AND. id_pb/=0 ) THEN
	275	rnpb = .TRUE.
	276	radio(id_rn)= .TRUE.
	277	radio(id_pb)= .TRUE.
	278	pbl_flg(id_rn) = 0 ! au lieu de clsol=true ! CL au sol calcule
	279	pbl_flg(id_pb) = 0 ! au lieu de clsol=true
	280	aerosol(id_rn) = .FALSE.
	281	aerosol(id_pb) = .TRUE. ! le Pb est un aerosol
[1191]	282
	283	CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
	284	END IF
	285
[1227]	286	!
	287	! Initialisation de module carbon_cycle_mod
	288	! ----------------------------------------------
	289	IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
[1454]	290	CALL carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
[1227]	291	END IF
	292
[1421]	293	! Check if all tracers have restart values
	294	! ----------------------------------------------
	295	DO it=1,nbtr
[2298]	296	!! iiq=niadv(it+2) ! jyg
	297	iiq=niadv(it+nqo) ! jyg
[1421]	298	! Test if tracer is zero everywhere.
	299	! Done by master process MPI and master thread OpenMP
	300	CALL gather(tr_seri(:,:,it),varglo)
	301	IF (is_mpi_root .AND. is_omp_root) THEN
	302	mintmp=MINVAL(varglo,dim=1)
	303	maxtmp=MAXVAL(varglo,dim=1)
	304	IF (MINVAL(mintmp,dim=1)==0. .AND. MAXVAL(maxtmp,dim=1)==0.) THEN
	305	! Tracer is zero everywhere
	306	zero=.TRUE.
	307	ELSE
	308	zero=.FALSE.
	309	END IF
[1403]	310	END IF
	311
[1421]	312	! Distribute variable at all processes
	313	CALL bcast(zero)
[1403]	314
[1421]	315	! Initalize tracer that was not found in restart file.
	316	IF (zero) THEN
	317	! The tracer was not found in restart file or it was equal zero everywhere.
	318	WRITE(lunout,*) "The tracer ",trim(tname(iiq))," will be initialized"
	319	IF (it==id_pcsat .OR. it==id_pcq .OR. &
	320	it==id_pcs0 .OR. it==id_pcq0) THEN
	321	tr_seri(:,:,it) = 100.
	322	ELSE IF (it==id_pcocsat .OR. it==id_pcos0) THEN
	323	DO i = 1, klon
	324	IF ( pctsrf (i, is_oce) == 0. ) THEN
	325	tr_seri(i,:,it) = 0.
	326	ELSE
	327	tr_seri(i,:,it) = 100.
	328	END IF
	329	END DO
	330	ELSE
	331	! No specific initialization exist for this tracer
	332	tr_seri(:,:,it) = 0.
	333	END IF
[1403]	334	END IF
[1421]	335	END DO
[1403]	336
[1191]	337	END SUBROUTINE traclmdz_init
	338
[1403]	339	SUBROUTINE traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
	340	cdragh, coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon, couchelimite, sh, &
[1864]	341	rh, pphi, ustar, wstar, ale_bl, ale_wake, zu10m, zv10m, &
[2187]	342	tr_seri, source, d_tr_cl,d_tr_dec, zmasse) !RomP
[1191]	343
	344	USE dimphy
[2408]	345	USE infotrac_phy
[1403]	346	USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
	347	USE o3_chem_m, ONLY: o3_chem
[1227]	348	USE carbon_cycle_mod, ONLY : carbon_cycle, carbon_cycle_tr, carbon_cycle_cpl
[1795]	349	USE indice_sol_mod
	350
[1191]	351	INCLUDE "YOMCST.h"
	352
	353	!==========================================================================
	354	! -- DESCRIPTION DES ARGUMENTS --
	355	!==========================================================================
	356
	357	! Input arguments
	358	!
	359	!Configuration grille,temps:
[1212]	360	INTEGER,INTENT(IN) :: nstep ! nombre d'appels de la physiq
[1403]	361	INTEGER,INTENT(IN) :: julien ! Jour julien
	362	REAL,INTENT(IN) :: gmtime
[1191]	363	REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
	364	REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point
[1403]	365	REAL, INTENT(IN):: xlon(:) ! dim(klon) longitude
[1191]	366
	367	!
	368	!Physique:
	369	!--------
	370	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
	371	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa)
	372	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
[1403]	373	REAL,intent(in):: zmasse (:, :) ! dim(klon,klev) density of air, in kg/m2
[1191]	374
	375
	376	!Couche limite:
	377	!--------------
	378	!
[1750]	379	REAL,DIMENSION(klon),INTENT(IN) :: cdragh ! coeff drag pour T et Q
	380	REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh ! coeff melange CL (m**2/s)
	381	REAL,DIMENSION(klon),INTENT(IN) :: yu1 ! vents au premier niveau
	382	REAL,DIMENSION(klon),INTENT(IN) :: yv1 ! vents au premier niveau
[1191]	383	LOGICAL,INTENT(IN) :: couchelimite
[1750]	384	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
[1707]	385	REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! Humidite relative
	386	REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentie
	387	REAL,DIMENSION(klon),INTENT(IN) :: ustar ! ustar (m/s)
[1864]	388	REAL,DIMENSION(klon),INTENT(IN) :: wstar,ale_bl,ale_wake ! wstar (m/s) and Avail. Lifti. Energ.
[1707]	389	REAL,DIMENSION(klon),INTENT(IN) :: zu10m ! vent zonal 10m (m/s)
	390	REAL,DIMENSION(klon),INTENT(IN) :: zv10m ! vent zonal 10m (m/s)
[1191]	391
	392	! Arguments necessaires pour les sources et puits de traceur:
	393	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol ! Temperature du sol (surf)(Kelvin)
	394	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
	395
	396	! InOutput argument
	397	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
	398
	399	! Output argument
	400	REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: source ! a voir lorsque le flux de surface est prescrit
	401	REAL,DIMENSION(klon,klev,nbtr), INTENT(OUT) :: d_tr_cl ! Td couche limite/traceur
	402
	403	!=======================================================================================
	404	! -- VARIABLES LOCALES TRACEURS --
	405	!=======================================================================================
	406
	407	INTEGER :: i, k, it
[1421]	408	INTEGER :: lmt_pas ! number of time steps of "physics" per day
[1191]	409
	410	REAL,DIMENSION(klon) :: d_trs ! Td dans le reservoir
[1421]	411	REAL,DIMENSION(klon,klev) :: qsat ! pression de la vapeur a saturation
[1191]	412	REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
[1421]	413	REAL :: zrho ! Masse Volumique de l'air KgA/m3
	414	REAL :: amn, amx
[1191]	415	!
	416	!=================================================================
	417	! Ajout de la production en Be7 (Beryllium) srcbe U/s/kgA
	418	!=================================================================
	419	!
	420	IF ( id_be /= 0 ) THEN
	421	DO k = 1, klev
	422	DO i = 1, klon
	423	tr_seri(i,k,id_be) = tr_seri(i,k,id_be)+srcbe(i,k)*pdtphys
	424	END DO
	425	END DO
	426	WRITE(,) 'Ajout srcBe dans tr_seri: OK'
	427	END IF
	428
[1421]	429
	430	!=================================================================
	431	! Update pseudo-vapor tracers
	432	!=================================================================
	433
	434	CALL q_sat(klon*klev,t_seri,pplay,qsat)
	435
[1403]	436	IF ( id_pcsat /= 0 ) THEN
	437	DO k = 1, klev
	438	DO i = 1, klon
[1421]	439	IF ( pplay(i,k).GE.85000.) THEN
	440	tr_seri(i,k,id_pcsat) = qsat(i,k)
	441	ELSE
	442	tr_seri(i,k,id_pcsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcsat))
	443	END IF
[1403]	444	END DO
	445	END DO
	446	END IF
[1191]	447
[1403]	448	IF ( id_pcocsat /= 0 ) THEN
	449	DO k = 1, klev
	450	DO i = 1, klon
[1421]	451	IF ( pplay(i,k).GE.85000.) THEN
	452	IF ( pctsrf (i, is_oce) > 0. ) THEN
	453	tr_seri(i,k,id_pcocsat) = qsat(i,k)
	454	ELSE
	455	tr_seri(i,k,id_pcocsat) = 0.
	456	END IF
	457	ELSE
	458	tr_seri(i,k,id_pcocsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcocsat))
[1403]	459	END IF
	460	END DO
	461	END DO
	462	END IF
	463
	464	IF ( id_pcq /= 0 ) THEN
	465	DO k = 1, klev
	466	DO i = 1, klon
[1421]	467	IF ( pplay(i,k).GE.85000.) THEN
	468	tr_seri(i,k,id_pcq) = sh(i,k)
	469	ELSE
	470	tr_seri(i,k,id_pcq) = MIN (qsat(i,k), tr_seri(i,k,id_pcq))
	471	END IF
[1403]	472	END DO
	473	END DO
	474	END IF
	475
[1421]	476
[1403]	477	IF ( id_pcs0 /= 0 ) THEN
	478	DO k = 1, klev
	479	DO i = 1, klon
[1421]	480	IF ( pplay(i,k).GE.85000.) THEN
	481	tr_seri(i,k,id_pcs0) = qsat(i,k)
	482	ELSE
	483	tr_seri(i,k,id_pcs0) = MIN (qsat(i,k), tr_seri(i,k,id_pcs0))
	484	END IF
[1403]	485	END DO
	486	END DO
	487	END IF
	488
[1421]	489
[1403]	490	IF ( id_pcos0 /= 0 ) THEN
	491	DO k = 1, klev
	492	DO i = 1, klon
[1421]	493	IF ( pplay(i,k).GE.85000.) THEN
	494	IF ( pctsrf (i, is_oce) > 0. ) THEN
	495	tr_seri(i,k,id_pcos0) = qsat(i,k)
	496	ELSE
	497	tr_seri(i,k,id_pcos0) = 0.
	498	END IF
	499	ELSE
	500	tr_seri(i,k,id_pcos0) = MIN (qsat(i,k), tr_seri(i,k,id_pcos0))
	501	END IF
[1403]	502	END DO
	503	END DO
	504	END IF
	505
[1421]	506
[1403]	507	IF ( id_pcq0 /= 0 ) THEN
	508	DO k = 1, klev
	509	DO i = 1, klon
[1421]	510	IF ( pplay(i,k).GE.85000.) THEN
	511	tr_seri(i,k,id_pcq0) = sh(i,k)
	512	ELSE
	513	tr_seri(i,k,id_pcq0) = MIN (qsat(i,k), tr_seri(i,k,id_pcq0))
	514	END IF
[1403]	515	END DO
	516	END DO
	517	END IF
	518
[1409]	519	!=================================================================
	520	! Update tracer : Age of stratospheric air
	521	!=================================================================
	522	IF (id_aga/=0) THEN
	523
	524	! Bottom layers
	525	DO k = 1, lev_1p5km
	526	tr_seri(:,k,id_aga) = 0.0
	527	END DO
	528
	529	! Layers above 1.5km
	530	DO k = lev_1p5km+1,klev-1
	531	tr_seri(:,k,id_aga) = tr_seri(:,k,id_aga) + pdtphys
	532	END DO
	533
	534	! Top layer
	535	tr_seri(:,klev,id_aga) = tr_seri(:,klev-1,id_aga)
	536
	537	END IF
	538
[1191]	539	!======================================================================
	540	! -- Calcul de l'effet de la couche limite --
	541	!======================================================================
	542
	543	IF (couchelimite) THEN
[1212]	544	source(:,:) = 0.0
	545
	546	IF (id_be /=0) THEN
	547	DO i=1, klon
	548	zrho = pplay(i,1)/t_seri(i,1)/RD
	549	source(i,id_be) = - vdeptr(id_be)tr_seri(i,1,id_be)zrho
	550	END DO
	551	END IF
	552
[1191]	553	END IF
	554
	555	DO it=1, nbtr
[1409]	556	IF (couchelimite .AND. pbl_flg(it) == 0 .AND. (it==id_rn .OR. it==id_pb)) THEN
	557	! couche limite avec quantite dans le sol calculee
[1191]	558	CALL cltracrn(it, pdtphys, yu1, yv1, &
	559	cdragh, coefh,t_seri,ftsol,pctsrf, &
	560	tr_seri(:,:,it),trs(:,it), &
[1403]	561	paprs, pplay, zmasse * rg, &
[1191]	562	masktr(:,it),fshtr(:,it),hsoltr(it),&
	563	tautr(it),vdeptr(it), &
	564	xlat,d_tr_cl(:,:,it),d_trs)
	565
	566	DO k = 1, klev
	567	DO i = 1, klon
	568	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
	569	END DO
	570	END DO
	571
	572	! Traceur dans le reservoir sol
	573	DO i = 1, klon
	574	trs(i,it) = trs(i,it) + d_trs(i)
	575	END DO
	576	END IF
	577	END DO
[1403]	578
	579
[1191]	580	!======================================================================
	581	! Calcul de l'effet du puits radioactif
	582	!======================================================================
	583	CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
[1421]	584
[1191]	585	DO it=1,nbtr
[1421]	586	WRITE(solsym(it),'(i2)') it
	587	END DO
	588
	589	DO it=1,nbtr
[1191]	590	IF(radio(it)) then
	591	DO k = 1, klev
	592	DO i = 1, klon
	593	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
	594	END DO
	595	END DO
	596	CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//solsym(it))
	597	END IF
	598	END DO
	599
[1227]	600	!======================================================================
[1403]	601	! Parameterization of ozone chemistry
	602	!======================================================================
	603
	604	IF (id_o3 /= 0) then
	605	lmt_pas = NINT(86400./pdtphys)
	606	IF (MOD(nstep - 1, lmt_pas) == 0) THEN
	607	! Once per day, update the coefficients for ozone chemistry:
	608	CALL regr_pr_comb_coefoz(julien, xlat, paprs, pplay)
	609	END IF
	610	CALL o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, xlat, &
	611	xlon, tr_seri(:, :, id_o3))
	612	END IF
	613
	614	!======================================================================
[1227]	615	! Calcul de cycle de carbon
	616	!======================================================================
	617	IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
[1454]	618	CALL carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
[1227]	619	END IF
	620
[1191]	621	END SUBROUTINE traclmdz
	622
	623
	624	SUBROUTINE traclmdz_to_restart(trs_out)
	625	! This subroutine is called from phyredem.F where the module
	626	! variable trs is written to restart file (restartphy.nc)
	627	USE dimphy
[2408]	628	USE infotrac_phy
[1191]	629
	630	REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out
[1203]	631	INTEGER :: ierr
[1212]	632
	633	IF ( ALLOCATED(trs) ) THEN
	634	trs_out(:,:) = trs(:,:)
[1203]	635	ELSE
[1212]	636	! No previous allocate of trs. This is the case for create_etat0_limit.
	637	trs_out(:,:) = 0.0
[1203]	638	END IF
[1212]	639
[1191]	640	END SUBROUTINE traclmdz_to_restart
	641
	642
	643	END MODULE traclmdz_mod

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

Download in other formats:

Original Format