Context Navigation

traclmdz_mod.F90 @ 2261

Last change on this file since 2261 was 2180, checked in by acozic, 10 years ago
After checking I reload commit 2169 and 2170
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.0 KB

Rev	Line
[1191]	1	!$Id $
	2	!
	3	MODULE traclmdz_mod
[1742]	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
	69	USE infotrac
	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)
[1212]	79	IF (ierr /= 0) CALL abort_gcm('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
[1579]	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
	91	USE infotrac
[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
[1785]	97	USE indice_sol_mod
[1191]	98
[1421]	99	INCLUDE "iniprint.h"
[1191]	100	! Input variables
[1227]	101	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
[1579]	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
[1742]	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)
[1212]	135	IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 9',1)
[1191]	136	scavtr(:)=1.
	137
	138	ALLOCATE( radio(nbtr), stat=ierr)
[1212]	139	IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 11',1)
[1191]	140	radio(:) = .false. ! Par defaut pas decroissance radioactive
	141
	142	ALLOCATE( masktr(klon,nbtr), stat=ierr)
[1212]	143	IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 2',1)
[1191]	144
	145	ALLOCATE( fshtr(klon,nbtr), stat=ierr)
[1212]	146	IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 3',1)
[1191]	147
	148	ALLOCATE( hsoltr(nbtr), stat=ierr)
[1212]	149	IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 4',1)
[1191]	150
	151	ALLOCATE( tautr(nbtr), stat=ierr)
[1212]	152	IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 5',1)
[1191]	153	tautr(:) = 0.
	154
	155	ALLOCATE( vdeptr(nbtr), stat=ierr)
[1212]	156	IF (ierr /= 0) CALL abort_gcm('traclmdz_init', 'pb in allocation 6',1)
[1191]	157	vdeptr(:) = 0.
	158
	159
	160	lessivage = .TRUE.
[1742]	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
	176	iiq=niadv(it+2)
[1409]	177	IF ( tname(iiq) == "RN" ) THEN
	178	id_rn=it ! radon
	179	ELSE IF ( tname(iiq) == "PB") THEN
	180	id_pb=it ! plomb
[1742]	181	! RomP >>> profil initial de PB210
	182	open (ilesfil2,file='prof.pb210',status='old',iostat=irr2)
	183	IF (irr2 == 0) THEN
	184	read(ilesfil2,*) kradio2
	185	print*,'number of levels for pb210 profile ',kradio2
	186	do k=kradio2,1,-1
	187	read (ilesfil2,*) plomb(:,k)
	188	enddo
	189	close(ilesfil2)
	190	do k=1,klev
	191	do i=1,klon
	192	tr_seri(i,k,id_pb)=plomb(i,k)
	193	!! print*, 'tr_seri',i,k,tr_seri(i,k,id_pb)
	194	enddo
	195	enddo
	196	ELSE
	197	print *, 'Prof.pb210 does not exist: use restart values'
	198	ENDIF
	199	! RomP <<<
[1409]	200	ELSE IF ( tname(iiq) == "Aga" .OR. tname(iiq) == "AGA" ) THEN
	201	! Age of stratospheric air
	202	id_aga=it
	203	radio(id_aga) = .FALSE.
	204	aerosol(id_aga) = .FALSE.
	205	pbl_flg(id_aga) = 0
	206
	207	! Find the first model layer above 1.5km from the surface
	208	IF (klev>=30) THEN
	209	lev_1p5km=6 ! NB! This value is for klev=39
	210	ELSE IF (klev>=10) THEN
	211	lev_1p5km=5 ! NB! This value is for klev=19
	212	ELSE
	213	lev_1p5km=klev/2
	214	END IF
	215	ELSE IF ( tname(iiq) == "BE" .OR. tname(iiq) == "Be" .OR. &
[1191]	216	tname(iiq) == "BE7" .OR. tname(iiq) == "Be7" ) THEN
[1212]	217	! Recherche du Beryllium 7
[1191]	218	id_be=it
	219	ALLOCATE( srcbe(klon,klev) )
	220	radio(id_be) = .TRUE.
	221	aerosol(id_be) = .TRUE. ! le Be est un aerosol
[1742]	222	!jyg le 13/03/2013 ; ajout de pplay en argument de init_be
	223	!!! CALL init_be(pctsrf,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
	224	CALL init_be(pctsrf,pplay,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
[1421]	225	WRITE(lunout,*) 'Initialisation srcBe: OK'
[1742]	226	! RomP >>> profil initial de Be7
	227	open (ilesfil,file='prof.be7',status='old',iostat=irr)
	228	IF (irr == 0) THEN
	229	read(ilesfil,*) kradio
	230	print*,'number of levels for Be7 profile ',kradio
	231	do k=kradio,1,-1
	232	read (ilesfil,*) beryllium(:,k)
	233	enddo
	234	close(ilesfil)
	235	do k=1,klev
	236	do i=1,klon
	237	tr_seri(i,k,id_be)=beryllium(i,k)
	238	!! print*, 'tr_seri',i,k,tr_seri(i,k,id_be)
	239	enddo
	240	enddo
	241	ELSE
	242	print *, 'Prof.Be7 does not exist: use restart values'
	243	ENDIF
	244	! RomP <<<
[1403]	245	ELSE IF (tname(iiq)=="O3" .OR. tname(iiq)=="o3") THEN
	246	! Recherche de l'ozone : parametrization de la chimie par Cariolle
	247	id_o3=it
	248	CALL alloc_coefoz ! allocate ozone coefficients
	249	CALL press_coefoz ! read input pressure levels
[1421]	250	ELSE IF ( tname(iiq) == "pcsat" .OR. tname(iiq) == "Pcsat" ) THEN
	251	id_pcsat=it
	252	ELSE IF ( tname(iiq) == "pcocsat" .OR. tname(iiq) == "Pcocsat" ) THEN
	253	id_pcocsat=it
	254	ELSE IF ( tname(iiq) == "pcq" .OR. tname(iiq) == "Pcq" ) THEN
	255	id_pcq=it
	256	ELSE IF ( tname(iiq) == "pcs0" .OR. tname(iiq) == "Pcs0" ) THEN
	257	id_pcs0=it
	258	conv_flg(it)=0 ! No transport by convection for this tracer
	259	ELSE IF ( tname(iiq) == "pcos0" .OR. tname(iiq) == "Pcos0" ) THEN
	260	id_pcos0=it
	261	conv_flg(it)=0 ! No transport by convection for this tracer
	262	ELSE IF ( tname(iiq) == "pcq0" .OR. tname(iiq) == "Pcq0" ) THEN
	263	id_pcq0=it
	264	conv_flg(it)=0 ! No transport by convection for this tracer
	265	ELSE
	266	WRITE(lunout,*) 'This is an unknown tracer in LMDZ : ', trim(tname(iiq))
[1403]	267	END IF
	268	END DO
	269
[1191]	270	!
	271	! Valeurs specifiques pour les traceurs Rn222 et Pb210
	272	! ----------------------------------------------
[1409]	273	IF ( id_rn/=0 .AND. id_pb/=0 ) THEN
	274	rnpb = .TRUE.
	275	radio(id_rn)= .TRUE.
	276	radio(id_pb)= .TRUE.
	277	pbl_flg(id_rn) = 0 ! au lieu de clsol=true ! CL au sol calcule
	278	pbl_flg(id_pb) = 0 ! au lieu de clsol=true
	279	aerosol(id_rn) = .FALSE.
	280	aerosol(id_pb) = .TRUE. ! le Pb est un aerosol
[1191]	281
	282	CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr)
	283	END IF
	284
[1227]	285	!
	286	! Initialisation de module carbon_cycle_mod
	287	! ----------------------------------------------
	288	IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
[1454]	289	CALL carbon_cycle_init(tr_seri, pdtphys, aerosol, radio)
[1227]	290	END IF
	291
[1421]	292	! Check if all tracers have restart values
	293	! ----------------------------------------------
	294	DO it=1,nbtr
	295	iiq=niadv(it+2)
	296	! Test if tracer is zero everywhere.
	297	! Done by master process MPI and master thread OpenMP
	298	CALL gather(tr_seri(:,:,it),varglo)
	299	IF (is_mpi_root .AND. is_omp_root) THEN
	300	mintmp=MINVAL(varglo,dim=1)
	301	maxtmp=MAXVAL(varglo,dim=1)
	302	IF (MINVAL(mintmp,dim=1)==0. .AND. MAXVAL(maxtmp,dim=1)==0.) THEN
	303	! Tracer is zero everywhere
	304	zero=.TRUE.
	305	ELSE
	306	zero=.FALSE.
	307	END IF
[1403]	308	END IF
	309
[1421]	310	! Distribute variable at all processes
	311	CALL bcast(zero)
[1403]	312
[1421]	313	! Initalize tracer that was not found in restart file.
	314	IF (zero) THEN
	315	! The tracer was not found in restart file or it was equal zero everywhere.
	316	WRITE(lunout,*) "The tracer ",trim(tname(iiq))," will be initialized"
	317	IF (it==id_pcsat .OR. it==id_pcq .OR. &
	318	it==id_pcs0 .OR. it==id_pcq0) THEN
	319	tr_seri(:,:,it) = 100.
	320	ELSE IF (it==id_pcocsat .OR. it==id_pcos0) THEN
	321	DO i = 1, klon
	322	IF ( pctsrf (i, is_oce) == 0. ) THEN
	323	tr_seri(i,:,it) = 0.
	324	ELSE
	325	tr_seri(i,:,it) = 100.
	326	END IF
	327	END DO
	328	ELSE
	329	! No specific initialization exist for this tracer
	330	tr_seri(:,:,it) = 0.
	331	END IF
[1403]	332	END IF
[1421]	333	END DO
[1403]	334
[1191]	335	END SUBROUTINE traclmdz_init
	336
[1403]	337	SUBROUTINE traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
	338	cdragh, coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon, couchelimite, sh, &
[1813]	339	rh, pphi, ustar, wstar, ale_bl, ale_wake, zu10m, zv10m, &
[2180]	340	tr_seri, source, d_tr_cl,d_tr_dec, zmasse) !RomP
[1191]	341
	342	USE dimphy
	343	USE infotrac
[1403]	344	USE regr_pr_comb_coefoz_m, ONLY: regr_pr_comb_coefoz
	345	USE o3_chem_m, ONLY: o3_chem
[1227]	346	USE carbon_cycle_mod, ONLY : carbon_cycle, carbon_cycle_tr, carbon_cycle_cpl
[1785]	347	USE indice_sol_mod
	348
[1191]	349	INCLUDE "YOMCST.h"
	350
	351	!==========================================================================
	352	! -- DESCRIPTION DES ARGUMENTS --
	353	!==========================================================================
	354
	355	! Input arguments
	356	!
	357	!Configuration grille,temps:
[1212]	358	INTEGER,INTENT(IN) :: nstep ! nombre d'appels de la physiq
[1403]	359	INTEGER,INTENT(IN) :: julien ! Jour julien
	360	REAL,INTENT(IN) :: gmtime
[1191]	361	REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
	362	REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point
[1403]	363	REAL, INTENT(IN):: xlon(:) ! dim(klon) longitude
[1191]	364
	365	!
	366	!Physique:
	367	!--------
	368	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
	369	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa)
	370	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
