Changeset 1409 for LMDZ4/trunk/libf/phylmd/traclmdz_mod.F90

Timestamp:

Jul 8, 2010, 1:46:29 PM (14 years ago)

Author:

jghattas

Message:

• Added tracer "Age of stratospheric air", activated by putting Aga in tracer.def.
• The logical rnpb is now controled during execution. rnpb is true only if both tracers RN and PB existe in tracer.def. RN and PB can now be removed from tracer.def
• In tracer.def, the 2 water traceurs (H2Ov and H2Ol) must still be the first 2 tracers. The following tracers have no specific order(RN and PB can now change places). Still a minimum of 3 tracers are required.

File:

• LMDZ4/trunk/libf/phylmd/traclmdz_mod.F90 (modified) (7 diffs)

LMDZ4/trunk/libf/phylmd/traclmdz_mod.F90

@@ -34 +34 @@
 !$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 +59 @@
   ! 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 +151 @@
 ! 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 +191 @@
 
     id_dry=0
-
     DO it=1,nbtr
        iiq=niadv(it+2)
@@ -224 +249 @@
 ! 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 +505 @@
        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 +543 @@
     
     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,  &