Index: LMDZ4/trunk/libf/phylmd/cltracrn.F90 =================================================================== --- LMDZ4/trunk/libf/phylmd/cltracrn.F90 (revision 1408) +++ LMDZ4/trunk/libf/phylmd/cltracrn.F90 (revision 1409) @@ -8,4 +8,5 @@ USE dimphy + USE traclmdz_mod, ONLY : id_rn, id_pb IMPLICIT NONE !====================================================================== @@ -20,5 +21,5 @@ !--------------------------------------------------------------------- ! Arguments: -! itr......input-R- le type de traceur 1- Rn 2 - Pb +! itr......input-R- le type de traceur : id_rn(radon), id_pb(plomb) ! dtime....input-R- intervalle du temps (en secondes) ~ pdtphys ! u1lay....input-R- vent u de la premiere couche (m/s) @@ -195,6 +196,6 @@ !-------------------------------------------------------- - IF ( (itr.eq.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60..AND.lat(i).LE.70.).OR. & - (itr.eq.2.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60..AND.lat(i).LE.70.) ) THEN + IF ( (itr.eq.id_rn.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60..AND.lat(i).LE.70.).OR. & + (itr.eq.id_pb.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GE.60..AND.lat(i).LE.70.) ) THEN zx_trs(i) = local_trs(i) zx_a = zx_trs(i) & @@ -221,6 +222,6 @@ !--------------------------------------------------------------- - IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0).OR. & - (itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).LT.-60.)) THEN + IF ( (itr.EQ.id_rn.AND.NINT(masktr(i)).EQ.0).OR. & + (itr.EQ.id_rn.AND.NINT(masktr(i)).EQ.1.AND.lat(i).LT.-60.)) THEN zx_trs(i) = 0. local_trs(i) = 0. @@ -230,6 +231,6 @@ ! des oceans et des continents !-------------------------------------------------------------- - IF ( (itr.EQ.1.AND.NINT(masktr(i)).EQ.0).OR. & - (itr.EQ.1.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GT.70.)) THEN + IF ( (itr.EQ.id_rn.AND.NINT(masktr(i)).EQ.0).OR. & + (itr.EQ.id_rn.AND.NINT(masktr(i)).EQ.1.AND.lat(i).GT.70.)) THEN zx_trs(i) = 0. local_trs(i) = 0. @@ -239,5 +240,5 @@ !-------------------------------------------- - IF (itr.eq.1.AND.NINT(masktr(i)).EQ.0) THEN + IF (itr.eq.id_rn.AND.NINT(masktr(i)).EQ.0) THEN zx_trs(i) = 0. local_trs(i) = 0. Index: LMDZ4/trunk/libf/phylmd/initrrnpb.F90 =================================================================== --- LMDZ4/trunk/libf/phylmd/initrrnpb.F90 (revision 1408) +++ LMDZ4/trunk/libf/phylmd/initrrnpb.F90 (revision 1409) @@ -5,12 +5,11 @@ USE dimphy USE infotrac, ONLY : nbtr + USE traclmdz_mod, ONLY : id_rn, id_pb IMPLICIT NONE !====================================================================== ! Auteur(s): AA + CG (LGGE/CNRS) Date 24-06-94 ! Objet: initialisation des constantes des traceurs -!AA Revison pour le controle avec la temperature du sol -!AA -!AA it = 1 radon ss controle de ts -!AA it = 2 plomb ss controle de ts +! id_rn : identificateur du traceur radon +! id_pb : identificateur du traceur plomb !====================================================================== ! Arguments: @@ -42,61 +41,56 @@ CHARACTER (LEN=80) :: abort_message +! +! Radon it = id_rn +!---------------- + IF (id_rn /= 0) THEN + it = id_rn + s = 1.E4 ! Source: atome par m2 + hsoltr(it) = 0.1 ! Hauteur equivalente du reservoir : + ! 1 m * porosite 0.1 + tautr(it) = 4.765E5 ! Decroissance du radon, secondes + vdeptr(it) = 0. ! Pas de depot sec pour le radon + scavtr(it) = 0. ! Pas de lessivage pour le radon + + WRITE(*,*)'-------------- SOURCE DU RADON ------------------------ ' + WRITE(*,*)'it = ',it + WRITE(*,*)'Source : ', s + WRITE(*,*)'Hauteur equivalente du reservoir de sol: ',hsoltr(it) + WRITE(*,*)'Decroissance (s): ', tautr(it) + WRITE(*,*)'Vitesse de depot sec: ',vdeptr(it) + WRITE(*,*)'Facteur de lessivage: ',scavtr(it) - WRITE(*,*)'PASSAGE initrrnpb ...' + DO i = 1,klon + masktr(i,it) = 0. + IF ( NINT(pctsrf(i,1)) .EQ. 1 ) masktr(i,it) = 1. + fshtr(i,it) = s * masktr(i,it) + END DO + + END IF ! id_rn /= 0 + ! -! Radon it = 1 +! 210Pb it = id_pb !---------------- - IF ( nbtr .LE. 0 ) then - abort_message = '**PHYTRAC:initrrnpb:** nbtr < 0; verifier RN dans traceur.def' - CALL abort_gcm (modname,abort_message,1) - ENDIF - it = 1 - s = 1.E4 ! Source: atome par m2 - hsoltr(it) = 0.1 ! Hauteur equivalente du reservoir : - ! 1 m * porosite 0.1 - tautr(it) = 4.765E5 ! Decroissance du radon, secondes - vdeptr(it) = 0. ! Pas de depot sec pour le radon - scavtr(it) = 0. ! Pas de lessivage pour le radon - - WRITE(*,*)'-------------- SOURCE DU RADON ------------------------ ' - WRITE(*,*)'it = ',it - WRITE(*,*)'Source : ', s - WRITE(*,*)'Hauteur equivalente du reservoir de sol: ',hsoltr(it) - WRITE(*,*)'Decroissance (s): ', tautr(it) - WRITE(*,*)'Vitesse de depot sec: ',vdeptr(it) - WRITE(*,*)'Facteur de lessivage: ',scavtr(it) - - DO i = 1,klon - masktr(i,it) = 0. - IF ( NINT(pctsrf(i,1)) .EQ. 1 ) masktr(i,it) = 1. - fshtr(i,it) = s * masktr(i,it) - END DO -! -! 210Pb it = 2 -!---------------- - IF ( nbtr .LE. 1 ) THEN - abort_message='**PHYTRAC**:initrrnpb:** nbtr <= 1; verifier PB dans traceur.def' - CALL abort_gcm (modname,abort_message,1) - ENDIF - it = 2 - s = 0. ! Pas de source - hsoltr(it) = 10. ! Hauteur equivalente du reservoir - ! a partir duquel le depot Brownien a lieu - tautr(it) = 1.028E9 ! Decroissance du Pb210, secondes - vdeptr(it) = 1.E-3 ! 1 mm/s pour le 210Pb - scavtr(it) = .5 ! Lessivage du Pb210 - DO i = 1,klon - masktr(i,it) = 1. ! Le depot sec peut avoir lieu partout - fshtr(i,it) = s * masktr(i,it) - END DO - WRITE(*,*)'-------------- SOURCE DU PLOMB ------------------------ ' - WRITE(*,*)'it = ',it - WRITE(*,*)'Source : ', s - WRITE(*,*)'Hauteur equivalente du reservoir : ',hsoltr(it) - WRITE(*,*)'Decroissance (s): ', tautr(it) - WRITE(*,*)'Vitesse de depot sec: ',vdeptr(it) - WRITE(*,*)'Facteur de lessivage: ',scavtr(it) - - WRITE(*,*) 'Initialisation RN et PB ok' - + IF (id_pb /= 0) THEN + it = id_pb + s = 0. ! Pas de source + hsoltr(it) = 10. ! Hauteur equivalente du reservoir + ! a partir duquel le depot Brownien a lieu + tautr(it) = 1.028E9 ! Decroissance du Pb210, secondes + vdeptr(it) = 1.E-3 ! 1 mm/s pour le 210Pb + scavtr(it) = .5 ! Lessivage du Pb210 + DO i = 1,klon + masktr(i,it) = 1. ! Le depot sec peut avoir lieu partout + fshtr(i,it) = s * masktr(i,it) + END DO + WRITE(*,*)'-------------- SOURCE DU PLOMB ------------------------ ' + WRITE(*,*)'it = ',it + WRITE(*,*)'Source : ', s + WRITE(*,*)'Hauteur equivalente du reservoir : ',hsoltr(it) + WRITE(*,*)'Decroissance (s): ', tautr(it) + WRITE(*,*)'Vitesse de depot sec: ',vdeptr(it) + WRITE(*,*)'Facteur de lessivage: ',scavtr(it) + + END IF + END SUBROUTINE initrrnpb Index: LMDZ4/trunk/libf/phylmd/radio_decay.F90 =================================================================== --- LMDZ4/trunk/libf/phylmd/radio_decay.F90 (revision 1408) +++ LMDZ4/trunk/libf/phylmd/radio_decay.F90 (revision 1409) @@ -8,10 +8,10 @@ USE dimphy USE infotrac, ONLY : nbtr + USE traclmdz_mod, ONLY : id_rn, id_pb IMPLICIT NONE !----------------------------------------------------------------------- ! Auteur(s): AA + CG (LGGE/CNRS) Date 24-06-94 ! Objet: Calcul de la tendance radioactive des traceurs type radioelements -!CG240694 : Pour un traceur, le radon -!CG161294 : Plus un 2eme traceur, le 210Pb. Le radon decroit en plomb. +! Cas particulier pour le couple radon-plomb : Le radon decroit en plomb !----------------------------------------------------------------------- ! @@ -34,4 +34,5 @@ DO it = 1,nbtr + d_tr(:,:,it) = 0. IF ( radio(it) ) THEN IF (tautr(it) .GT. 0.) THEN @@ -41,18 +42,15 @@ END DO END DO - ELSE - d_tr(:,:,it) = 0. END IF - ELSE - d_tr(:,:,it) = 0. END IF END DO + !------------------------------------------------------- -!CG161294 : Cas particulier radon [it=1] => plomb [it=2] +! Cas particulier radon (id_rn) => plomb (id_pb) !------------------------------------------------------- IF ( rnpb ) THEN DO k = 1,klev DO i = 1,klon - d_tr(i,k,2) = d_tr(i,k,2) - d_tr(i,k,1) + d_tr(i,k,id_pb) = d_tr(i,k,id_pb) - d_tr(i,k,id_rn) ENDDO ENDDO Index: LMDZ4/trunk/libf/phylmd/traclmdz_mod.F90 =================================================================== --- LMDZ4/trunk/libf/phylmd/traclmdz_mod.F90 (revision 1408) +++ LMDZ4/trunk/libf/phylmd/traclmdz_mod.F90 (revision 1409) @@ -34,4 +34,12 @@ !$OMP THREADPRIVATE(trs) + INTEGER,SAVE :: id_aga ! Identification number for tracer : Age of stratospheric air +!$OMP THREADPRIVATE(id_aga) + 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. +!$OMP THREADPRIVATE(lev_1p5km) + + INTEGER,SAVE :: id_rn, id_pb ! Identification number for tracer : radon (Rn222), lead (Pb210) +!$OMP THREADPRIVATE(id_rn, id_pb) + INTEGER,SAVE :: id_be ! Activation et position du traceur Be7 [ id_be=0 -> desactive ] !$OMP THREADPRIVATE(id_be) @@ -51,5 +59,5 @@ ! 0 means no ozone tracer - LOGICAL,SAVE :: rnpb=.TRUE. ! Presence du couple Rn222, Pb210 + LOGICAL,SAVE :: rnpb=.FALSE. ! Presence du couple Rn222, Pb210 !$OMP THREADPRIVATE(rnpb) @@ -143,9 +151,27 @@ ! Recherche des traceurs connus : Be7, O3, CO2,... ! -------------------------------------------- - id_be=0 - id_o3=0 - DO it=1,nbtr - iiq=niadv(it+2) - IF ( tname(iiq) == "BE" .OR. tname(iiq) == "Be" .OR. & + id_rn=0; id_pb=0; id_aga=0; id_be=0; id_o3=0 + DO it=1,nbtr + iiq=niadv(it+2) + IF ( tname(iiq) == "RN" ) THEN + id_rn=it ! radon + ELSE IF ( tname(iiq) == "PB") THEN + id_pb=it ! plomb + ELSE IF ( tname(iiq) == "Aga" .OR. tname(iiq) == "AGA" ) THEN + ! Age of stratospheric air + id_aga=it + radio(id_aga) = .FALSE. + aerosol(id_aga) = .FALSE. + pbl_flg(id_aga) = 0 + + ! Find the first model layer above 1.5km from the surface + IF (klev>=30) THEN + lev_1p5km=6 ! NB! This value is for klev=39 + ELSE IF (klev>=10) THEN + lev_1p5km=5 ! NB! This value is for klev=19 + ELSE + lev_1p5km=klev/2 + END IF + ELSE IF ( tname(iiq) == "BE" .OR. tname(iiq) == "Be" .OR. & tname(iiq) == "BE7" .OR. tname(iiq) == "Be7" ) THEN ! Recherche du Beryllium 7 @@ -165,5 +191,4 @@ id_dry=0 - DO it=1,nbtr iiq=niadv(it+2) @@ -224,12 +249,12 @@ ! Valeurs specifiques pour les traceurs Rn222 et Pb210 ! ---------------------------------------------- - IF (rnpb) THEN - - radio(1)= .TRUE. - radio(2)= .TRUE. - pbl_flg(1) = 0 ! au lieu de clsol=true ! CL au sol calcule - pbl_flg(2) = 0 ! au lieu de clsol=true - - aerosol(2) = .TRUE. ! le Pb est un aerosol + IF ( id_rn/=0 .AND. id_pb/=0 ) THEN + rnpb = .TRUE. + radio(id_rn)= .TRUE. + radio(id_pb)= .TRUE. + pbl_flg(id_rn) = 0 ! au lieu de clsol=true ! CL au sol calcule + pbl_flg(id_pb) = 0 ! au lieu de clsol=true + aerosol(id_rn) = .FALSE. + aerosol(id_pb) = .TRUE. ! le Pb est un aerosol CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr) @@ -480,4 +505,25 @@ WRITE(solsym(it),'(i2)') it END DO + +!================================================================= +! Update tracer : Age of stratospheric air +!================================================================= + IF (id_aga/=0) THEN + + ! Bottom layers + DO k = 1, lev_1p5km + tr_seri(:,k,id_aga) = 0.0 + END DO + + ! Layers above 1.5km + DO k = lev_1p5km+1,klev-1 + tr_seri(:,k,id_aga) = tr_seri(:,k,id_aga) + pdtphys + END DO + + ! Top layer + tr_seri(:,klev,id_aga) = tr_seri(:,klev-1,id_aga) + + END IF + !====================================================================== ! -- Calcul de l'effet de la couche limite -- @@ -497,5 +543,6 @@ DO it=1, nbtr - IF (couchelimite .AND. pbl_flg(it) == 0 ) THEN ! couche limite avec quantite dans le sol calculee + IF (couchelimite .AND. pbl_flg(it) == 0 .AND. (it==id_rn .OR. it==id_pb)) THEN + ! couche limite avec quantite dans le sol calculee CALL cltracrn(it, pdtphys, yu1, yv1, & cdragh, coefh,t_seri,ftsol,pctsrf, &