[1403]	371	REAL,intent(in):: zmasse (:, :) ! dim(klon,klev) density of air, in kg/m2
[1191]	372
	373
	374	!Couche limite:
	375	!--------------
	376	!
[1742]	377	REAL,DIMENSION(klon),INTENT(IN) :: cdragh ! coeff drag pour T et Q
	378	REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh ! coeff melange CL (m**2/s)
	379	REAL,DIMENSION(klon),INTENT(IN) :: yu1 ! vents au premier niveau
	380	REAL,DIMENSION(klon),INTENT(IN) :: yv1 ! vents au premier niveau
[1191]	381	LOGICAL,INTENT(IN) :: couchelimite
[1742]	382	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
[1670]	383	REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! Humidite relative
	384	REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentie
	385	REAL,DIMENSION(klon),INTENT(IN) :: ustar ! ustar (m/s)
[1813]	386	REAL,DIMENSION(klon),INTENT(IN) :: wstar,ale_bl,ale_wake ! wstar (m/s) and Avail. Lifti. Energ.
[1670]	387	REAL,DIMENSION(klon),INTENT(IN) :: zu10m ! vent zonal 10m (m/s)
	388	REAL,DIMENSION(klon),INTENT(IN) :: zv10m ! vent zonal 10m (m/s)
