!$Id $ SUBROUTINE cltracrn( itr, dtime,u1lay, v1lay, & cdrag,coef,t,ftsol,pctsrf, & tr,trs,paprs,pplay,delp, & masktr,fshtr,hsoltr,tautr,vdeptr, & lat,d_tr,d_trs ) USE dimphy USE traclmdz_mod, ONLY : id_rn, id_pb USE indice_sol_mod USE yomcst_mod_h IMPLICIT NONE !====================================================================== ! Auteur(s): Alex/LMD) date: fev 99 ! inspire de clqh + clvent ! Objet: diffusion verticale de traceurs avec quantite de traceur ds ! le sol ( reservoir de sol de radon ) ! ! note : pour l'instant le traceur dans le sol et le flux sont ! calcules mais ils ne servent que de diagnostiques ! seule la tendance sur le traceur est sortie (d_tr) !--------------------------------------------------------------------- ! Arguments: ! 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) ! v1lay....input-R- vent v de la premiere couche (m/s) ! cdrag....input-R- cdrag ! coef.....input-R- le coefficient d'echange (m**2/s) l>1, valable uniquement pour k entre 2 et klev ! t........input-R- temperature (K) ! paprs....input-R- pression a inter-couche (Pa) ! pplay....input-R- pression au milieu de couche (Pa) ! delp.....input-R- epaisseur de couche (Pa) ! ftsol....input-R- temperature du sol (en Kelvin) ! tr.......input-R- traceurs ! trs......input-R- traceurs dans le sol ! masktr...input-R- Masque reservoir de sol traceur (1 = reservoir) ! fshtr....input-R- Flux surfacique de production dans le sol ! tautr....input-R- Constante de decroissance du traceur ! vdeptr...input-R- Vitesse de depot sec dans la couche brownienne ! hsoltr...input-R- Epaisseur equivalente du reservoir de sol ! lat......input-R- latitude en degree ! d_tr.....output-R- le changement de "tr" ! d_trs....output-R- le changement de "trs" !====================================================================== ! !Entrees INTEGER,INTENT(IN) :: itr REAL,INTENT(IN) :: dtime REAL,DIMENSION(klon),INTENT(IN) :: u1lay, v1lay REAL,DIMENSION(klon),INTENT(IN) :: cdrag REAL,DIMENSION(klon,klev),INTENT(IN) :: coef, t REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol, pctsrf REAL,DIMENSION(klon,klev),INTENT(IN) :: tr REAL,DIMENSION(klon),INTENT(IN) :: trs REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay, delp REAL,DIMENSION(klon),INTENT(IN) :: masktr REAL,DIMENSION(klon),INTENT(IN) :: fshtr REAL,INTENT(IN) :: hsoltr REAL,INTENT(IN) :: tautr REAL,INTENT(IN) :: vdeptr REAL,DIMENSION(klon),INTENT(IN) :: lat !Sorties REAL,DIMENSION(klon,klev),INTENT(OUT) :: d_tr REAL,DIMENSION(klon),INTENT(OUT) :: d_trs ! (diagnostic) traceur ds le sol !Locales REAL,DIMENSION(klon,klev) :: flux_tr ! (diagnostic) flux de traceur INTEGER :: i, k, n, l REAL,DIMENSION(klon) :: rotrhi REAL,DIMENSION(klon,klev) :: zx_coef REAL,DIMENSION(klon) :: zx_buf REAL,DIMENSION(klon,klev) :: zx_ctr REAL,DIMENSION(klon,klev) :: zx_dtr REAL,DIMENSION(klon) :: zx_trs REAL :: zx_a, zx_b REAL,DIMENSION(klon,klev) :: local_tr REAL,DIMENSION(klon) :: local_trs REAL,DIMENSION(klon) :: zts ! champ de temperature du sol REAL,DIMENSION(klon) :: zx_alpha1, zx_alpha2 !====================================================================== !AA Pour l'instant les 4 types de surface ne sont pas pris en compte !AA On fabrique avec zts un champ de temperature de sol !AA que le pondere par la fraction de nature de sol. DO i = 1,klon zts(i) = 0. ENDDO DO n=1,nbsrf DO i = 1,klon zts(i) = zts(i) + ftsol(i,n)*pctsrf(i,n) ENDDO ENDDO DO i = 1,klon rotrhi(i) = RD * zts(i) / hsoltr ENDDO DO k = 1, klev DO i = 1, klon local_tr(i,k) = tr(i,k) ENDDO ENDDO DO i = 1, klon local_trs(i) = trs(i) ENDDO !====================================================================== !AA Attention si dans clmain zx_alf1(i) = 1.0 !AA Il doit y avoir coherence (dc la meme chose ici) DO i = 1, klon !AA zx_alpha1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2)) zx_alpha1(i) = 1.0 zx_alpha2(i) = 1.0 - zx_alpha1(i) ENDDO !====================================================================== DO i = 1, klon zx_coef(i,1) = cdrag(i)*(1.0+SQRT(u1lay(i)**2+v1lay(i)**2)) & *pplay(i,1)/(RD*t(i,1)) zx_coef(i,1) = zx_coef(i,1) * dtime*RG ENDDO DO k = 2, klev DO i = 1, klon zx_coef(i,k) = coef(i,k)*RG/(pplay(i,k-1)-pplay(i,k)) & *(paprs(i,k)*2/(t(i,k)+t(i,k-1))/RD)**2 zx_coef(i,k) = zx_coef(i,k) * dtime*RG ENDDO ENDDO !====================================================================== DO i = 1, klon zx_buf(i) = delp(i,klev) + zx_coef(i,klev) zx_ctr(i,klev) = local_tr(i,klev)*delp(i,klev)/zx_buf(i) zx_dtr(i,klev) = zx_coef(i,klev) / zx_buf(i) ENDDO DO l = klev-1, 2 , -1 DO i = 1, klon zx_buf(i) = delp(i,l)+zx_coef(i,l) & +zx_coef(i,l+1)*(1.-zx_dtr(i,l+1)) zx_ctr(i,l) = ( local_tr(i,l)*delp(i,l) & + zx_coef(i,l+1)*zx_ctr(i,l+1) )/zx_buf(i) zx_dtr(i,l) = zx_coef(i,l) / zx_buf(i) ENDDO ENDDO DO i = 1, klon zx_buf(i) = delp(i,1) + zx_coef(i,2)*(1.-zx_dtr(i,2)) & + masktr(i) * zx_coef(i,1) & *( zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2) ) zx_ctr(i,1) = ( local_tr(i,1)*delp(i,1) & + zx_ctr(i,2) & *(zx_coef(i,2) & - masktr(i) * zx_coef(i,1) & *zx_alpha2(i) ) ) / zx_buf(i) zx_dtr(i,1) = masktr(i) * zx_coef(i,1) / zx_buf(i) ENDDO !====================================================================== ! Calculer d'abord local_trs nouvelle quantite dans le reservoir ! de sol !===================================================================== DO i = 1, klon !------------------------- ! Au dessus des continents !-- ! Le pb peut se deposer partout : vdeptr = 10-3 m/s ! Le Rn est traiter commme une couche Brownienne puisque vdeptr = 0. !------------------------------------------------------------------- IF ( NINT(masktr(i)) .EQ. 1 ) THEN zx_trs(i) = local_trs(i) zx_a = zx_trs(i) & +fshtr(i)*dtime*rotrhi(i) & +rotrhi(i)*masktr(i)*zx_coef(i,1)/RG & *(zx_ctr(i,1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) & +zx_alpha2(i)*zx_ctr(i,2)) ! Pour l'instant, pour aller vite, le depot sec est traite comme une decroissance zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i,1)/RG & * (1.-zx_dtr(i,1) & *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2))) & + dtime / tautr & + dtime * vdeptr / hsoltr zx_trs(i) = zx_a / zx_b local_trs(i) = zx_trs(i) ENDIF !-------------------------------------------------------- ! Si on est entre 60N et 70N on divise par 2 l'emanation !-------------------------------------------------------- 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) & +(fshtr(i)/2.)*dtime*rotrhi(i) & +rotrhi(i)*masktr(i)*zx_coef(i,1)/RG & *(zx_ctr(i,1)*(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2)) & +zx_alpha2(i)*zx_ctr(i,2)) ! zx_b = 1. + rotrhi(i)*masktr(i)*zx_coef(i,1)/RG & * (1.-zx_dtr(i,1) & *(zx_alpha1(i)+zx_alpha2(i)*zx_dtr(i,2))) & + dtime / tautr & + dtime * vdeptr / hsoltr ! zx_trs(i) = zx_a / zx_b local_trs(i) = zx_trs(i) ENDIF !---------------------------------------------- ! Au dessus des oceans et aux hautes latitudes !-- ! au dessous de -60S pas d'emission de radon au dessus ! des oceans et des continents !--------------------------------------------------------------- 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. END IF !-- ! au dessus de 70 N pas d'emission de radon au dessus ! des oceans et des continents !-------------------------------------------------------------- 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. END IF !--------------------------------------------- ! Au dessus des oceans la source est nulle !-------------------------------------------- IF (itr.eq.id_rn.AND.NINT(masktr(i)).EQ.0) THEN zx_trs(i) = 0. local_trs(i) = 0. END IF ENDDO ! sur le i=1,klon ! !====================================================================== ! Une fois on a zx_trs, on peut faire l'iteration !====================================================================== DO i = 1, klon local_tr(i,1) = zx_ctr(i,1)+zx_dtr(i,1)*zx_trs(i) ENDDO DO l = 2, klev DO i = 1, klon local_tr(i,l) = zx_ctr(i,l) + zx_dtr(i,l)*local_tr(i,l-1) ENDDO ENDDO !====================================================================== ! Calcul du flux de traceur (flux_tr): UA/(m**2 s) !====================================================================== DO i = 1, klon flux_tr(i,1) = masktr(i)*zx_coef(i,1)/RG & * (zx_alpha1(i)*local_tr(i,1)+zx_alpha2(i)*local_tr(i,2) & -zx_trs(i)) / dtime ENDDO DO l = 2, klev DO i = 1, klon flux_tr(i,l) = zx_coef(i,l)/RG & * (local_tr(i,l)-local_tr(i,l-1)) / dtime ENDDO ENDDO !====================================================================== ! Calcul des tendances du traceur ds le sol et dans l'atmosphere !====================================================================== DO l = 1, klev DO i = 1, klon d_tr(i,l) = local_tr(i,l) - tr(i,l) ENDDO ENDDO DO i = 1, klon d_trs(i) = local_trs(i) - trs(i) ENDDO END SUBROUTINE cltracrn