[1191]	389
	390	! Arguments necessaires pour les sources et puits de traceur:
	391	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol ! Temperature du sol (surf)(Kelvin)
	392	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol)
	393
	394	! InOutput argument
	395	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
	396
	397	! Output argument
	398	REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: source ! a voir lorsque le flux de surface est prescrit
	399	REAL,DIMENSION(klon,klev,nbtr), INTENT(OUT) :: d_tr_cl ! Td couche limite/traceur
	400
	401	!=======================================================================================
	402	! -- VARIABLES LOCALES TRACEURS --
	403	!=======================================================================================
	404
	405	INTEGER :: i, k, it
[1421]	406	INTEGER :: lmt_pas ! number of time steps of "physics" per day
[1191]	407
	408	REAL,DIMENSION(klon) :: d_trs ! Td dans le reservoir
[1421]	409	REAL,DIMENSION(klon,klev) :: qsat ! pression de la vapeur a saturation
[1191]	410	REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive
[1421]	411	REAL :: zrho ! Masse Volumique de l'air KgA/m3
	412	REAL :: amn, amx
[1191]	413	!
	414	!=================================================================
	415	! Ajout de la production en Be7 (Beryllium) srcbe U/s/kgA
	416	!=================================================================
	417	!
	418	IF ( id_be /= 0 ) THEN
	419	DO k = 1, klev
	420	DO i = 1, klon
	421	tr_seri(i,k,id_be) = tr_seri(i,k,id_be)+srcbe(i,k)*pdtphys
	422	END DO
	423	END DO
	424	WRITE(,) 'Ajout srcBe dans tr_seri: OK'
	425	END IF
	426
[1421]	427
	428	!=================================================================
	429	! Update pseudo-vapor tracers
	430	!=================================================================
	431
	432	CALL q_sat(klon*klev,t_seri,pplay,qsat)
	433
[1403]	434	IF ( id_pcsat /= 0 ) THEN
	435	DO k = 1, klev
	436	DO i = 1, klon
[1421]	437	IF ( pplay(i,k).GE.85000.) THEN
	438	tr_seri(i,k,id_pcsat) = qsat(i,k)
	439	ELSE
	440	tr_seri(i,k,id_pcsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcsat))
	441	END IF
[1403]	442	END DO
	443	END DO
	444	END IF
[1191]	445
[1403]	446	IF ( id_pcocsat /= 0 ) THEN
	447	DO k = 1, klev
	448	DO i = 1, klon
[1421]	449	IF ( pplay(i,k).GE.85000.) THEN
	450	IF ( pctsrf (i, is_oce) > 0. ) THEN
	451	tr_seri(i,k,id_pcocsat) = qsat(i,k)
	452	ELSE
	453	tr_seri(i,k,id_pcocsat) = 0.
	454	END IF
	455	ELSE
	456	tr_seri(i,k,id_pcocsat) = MIN (qsat(i,k), tr_seri(i,k,id_pcocsat))
[1403]	457	END IF
	458	END DO
	459	END DO
	460	END IF
	461
	462	IF ( id_pcq /= 0 ) THEN
	463	DO k = 1, klev
	464	DO i = 1, klon
[1421]	465	IF ( pplay(i,k).GE.85000.) THEN
	466	tr_seri(i,k,id_pcq) = sh(i,k)
	467	ELSE
	468	tr_seri(i,k,id_pcq) = MIN (qsat(i,k), tr_seri(i,k,id_pcq))
	469	END IF
[1403]	470	END DO
	471	END DO
	472	END IF
	473
[1421]	474
[1403]	475	IF ( id_pcs0 /= 0 ) THEN
	476	DO k = 1, klev
	477	DO i = 1, klon
[1421]	478	IF ( pplay(i,k).GE.85000.) THEN
	479	tr_seri(i,k,id_pcs0) = qsat(i,k)
	480	ELSE
	481	tr_seri(i,k,id_pcs0) = MIN (qsat(i,k), tr_seri(i,k,id_pcs0))
	482	END IF
[1403]	483	END DO
	484	END DO
	485	END IF
	486
[1421]	487
[1403]	488	IF ( id_pcos0 /= 0 ) THEN
	489	DO k = 1, klev
	490	DO i = 1, klon
[1421]	491	IF ( pplay(i,k).GE.85000.) THEN
	492	IF ( pctsrf (i, is_oce) > 0. ) THEN
	493	tr_seri(i,k,id_pcos0) = qsat(i,k)
	494	ELSE
	495	tr_seri(i,k,id_pcos0) = 0.
	496	END IF
	497	ELSE
	498	tr_seri(i,k,id_pcos0) = MIN (qsat(i,k), tr_seri(i,k,id_pcos0))
	499	END IF
[1403]	500	END DO
	501	END DO
	502	END IF
	503
[1421]	504
[1403]	505	IF ( id_pcq0 /= 0 ) THEN
	506	DO k = 1, klev
	507	DO i = 1, klon
[1421]	508	IF ( pplay(i,k).GE.85000.) THEN
	509	tr_seri(i,k,id_pcq0) = sh(i,k)
	510	ELSE
	511	tr_seri(i,k,id_pcq0) = MIN (qsat(i,k), tr_seri(i,k,id_pcq0))
	512	END IF
[1403]	513	END DO
	514	END DO
	515	END IF
	516
[1409]	517	!=================================================================
	518	! Update tracer : Age of stratospheric air
	519	!=================================================================
	520	IF (id_aga/=0) THEN
	521
	522	! Bottom layers
	523	DO k = 1, lev_1p5km
	524	tr_seri(:,k,id_aga) = 0.0
	525	END DO
	526
	527	! Layers above 1.5km
	528	DO k = lev_1p5km+1,klev-1
	529	tr_seri(:,k,id_aga) = tr_seri(:,k,id_aga) + pdtphys
	530	END DO
	531
	532	! Top layer
	533	tr_seri(:,klev,id_aga) = tr_seri(:,klev-1,id_aga)
	534
	535	END IF
	536
[1191]	537	!======================================================================
	538	! -- Calcul de l'effet de la couche limite --
	539	!======================================================================
	540
	541	IF (couchelimite) THEN
[1212]	542	source(:,:) = 0.0
	543
	544	IF (id_be /=0) THEN
	545	DO i=1, klon
	546	zrho = pplay(i,1)/t_seri(i,1)/RD
	547	source(i,id_be) = - vdeptr(id_be)tr_seri(i,1,id_be)zrho
	548	END DO
	549	END IF
	550
[1191]	551	END IF
	552
	553	DO it=1, nbtr
[1409]	554	IF (couchelimite .AND. pbl_flg(it) == 0 .AND. (it==id_rn .OR. it==id_pb)) THEN
	555	! couche limite avec quantite dans le sol calculee
[1191]	556	CALL cltracrn(it, pdtphys, yu1, yv1, &
	557	cdragh, coefh,t_seri,ftsol,pctsrf, &
	558	tr_seri(:,:,it),trs(:,it), &
[1403]	559	paprs, pplay, zmasse * rg, &
[1191]	560	masktr(:,it),fshtr(:,it),hsoltr(it),&
	561	tautr(it),vdeptr(it), &
	562	xlat,d_tr_cl(:,:,it),d_trs)
	563
	564	DO k = 1, klev
	565	DO i = 1, klon
	566	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
	567	END DO
	568	END DO
	569
	570	! Traceur dans le reservoir sol
	571	DO i = 1, klon
	572	trs(i,it) = trs(i,it) + d_trs(i)
	573	END DO
	574	END IF
	575	END DO
[1403]	576
	577
[1191]	578	!======================================================================
	579	! Calcul de l'effet du puits radioactif
	580	!======================================================================
	581	CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec)
[1421]	582
[1191]	583	DO it=1,nbtr
[1421]	584	WRITE(solsym(it),'(i2)') it
	585	END DO
	586
	587	DO it=1,nbtr
[1191]	588	IF(radio(it)) then
	589	DO k = 1, klev
	590	DO i = 1, klon
	591	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
	592	END DO
	593	END DO
	594	CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//solsym(it))
	595	END IF
	596	END DO
	597
[1227]	598	!======================================================================
[1403]	599	! Parameterization of ozone chemistry
	600	!======================================================================
	601
	602	IF (id_o3 /= 0) then
	603	lmt_pas = NINT(86400./pdtphys)
	604	IF (MOD(nstep - 1, lmt_pas) == 0) THEN
	605	! Once per day, update the coefficients for ozone chemistry:
	606	CALL regr_pr_comb_coefoz(julien, xlat, paprs, pplay)
	607	END IF
	608	CALL o3_chem(julien, gmtime, t_seri, zmasse, pdtphys, xlat, &
	609	xlon, tr_seri(:, :, id_o3))
	610	END IF
	611
	612	!======================================================================
[1227]	613	! Calcul de cycle de carbon
	614	!======================================================================
	615	IF (carbon_cycle_tr .OR. carbon_cycle_cpl) THEN
[1454]	616	CALL carbon_cycle(nstep, pdtphys, pctsrf, tr_seri, source)
[1227]	617	END IF
	618
[1191]	619	END SUBROUTINE traclmdz
	620
	621
	622	SUBROUTINE traclmdz_to_restart(trs_out)
	623	! This subroutine is called from phyredem.F where the module
	624	! variable trs is written to restart file (restartphy.nc)
	625	USE dimphy
	626	USE infotrac
	627
	628	REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out
[1203]	629	INTEGER :: ierr
[1212]	630
	631	IF ( ALLOCATED(trs) ) THEN
	632	trs_out(:,:) = trs(:,:)
[1203]	633	ELSE
[1212]	634	! No previous allocate of trs. This is the case for create_etat0_limit.
	635	trs_out(:,:) = 0.0
[1203]	636	END IF
[1212]	637
[1191]	638	END SUBROUTINE traclmdz_to_restart
	639
	640
	641	END MODULE traclmdz_mod

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

Download in other formats:

Original Format