Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/infotrac.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/dyn3d/infotrac.F90 (revision 1190) +++ LMDZ4/branches/LMDZ4-dev/libf/dyn3d/infotrac.F90 (revision 1191) @@ -21,8 +21,11 @@ INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: niadv ! equivalent dyn / physique -! Variables for INCA +! conv_flg(it)=0 : convection desactivated for tracer number it INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: conv_flg +! pbl_flg(it)=0 : boundary layer diffusion desactivaded for tracer number it INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: pbl_flg + CHARACTER(len=4),SAVE :: type_trac + CONTAINS @@ -80,4 +83,11 @@ descrq(20)='SLP' descrq(30)='PRA' + + + IF (config_inca=='none') THEN + type_trac='lmdz' + ELSE + type_trac='inca' + END IF !----------------------------------------------------------------------- @@ -87,5 +97,5 @@ ! !----------------------------------------------------------------------- - IF (config_inca == 'none') THEN + IF (type_trac == 'lmdz') THEN OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr) IF(ierr.EQ.0) THEN @@ -109,13 +119,11 @@ END IF ! -! Allocate variables depending on nqtrue +! Allocate variables depending on nqtrue and nbtr ! ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue)) - - IF (config_inca /= 'none') THEN - ! Varaibles only needed in case of INCA - ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), tracnam(nbtr)) - END IF - + ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), tracnam(nbtr)) + conv_flg(:) = 1 ! convection activated for all tracers + pbl_flg(:) = 1 ! boundary layer activated for all tracers + !----------------------------------------------------------------------- ! 2) Choix des schemas d'advection pour l'eau et les traceurs @@ -144,5 +152,5 @@ ! Get choice of advection schema from file tracer.def or from INCA !--------------------------------------------------------------------- - IF (config_inca == 'none') THEN + IF (type_trac == 'lmdz') THEN IF(ierr.EQ.0) THEN ! Continue to read tracer.def @@ -172,5 +180,5 @@ END DO - ELSE ! config_inca='aero' ou 'chem' + ELSE ! type_trac=inca : config_inca='aero' ou 'chem' ! le module de chimie fournit les noms des traceurs ! et les schemas d'advection associes. @@ -191,5 +199,5 @@ END DO - END IF ! config_inca + END IF ! type_trac !----------------------------------------------------------------------- @@ -319,5 +327,5 @@ ! DEALLOCATE(tnom_0, hadv, vadv) - IF (config_inca /= 'none') DEALLOCATE(tracnam) + DEALLOCATE(tracnam) 999 FORMAT (i2,1x,i2,1x,a8) Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/infotrac.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/infotrac.F90 (revision 1190) +++ LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/infotrac.F90 (revision 1191) @@ -21,8 +21,11 @@ INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: niadv ! equivalent dyn / physique -! Variables for INCA +! conv_flg(it)=0 : convection desactivated for tracer number it INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: conv_flg +! pbl_flg(it)=0 : boundary layer diffusion desactivaded for tracer number it INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: pbl_flg + CHARACTER(len=4),SAVE :: type_trac + CONTAINS @@ -80,4 +83,11 @@ descrq(20)='SLP' descrq(30)='PRA' + + + IF (config_inca=='none') THEN + type_trac='lmdz' + ELSE + type_trac='inca' + END IF !----------------------------------------------------------------------- @@ -87,5 +97,5 @@ ! !----------------------------------------------------------------------- - IF (config_inca == 'none') THEN + IF (type_trac == 'lmdz') THEN OPEN(90,file='traceur.def',form='formatted',status='old', iostat=ierr) IF(ierr.EQ.0) THEN @@ -109,13 +119,11 @@ END IF ! -! Allocate variables depending on nqtrue +! Allocate variables depending on nqtrue and nbtr ! ALLOCATE(tnom_0(nqtrue), hadv(nqtrue), vadv(nqtrue)) - - IF (config_inca /= 'none') THEN - ! Varaibles only needed in case of INCA - ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), tracnam(nbtr)) - END IF - + ALLOCATE(conv_flg(nbtr), pbl_flg(nbtr), tracnam(nbtr)) + conv_flg(:) = 1 ! convection activated for all tracers + pbl_flg(:) = 1 ! boundary layer activated for all tracers + !----------------------------------------------------------------------- ! 2) Choix des schemas d'advection pour l'eau et les traceurs @@ -144,5 +152,5 @@ ! Get choice of advection schema from file tracer.def or from INCA !--------------------------------------------------------------------- - IF (config_inca == 'none') THEN + IF (type_trac == 'lmdz') THEN IF(ierr.EQ.0) THEN ! Continue to read tracer.def @@ -172,5 +180,5 @@ END DO - ELSE ! config_inca='aero' ou 'chem' + ELSE ! type_trac=inca : config_inca='aero' ou 'chem' ! le module de chimie fournit les noms des traceurs ! et les schemas d'advection associes. @@ -191,5 +199,5 @@ END DO - END IF ! config_inca + END IF ! type_trac !----------------------------------------------------------------------- @@ -319,5 +327,5 @@ ! DEALLOCATE(tnom_0, hadv, vadv) - IF (config_inca /= 'none') DEALLOCATE(tracnam) + DEALLOCATE(tracnam) 999 FORMAT (i2,1x,i2,1x,a8) Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/YOECUMF.h =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/YOECUMF.h (revision 1190) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/YOECUMF.h (revision 1191) @@ -1,43 +1,48 @@ ! -! $Header$ +! $Id$ ! -C ---------------------------------------------------------------- -C* *COMMON* *YOECUMF* - PARAMETERS FOR CUMULUS MASSFLUX SCHEME -C ---------------------------------------------------------------- -C - COMMON /YOECUMF/ - L LMFPEN,LMFSCV,LMFMID,LMFDD,LMFDUDV, - R ENTRPEN,ENTRSCV,ENTRMID,ENTRDD,CMFCTOP, - R CMFCMAX,CMFCMIN,CMFDEPS,RHCDD,CPRCON -C +! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre +! veillez n'utiliser que des ! pour les commentaires +! et bien positionner les & des lignes de continuation +! (les placer en colonne 6 et en colonne 73) +! +! ---------------------------------------------------------------- +!* *COMMON* *YOECUMF* - PARAMETERS FOR CUMULUS MASSFLUX SCHEME +! ---------------------------------------------------------------- +! + COMMON /YOECUMF/ & + & LMFPEN,LMFSCV,LMFMID,LMFDD,LMFDUDV, & + & ENTRPEN,ENTRSCV,ENTRMID,ENTRDD,CMFCTOP, & + & CMFCMAX,CMFCMIN,CMFDEPS,RHCDD,CPRCON + LOGICAL LMFPEN,LMFSCV,LMFMID,LMFDD,LMFDUDV REAL ENTRPEN, ENTRSCV, ENTRMID, ENTRDD REAL CMFCTOP, CMFCMAX, CMFCMIN, CMFDEPS, RHCDD, CPRCON -c$OMP THREADPRIVATE(/YOECUMF/) -C -*if (DOC,declared) <> 'UNKNOWN' -C* *COMMON* *YOECUMF* - PARAMETERS FOR CUMULUS MASSFLUX SCHEME -C -C M.TIEDTKE E. C. M. W. F. 18/1/89 -C -C NAME TYPE PURPOSE -C ---- ---- ------- -C -C LMFPEN LOGICAL TRUE IF PENETRATIVE CONVECTION IS SWITCHED ON -C LMFSCV LOGICAL TRUE IF SHALLOW CONVECTION IS SWITCHED ON -C LMFMID LOGICAL TRUE IF MIDLEVEL CONVECTION IS SWITCHED ON -C LMFDD LOGICAL TRUE IF CUMULUS DOWNDRAFT IS SWITCHED ON -C LMFDUDV LOGICAL TRUE IF CUMULUS FRICTION IS SWITCHED ON -C ENTRPEN REAL ENTRAINMENT RATE FOR PENETRATIVE CONVECTION -C ENTRSCV REAL ENTRAINMENT RATE FOR SHALLOW CONVECTION -C ENTRMID REAL ENTRAINMENT RATE FOR MIDLEVEL CONVECTION -C ENTRDD REAL ENTRAINMENT RATE FOR CUMULUS DOWNDRAFTS -C CMFCTOP REAL RELAT. CLOUD MASSFLUX AT LEVEL ABOVE NONBUOYANC -C CMFCMAX REAL MAXIMUM MASSFLUX VALUE ALLOWED FOR -C CMFCMIN REAL MINIMUM MASSFLUX VALUE (FOR SAFETY) -C CMFDEPS REAL FRACTIONAL MASSFLUX FOR DOWNDRAFTS AT LFS -C RHCDD REAL RELATIVE SATURATION IN DOWNDRAFTS -C CPRCON REAL COEFFICIENTS FOR DETERMINING CONVERSION -C FROM CLOUD WATER TO RAIN -*ifend -C ---------------------------------------------------------------- +!$OMP THREADPRIVATE(/YOECUMF/) +! +!*if (DOC,declared) <> 'UNKNOWN' +!* *COMMON* *YOECUMF* - PARAMETERS FOR CUMULUS MASSFLUX SCHEME +! +! M.TIEDTKE E. C. M. W. F. 18/1/89 +! +! NAME TYPE PURPOSE +! ---- ---- ------- +! +! LMFPEN LOGICAL TRUE IF PENETRATIVE CONVECTION IS SWITCHED ON +! LMFSCV LOGICAL TRUE IF SHALLOW CONVECTION IS SWITCHED ON +! LMFMID LOGICAL TRUE IF MIDLEVEL CONVECTION IS SWITCHED ON +! LMFDD LOGICAL TRUE IF CUMULUS DOWNDRAFT IS SWITCHED ON +! LMFDUDV LOGICAL TRUE IF CUMULUS FRICTION IS SWITCHED ON +! ENTRPEN REAL ENTRAINMENT RATE FOR PENETRATIVE CONVECTION +! ENTRSCV REAL ENTRAINMENT RATE FOR SHALLOW CONVECTION +! ENTRMID REAL ENTRAINMENT RATE FOR MIDLEVEL CONVECTION +! ENTRDD REAL ENTRAINMENT RATE FOR CUMULUS DOWNDRAFTS +! CMFCTOP REAL RELAT. CLOUD MASSFLUX AT LEVEL ABOVE NONBUOYANC +! CMFCMAX REAL MAXIMUM MASSFLUX VALUE ALLOWED FOR +! CMFCMIN REAL MINIMUM MASSFLUX VALUE (FOR SAFETY) +! CMFDEPS REAL FRACTIONAL MASSFLUX FOR DOWNDRAFTS AT LFS +! RHCDD REAL RELATIVE SATURATION IN DOWNDRAFTS +! CPRCON REAL COEFFICIENTS FOR DETERMINING CONVERSION +! FROM CLOUD WATER TO RAIN +!*ifend +! ---------------------------------------------------------------- Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltrac.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltrac.F (revision 1190) +++ (revision ) @@ -1,122 +1,0 @@ -! -! $Header$ -! - SUBROUTINE cltrac(dtime,coef,t,tr,flux,paprs,pplay,delp, - s d_tr) - USE dimphy - IMPLICIT none -c====================================================================== -c Auteur(s): O. Boucher (LOA/LMD) date: 19961127 -c inspire de clvent -c Objet: diffusion verticale de traceurs avec flux fixe a la surface -c ou/et flux du type c-drag -c====================================================================== -c Arguments: -c dtime----input-R- intervalle du temps (en second) -c coef-----input-R- le coefficient d'echange (m**2/s) l>1 -c t--------input-R- temperature (K) -c tr-------input-R- la q. de traceurs -c flux-----input-R- le flux de traceurs a la surface -c paprs----input-R- pression a inter-couche (Pa) -c pplay----input-R- pression au milieu de couche (Pa) -c delp-----input-R- epaisseur de couche (Pa) -c cdrag----input-R- cdrag pour le flux de surface (non active) -c tr0------input-R- traceurs a la surface ou dans l'ocean (non active) -c d_tr-----output-R- le changement de tr -c flux_tr--output-R- flux de tr -c====================================================================== -cym#include "dimensions.h" -cym#include "dimphy.h" - REAL dtime - REAL coef(klon,klev) - REAL t(klon,klev), tr(klon,klev) - REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev) - REAL d_tr(klon,klev) - REAL flux(klon), cdrag(klon), tr0(klon) -c REAL flux_tr(klon,klev) -c====================================================================== -#include "YOMCST.h" -c====================================================================== - INTEGER i, k - REAL zx_ctr(klon,2:klev) - REAL zx_dtr(klon,2:klev) - REAL zx_buf(klon) - REAL zx_coef(klon,klev) - REAL local_tr(klon,klev) - REAL zx_alf1(klon), zx_alf2(klon), zx_flux(klon) -c====================================================================== - DO k = 1, klev - DO i = 1, klon - local_tr(i,k) = tr(i,k) - ENDDO - ENDDO -c - -c====================================================================== - DO i = 1, klon - zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2)) - zx_alf2(i) = 1.0 - zx_alf1(i) - zx_flux(i) = -flux(i)*dtime*RG -c--pour le moment le flux est prescrit -c--cdrag et zx_coef(1) vaut 0 - cdrag(i) = 0.0 - tr0(i) = 0.0 - zx_coef(i,1) = cdrag(i)*dtime*RG - ENDDO -c====================================================================== - 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 -c====================================================================== - DO i = 1, klon - zx_buf(i) = delp(i,1) + zx_coef(i,1)*zx_alf1(i) + zx_coef(i,2) - zx_ctr(i,2) = (local_tr(i,1)*delp(i,1)+ - . zx_coef(i,1)*tr0(i)-zx_flux(i))/zx_buf(i) - zx_dtr(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1)) / - . zx_buf(i) - ENDDO -c - DO k = 3, klev - DO i = 1, klon - zx_buf(i) = delp(i,k-1) + zx_coef(i,k) - . + zx_coef(i,k-1)*(1.-zx_dtr(i,k-1)) - zx_ctr(i,k) = (local_tr(i,k-1)*delp(i,k-1) - . +zx_coef(i,k-1)*zx_ctr(i,k-1) )/zx_buf(i) - zx_dtr(i,k) = zx_coef(i,k)/zx_buf(i) - ENDDO - ENDDO - DO i = 1, klon - local_tr(i,klev) = ( local_tr(i,klev)*delp(i,klev) - . +zx_coef(i,klev)*zx_ctr(i,klev) ) - . / ( delp(i,klev) + zx_coef(i,klev) - . -zx_coef(i,klev)*zx_dtr(i,klev) ) - ENDDO - DO k = klev-1, 1, -1 - DO i = 1, klon - local_tr(i,k) = zx_ctr(i,k+1) + zx_dtr(i,k+1)*local_tr(i,k+1) - ENDDO - ENDDO -c====================================================================== -c== flux_tr est le flux de traceur (positif vers bas) -c DO i = 1, klon -c flux_tr(i,1) = zx_coef(i,1)/(RG*dtime) -c ENDDO -c DO k = 2, klev -c DO i = 1, klon -c flux_tr(i,k) = zx_coef(i,k)/(RG*dtime) -c . * (local_tr(i,k)-local_tr(i,k-1)) -c ENDDO -c ENDDO -c====================================================================== - DO k = 1, klev - DO i = 1, klon - d_tr(i,k) = local_tr(i,k) - tr(i,k) - ENDDO - ENDDO -c - RETURN - END Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltrac.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltrac.F90 (revision 1191) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltrac.F90 (revision 1191) @@ -0,0 +1,140 @@ +! +! $Id $ +! +SUBROUTINE cltrac(dtime,coef,t,tr,flux,paprs,pplay,delp,d_tr) + USE dimphy + IMPLICIT NONE +!====================================================================== +! Auteur(s): O. Boucher (LOA/LMD) date: 19961127 +! inspire de clvent +! Objet: diffusion verticale de traceurs avec flux fixe a la surface +! ou/et flux du type c-drag +! +! Arguments: +!----------- +! dtime....input-R- intervalle du temps (en secondes) +! coef.....input-R- le coefficient d'echange (m**2/s) l>1 +! t........input-R- temperature (K) +! tr.......input-R- la q. de traceurs +! flux.....input-R- le flux de traceurs a la surface +! paprs....input-R- pression a inter-couche (Pa) +! pplay....input-R- pression au milieu de couche (Pa) +! delp.....input-R- epaisseur de couche (Pa) +! cdrag....input-R- cdrag pour le flux de surface (non active) +! tr0......input-R- traceurs a la surface ou dans l'ocean (non active) +! d_tr.....output-R- le changement de tr +! flux_tr..output-R- flux de tr +!====================================================================== + include "YOMCST.h" +! +! Entree +! + REAL,INTENT(IN) :: dtime + REAL,DIMENSION(klon,klev),INTENT(IN) :: coef + REAL,DIMENSION(klon,klev),INTENT(IN) :: t, tr + REAL,DIMENSION(klon),INTENT(IN) :: flux !(at/s/m2) + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs + REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay, delp +! +! Sorties +! + REAL ,DIMENSION(klon,klev),INTENT(OUT) :: d_tr +! REAL ,DIMENSION(klon,klev),INTENT(OUT) :: flux_tr +! +! Local +! + INTEGER :: i, k + REAL,DIMENSION(klon) :: cdrag, tr0 + REAL,DIMENSION(klon,klev) :: zx_ctr + REAL,DIMENSION(klon,klev) :: zx_dtr + REAL,DIMENSION(klon) :: zx_buf + REAL,DIMENSION(klon,klev) :: zx_coef + REAL,DIMENSION(klon,klev) :: local_tr + REAL,DIMENSION(klon) :: zx_alf1,zx_alf2,zx_flux + +!====================================================================== + + DO k = 1, klev + DO i = 1, klon + local_tr(i,k) = tr(i,k) + ENDDO + ENDDO + +!====================================================================== + + DO i = 1, klon + zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2)) + zx_alf2(i) = 1.0 - zx_alf1(i) + zx_flux(i) = -flux(i)*dtime*RG +! Pour le moment le flux est prescrit cdrag et zx_coef(1) vaut 0 + cdrag(i) = 0.0 + tr0(i) = 0.0 + zx_coef(i,1) = cdrag(i)*dtime*RG + zx_ctr(i,1)=0. + zx_dtr(i,1)=0. + 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,1) + zx_coef(i,1)*zx_alf1(i) + zx_coef(i,2) + ! + zx_ctr(i,2) = (local_tr(i,1)*delp(i,1)+ & + zx_coef(i,1)*tr0(i)-zx_flux(i))/zx_buf(i) + ! + zx_dtr(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1)) / & + zx_buf(i) + ENDDO + + DO k = 3, klev + DO i = 1, klon + zx_buf(i) = delp(i,k-1) + zx_coef(i,k) & + + zx_coef(i,k-1)*(1.-zx_dtr(i,k-1)) + zx_ctr(i,k) = (local_tr(i,k-1)*delp(i,k-1) & + +zx_coef(i,k-1)*zx_ctr(i,k-1) )/zx_buf(i) + zx_dtr(i,k) = zx_coef(i,k)/zx_buf(i) + ENDDO + ENDDO + + DO i = 1, klon + local_tr(i,klev) = ( local_tr(i,klev)*delp(i,klev) & + +zx_coef(i,klev)*zx_ctr(i,klev) ) & + / ( delp(i,klev) + zx_coef(i,klev) & + -zx_coef(i,klev)*zx_dtr(i,klev) ) + ENDDO + + DO k = klev-1, 1, -1 + DO i = 1, klon + local_tr(i,k) = zx_ctr(i,k+1) + zx_dtr(i,k+1)*local_tr(i,k+1) + ENDDO + ENDDO + +!====================================================================== +!== flux_tr est le flux de traceur (positif vers bas) +! DO i = 1, klon +! flux_tr(i,1) = zx_coef(i,1)/(RG*dtime) +! ENDDO +! DO k = 2, klev +! DO i = 1, klon +! flux_tr(i,k) = zx_coef(i,k)/(RG*dtime) +! . * (local_tr(i,k)-local_tr(i,k-1)) +! ENDDO +! ENDDO +!====================================================================== + DO k = 1, klev + DO i = 1, klon + d_tr(i,k) = local_tr(i,k) - tr(i,k) + ENDDO + ENDDO + +END SUBROUTINE cltrac Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltracrn.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltracrn.F (revision 1190) +++ (revision ) @@ -1,294 +1,0 @@ -! - SUBROUTINE cltracrn( itr, dtime,u1lay, v1lay, - e cdrag,coef,t,ftsol,pctsrf, - e tr,trs,paprs,pplay,delp, - e masktr,fshtr,hsoltr,tautr,vdeptr, - e lat, - s d_tr,d_trs ) - - USE dimphy - IMPLICIT none -c====================================================================== -c Auteur(s): Alex/LMD) date: fev 99 -c inspire de clqh + clvent -c Objet: diffusion verticale de traceurs avec quantite de traceur ds -c le sol ( reservoir de sol de radon ) -c -c note : pour l'instant le traceur dans le sol et le flux sont -c calcules mais ils ne servent que de diagnostiques -c seule la tendance sur le traceur est sortie (d_tr) -c====================================================================== -c Arguments: -c itr--- -input-R- le type de traceur 1- Rn 2 - Pb -c dtime----input-R- intervalle du temps (en second) -c u1lay----input-R- vent u de la premiere couche (m/s) -c v1lay----input-R- vent v de la premiere couche (m/s) -c cdrag----input-R- cdrag -c coef-----input-R- le coefficient d'echange (m**2/s) l>1, valable uniquement pour k entre 2 et klev -c t--------input-R- temperature (K) -c paprs----input-R- pression a inter-couche (Pa) -c pplay----input-R- pression au milieu de couche (Pa) -c delp-----input-R- epaisseur de couche (Pa) -c ftsol----input-R- temperature du sol (en Kelvin) -c tr-------input-R- traceurs -c trs------input-R- traceurs dans le sol -c masktr---input-R- Masque reservoir de sol traceur (1 = reservoir) -c fshtr----input-R- Flux surfacique de production dans le sol -c tautr----input-R- Constante de decroissance du traceur -c vdeptr---input-R- Vitesse de depot sec dans la couche brownienne -c hsoltr---input-R- Epaisseur equivalente du reservoir de sol -c lat-----input-R- latitude en degree -c d_tr-----output-R- le changement de "tr" -c d_trs----output-R- le changement de "trs" -c====================================================================== -cym#include "dimensions.h" -cym#include "dimphy.h" -#include "YOMCST.h" -#include "indicesol.h" -c====================================================================== - REAL dtime - REAL u1lay(klon), v1lay(klon) - REAL cdrag(klon) - REAL coef(klon,klev) - REAL t(klon,klev), ftsol(klon,nbsrf), pctsrf(klon,nbsrf) - REAL tr(klon,klev), trs(klon) - REAL paprs(klon,klev+1), pplay(klon,klev), delp(klon,klev) - REAL masktr(klon) - REAL fshtr(klon) - REAL hsoltr - REAL tautr - REAL vdeptr - REAL lat(klon) - REAL d_tr(klon,klev) -c====================================================================== - REAL flux_tr(klon,klev) ! (diagnostic) flux de traceur - REAL d_trs(klon) ! (diagnostic) traceur ds le sol -c====================================================================== - INTEGER i, k, itr, n, l - REAL rotrhi(klon) - REAL zx_coef(klon,klev) - REAL zx_buf(klon) - REAL zx_ctr(klon,klev) - REAL zx_dtr(klon,klev) - REAL zx_trs(klon) - REAL zx_a, zx_b - - REAL local_tr(klon,klev) - REAL local_trs(klon) - REAL zts(klon) - REAL zx_alpha1(klon), zx_alpha2(klon) -c====================================================================== -cAA Pour l'instant les 4 types de surface ne sont pas pris en compte -cAA On fabrique avec zts un champ de temperature de sol -cAA que le pondere par la fraction de nature de sol. -c -c print*,'PASSAGE DANS CLTRACRN' - - DO i = 1,klon - zts(i) = 0. - ENDDO -c - DO n=1,nbsrf - DO i = 1,klon - zts(i) = zts(i) + ftsol(i,n)*pctsrf(i,n) - ENDDO - ENDDO -c - DO i = 1,klon - rotrhi(i) = RD * zts(i) / hsoltr - END DO -c - DO k = 1, klev - DO i = 1, klon - local_tr(i,k) = tr(i,k) - ENDDO - ENDDO -c - DO i = 1, klon - local_trs(i) = trs(i) - ENDDO -c====================================================================== -cAA Attention si dans clmain zx_alf1(i) = 1.0 -cAA Il doit y avoir coherence (dc la meme chose ici) - - DO i = 1, klon -cAA 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 -c====================================================================== - 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 -c - 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 -c====================================================================== - 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 -c - 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 -c - 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 -c====================================================================== -c Calculer d'abord local_trs nouvelle quantite dans le reservoir -c de sol -c -c------------------------- -c Au dessus des continents -c------------------------- -c Le pb peut se deposer partout : vdeptr = 10-3 m/s -c Le Rn est traiter commme une couche Brownienne puisque vdeptr = 0. -c - DO i = 1, klon -c - 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)) - 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 -cAA: Pour l'instant, pour aller vite, le depot sec est traite -C comme une decroissance - . + dtime * vdeptr / hsoltr - zx_trs(i) = zx_a / zx_b - local_trs(i) = zx_trs(i) - ENDIF -c -c Si on est entre 60N et 70N on divise par 2 l'emanation -c-------------------------------------------------------- -c - 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 - 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 -c -c---------------------------------------------- -c Au dessus des oceans et aux hautes latitudes -c---------------------------------------------- -c -c au dessous de -60S pas d'emission de radon au dessus -c des oceans et des continents -c--------------------------------------------------------------- - - 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 - zx_trs(i) = 0. - local_trs(i) = 0. - END IF - -c au dessus de 70 N pas d'emission de radon au dessus -c des oceans et des continents -c-------------------------------------------------------------- - 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 - zx_trs(i) = 0. - local_trs(i) = 0. - END IF - -c Au dessus des oceans la source est nulle -c----------------------------------------- -c - IF (itr.eq.1.AND.NINT(masktr(i)).EQ.0) THEN - zx_trs(i) = 0. - local_trs(i) = 0. - END IF -c - ENDDO ! sur le i=1,klon -c -c====================================================================== -c==== une fois on a zx_trs, on peut faire l'iteration ======== -c - 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 -c====================================================================== -c== Calcul du flux de traceur (flux_tr): UA/(m**2 s) -c - 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 -c====================================================================== -c== Calcul des tendances du traceur ds le sol et dans l'atmosphere -c - 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 -c====================================================================== -c - RETURN - END Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltracrn.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltracrn.F90 (revision 1191) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/cltracrn.F90 (revision 1191) @@ -0,0 +1,286 @@ +!$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 + 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 1- Rn 2 - Pb +! 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" +!====================================================================== + include "YOMCST.h" + include "indicesol.h" +! +!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.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 + 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.1.AND.NINT(masktr(i)).EQ.0).OR. & + (itr.EQ.1.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.1.AND.NINT(masktr(i)).EQ.0).OR. & + (itr.EQ.1.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.1.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 Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/cvltr.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/cvltr.F (revision 1190) +++ (revision ) @@ -1,145 +1,0 @@ -c -c $Header$ -c - SUBROUTINE cvltr(pdtime,da, phi, mp,paprs,pplay,x,upd,dnd,dx) - USE dimphy - IMPLICIT NONE -c===================================================================== -c Objet : convection des traceurs / KE -c Auteurs: M-A Filiberti and J-Y Grandpeix -c===================================================================== -c -cym#include "dimensions.h" -cym#include "dimphy.h" -#include "YOMCST.h" -#include "YOECUMF.h" -c - REAL pdtime - REAL paprs(klon,klev+1) ! pression aux 1/2 couches (bas en haut) - REAL pplay(klon,klev) ! pression pour le milieu de chaque couche - REAL x(klon,klev) ! q de traceur (bas en haut) - REAL dx(klon,klev) ! tendance de traceur (bas en haut) - real da(klon,klev),phi(klon,klev,klev),mp(klon,klev) - REAL upd(klon,klev) ! saturated updraft mass flux - REAL dnd(klon,klev) ! saturated downdraft mass flux -c -c--variables locales - real zed(klon,klev),zmd(klon,klev,klev) - real za(klon,klev,klev) - real zmfd(klon,klev),zmfa(klon,klev) - real zmfp(klon,klev),zmfu(klon,klev) - integer i,k,j -c test conservation -c real conserv -c ========================================= -c calcul des tendances liees au downdraft -c ========================================= - zed(:,:)=0. - zmfd(:,:)=0. - zmfa(:,:)=0. - zmfu(:,:)=0. - zmfp(:,:)=0. - zmd(:,:,:)=0. - za(:,:,:)=0. -c entrainement - do k=1,klev-1 - do i=1,klon - zed(i,k)=max(0.,mp(i,k)-mp(i,k+1)) - end do - end do -c -c calcul de la matrice d echange -c matrice de distribution de la masse entrainee en k -c - do k=1,klev - do i=1,klon - zmd(i,k,k)=zed(i,k) - end do - end do - do k=2,klev - do j=k-1,1,-1 - do i=1,klon - if(mp(i,j+1).ne.0) then - zmd(i,j,k)=zmd(i,j+1,k)*min(1.,mp(i,j)/mp(i,j+1)) - endif - end do - end do - end do - do k=1,klev - do j=1,klev-1 - do i=1,klon - za(i,j,k)=max(0.,zmd(i,j+1,k)-zmd(i,j,k)) - end do - end do - end do -c -c rajout du terme lie a l ascendance induite -c - do j=2,klev - do i=1,klon - za(i,j,j-1)=za(i,j,j-1)+mp(i,j) - end do - end do -C -c tendances -c - do k=1,klev - do j=1,klev - do i=1,klon - zmfd(i,j)=zmfd(i,j)+za(i,j,k)*(x(i,k)-x(i,j)) - end do - end do - end do -c -c ========================================= -c calcul des tendances liees aux flux satures -c ========================================= - do j=1,klev - do i=1,klon - zmfa(i,j)=da(i,j)*(x(i,1)-x(i,j)) - end do - end do - do k=1,klev - do j=1,klev - do i=1,klon - zmfp(i,j)=zmfp(i,j)+phi(i,j,k)*(x(i,k)-x(i,j)) - end do - end do - end do - do j=1,klev-1 - do i=1,klon - zmfu(i,j)=max(0.,upd(i,j+1)+dnd(i,j+1))*(x(i,j+1)-x(i,j)) - end do - end do - do j=2,klev - do i=1,klon - zmfu(i,j)=zmfu(i,j) - . +min(0.,upd(i,j)+dnd(i,j))*(x(i,j)-x(i,j-1)) - end do - end do - -c ========================================= -c--calcul final des tendances -c ========================================= - do k=1, klev - do i=1, klon - dx(i,k)=(zmfd(i,k)+zmfu(i,k) - . +zmfa(i,k)+zmfp(i,k))*pdtime - . *RG/(paprs(i,k)-paprs(i,k+1)) -c print*,'dx',k,dx(i,k) - enddo - enddo -c -c test de conservation du traceur -c conserv=0. -c do k=1, klev -c do i=1, klon -c conserv=conserv+dx(i,k)* -c . (paprs(i,k)-paprs(i,k+1))/RG -C -c enddo -c enddo -c print *,'conserv',conserv - - return - end Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/cvltr.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/cvltr.F90 (revision 1191) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/cvltr.F90 (revision 1191) @@ -0,0 +1,145 @@ +! +! $Id $ +! +SUBROUTINE cvltr(pdtime,da, phi, mp,paprs,pplay,x,upd,dnd,dx) + USE dimphy + IMPLICIT NONE +!===================================================================== +! Objet : convection des traceurs / KE +! Auteurs: M-A Filiberti and J-Y Grandpeix +!===================================================================== + + include "YOMCST.h" + include "YOECUMF.h" + +! Entree + REAL,INTENT(IN) :: pdtime + REAL,DIMENSION(klon,klev),INTENT(IN) :: da + REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi + REAL,DIMENSION(klon,klev),INTENT(IN) :: mp + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression aux 1/2 couches (bas en haut) + REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le milieu de chaque couche + REAL,DIMENSION(klon,klev),INTENT(IN) :: x ! q de traceur (bas en haut) + REAL,DIMENSION(klon,klev),INTENT(IN) :: upd ! saturated updraft mass flux + REAL,DIMENSION(klon,klev),INTENT(IN) :: dnd ! saturated downdraft mass flux + +! Sortie + REAL,DIMENSION(klon,klev),INTENT(OUT) :: dx ! tendance de traceur (bas en haut) + +! Variables locales + REAL,DIMENSION(klon,klev) :: zed + REAL,DIMENSION(klon,klev,klev) :: zmd + REAL,DIMENSION(klon,klev,klev) :: za + REAL,DIMENSION(klon,klev) :: zmfd,zmfa + REAL,DIMENSION(klon,klev) :: zmfp,zmfu + INTEGER :: i,k,j + +! ========================================= +! calcul des tendances liees au downdraft +! ========================================= + zed(:,:)=0. + zmfd(:,:)=0. + zmfa(:,:)=0. + zmfu(:,:)=0. + zmfp(:,:)=0. + zmd(:,:,:)=0. + za(:,:,:)=0. +! entrainement + DO k=1,klev-1 + DO i=1,klon + zed(i,k)=max(0.,mp(i,k)-mp(i,k+1)) + END DO + END DO + +! calcul de la matrice d echange +! matrice de distribution de la masse entrainee en k + + DO k=1,klev + DO i=1,klon + zmd(i,k,k)=zed(i,k) + END DO + END DO + DO k=2,klev + DO j=k-1,1,-1 + DO i=1,klon + if(mp(i,j+1).ne.0) then + zmd(i,j,k)=zmd(i,j+1,k)*min(1.,mp(i,j)/mp(i,j+1)) + ENDif + END DO + END DO + END DO + DO k=1,klev + DO j=1,klev-1 + DO i=1,klon + za(i,j,k)=max(0.,zmd(i,j+1,k)-zmd(i,j,k)) + END DO + END DO + END DO +! +! rajout du terme lie a l ascendance induite +! + DO j=2,klev + DO i=1,klon + za(i,j,j-1)=za(i,j,j-1)+mp(i,j) + END DO + END DO +! +! tendances +! + DO k=1,klev + DO j=1,klev + DO i=1,klon + zmfd(i,j)=zmfd(i,j)+za(i,j,k)*(x(i,k)-x(i,j)) + END DO + END DO + END DO +! +! ========================================= +! calcul des tendances liees aux flux satures +! ========================================= + DO j=1,klev + DO i=1,klon + zmfa(i,j)=da(i,j)*(x(i,1)-x(i,j)) + END DO + END DO + DO k=1,klev + DO j=1,klev + DO i=1,klon + zmfp(i,j)=zmfp(i,j)+phi(i,j,k)*(x(i,k)-x(i,j)) + END DO + END DO + END DO + DO j=1,klev-1 + DO i=1,klon + zmfu(i,j)=max(0.,upd(i,j+1)+dnd(i,j+1))*(x(i,j+1)-x(i,j)) + END DO + END DO + DO j=2,klev + DO i=1,klon + zmfu(i,j)=zmfu(i,j)+min(0.,upd(i,j)+dnd(i,j))*(x(i,j)-x(i,j-1)) + END DO + END DO + +! ========================================= +! calcul final des tendances +! ========================================= + DO k=1, klev + DO i=1, klon + dx(i,k)=(zmfd(i,k)+zmfu(i,k) & + +zmfa(i,k)+zmfp(i,k))*pdtime & + *RG/(paprs(i,k)-paprs(i,k+1)) + ! print*,'dx',k,dx(i,k) + ENDDO + ENDDO + +! test de conservation du traceur +! conserv=0. +! DO k=1, klev +! DO i=1, klon +! conserv=conserv+dx(i,k)* & +! (paprs(i,k)-paprs(i,k+1))/RG +! ENDDO +! ENDDO +! print *,'conserv',conserv + +END SUBROUTINE cvltr Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/ini_histrac.h =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/ini_histrac.h (revision 1190) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/ini_histrac.h (revision 1191) @@ -1,123 +1,128 @@ ! -! $Header$ +! $Id $ ! - IF (ecrit_tra>0. .AND. config_inca == 'none') THEN -c$OMP MASTER - CALL ymds2ju(annee_ref, 1, day_ref, 0.0, zjulian) -c - CALL histbeg_phy("histrac", itau_phy, zjulian, pdtphys, - . nhori, nid_tra) - CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb", - . klev, presnivs, nvert) +#ifdef CPP_IOIPSL + IF (ecrit_tra>0. .AND. config_inca == 'none') THEN +!$OMP MASTER + CALL ymds2ju(annee_ref, 1, day_ref, 0.0, zjulian) + CALL histbeg_phy("histrac", itau_phy, zjulian, pdtphys,nhori, nid_tra) + CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb",klev, presnivs, nvert) + zsto = pdtphys + zout = ecrit_tra + CALL histdef(nid_tra, "phis", "Surface geop. height", "-", & + iim,jj_nb,nhori, 1,1,1, -99, 32,"once", zsto,zout) + CALL histdef(nid_tra, "aire", "Grid area", "-", & + iim,jj_nb,nhori, 1,1,1, -99, 32,"once", zsto,zout) - zsto = pdtphys - zout = ecrit_tra -c - CALL histdef(nid_tra, "phis", "Surface geop. height", "-", - . iim,jj_nb,nhori, 1,1,1, -99, 32, - . "once", zsto,zout) -c - CALL histdef(nid_tra, "aire", "Grid area", "-", - . iim,jj_nb,nhori, 1,1,1, -99, 32, - . "once", zsto,zout) - DO it=1,nbtr -C champ 2D - iq=it+2 - iiq=niadv(iq) - CALL histdef(nid_tra, tname(iiq), ttext(iiq), "U/kga", - . iim,jj_nb,nhori, klev,1,klev,nvert, 32, - . "ave(X)", zsto,zout) - if (lessivage) THEN - CALL histdef(nid_tra, "fl"//tname(iiq),"Flux "//ttext(iiq), - . "U/m2/s",iim,jj_nb,nhori, klev,1,klev,nvert, 32, - . "ave(X)", zsto,zout) - endif +!TRACEURS +!---------------- + DO it = 1,nbtr + iiq = niadv(it+2) -c---Ajout Olivia - CALL histdef(nid_tra, "d_tr_th_"//tname(iiq), - . "tendance thermique"// ttext(iiq), "?", - . iim,jj_nb,nhori, klev,1,klev,nvert, 32, - . "ave(X)", zsto,zout) -c - if(iflag_con.GE.2) then - CALL histdef(nid_tra, "d_tr_cv_"//tname(iiq), - . "tendance convection"// ttext(iiq), "?", - . iim,jj_nb,nhori, klev,1,klev,nvert, 32, - . "ave(X)", zsto,zout) - endif !(iflag_con.GE.2) then - CALL histdef(nid_tra, "d_tr_cl_"//tname(iiq), - . "tendance couche limite"// ttext(iiq), "?", - . iim,jj_nb,nhori, klev,1,klev,nvert, 32, - . "ave(X)", zsto,zout) -c---fin Olivia +! CONCENTRATIONS + CALL histdef(nid_tra, tname(iiq), ttext(iiq), "U/kga", & + iim,jj_nb,nhori, klev,1,klev,nvert, 32,"ave(X)", zsto,zout) - ENDDO +! TD LESSIVAGE + IF (lessivage .AND. aerosol(it)) THEN + CALL histdef(nid_tra, "fl"//tname(iiq),"Flux "//ttext(iiq), & + "U/m2/s",iim,jj_nb,nhori, klev,1,klev,nvert, 32, & + "ave(X)", zsto,zout) + END IF - CALL histdef(nid_tra, "pyu1", "Vent niv 1", "-", - . iim,jj_nb,nhori, 1,1,1, -99, 32, - . "inst(X)", zout,zout) +! TD THERMIQUES + IF (iflag_thermals.gt.0) THEN + CALL histdef(nid_tra, "d_tr_th_"//tname(iiq), & + "tendance thermique"// ttext(iiq), "?", & + iim,jj_nb,nhori, klev,1,klev,nvert, 32, & + "ave(X)", zsto,zout) + ENDIF - CALL histdef(nid_tra, "pyv1", "Vent niv 1", "-", - . iim,jj_nb,nhori, 1,1,1, -99, 32, - . "inst(X)", zout,zout) - CALL histdef(nid_tra, "psrf1", "nature sol", "-", - . iim,jj_nb,nhori, 1,1,1, -99, 32, - . "inst(X)", zout,zout) - CALL histdef(nid_tra, "psrf2", "nature sol", "-", - . iim,jj_nb,nhori, 1,1,1, -99, 32, - . "inst(X)", zout,zout) - CALL histdef(nid_tra, "psrf3", "nature sol", "-", - . iim,jj_nb,nhori, 1,1,1, -99, 32, - . "inst(X)", zout,zout) - CALL histdef(nid_tra, "psrf4", "nature sol", "-", - . iim,jj_nb,nhori, 1,1,1, -99, 32, - . "inst(X)", zout,zout) - CALL histdef(nid_tra, "ftsol1", "temper sol", "-", - . iim,jj_nb,nhori, 1,1,1, -99, 32, - . "inst(X)", zout,zout) - CALL histdef(nid_tra, "ftsol2", "temper sol", "-", - . iim,jj_nb,nhori, 1,1,1, -99, 32, - . "inst(X)", zout,zout) - CALL histdef(nid_tra, "ftsol3", "temper sol", "-", - . iim,jj_nb,nhori, 1,1,1, -99, 32, - . "inst", zout,zout) - CALL histdef(nid_tra, "ftsol4", "temper sol", "-", - . iim,jj_nb,nhori, 1,1,1, -99, 32, - . "inst(X)", zout,zout) - CALL histdef(nid_tra, "pplay", "flux u mont","-", - . iim,jj_nb,nhori, klev,1,klev,nvert, 32, - . "inst(X)", zout,zout) - CALL histdef(nid_tra, "t", "flux u mont","-", - . iim,jj_nb,nhori, klev,1,klev,nvert, 32, - . "inst(X)", zout,zout) - CALL histdef(nid_tra, "mfu", "flux u mont","-", - . iim,jj_nb,nhori, klev,1,klev,nvert, 32, - . "ave(X)", zsto,zout) - CALL histdef(nid_tra, "mfd", "flux u decen","-", - . iim,jj_nb,nhori, klev,1,klev,nvert, 32, - . "ave(X)", zsto,zout) - CALL histdef(nid_tra, "en_u", "flux u mont","-", - . iim,jj_nb,nhori, klev,1,klev,nvert, 32, - . "ave(X)", zsto,zout) - CALL histdef(nid_tra, "en_d", "flux u mont","-", - . iim,jj_nb,nhori, klev,1,klev,nvert, 32, - . "ave(X)", zsto,zout) - CALL histdef(nid_tra, "de_d", "flux u mont","-", - . iim,jj_nb,nhori, klev,1,klev,nvert, 32, - . "ave(X)", zsto,zout) - CALL histdef(nid_tra, "de_u", "flux u decen","-", - . iim,jj_nb,nhori, klev,1,klev,nvert, 32, - . "ave(X)", zsto,zout) - CALL histdef(nid_tra, "coefh", "turbulent coef","-", - . iim,jj_nb,nhori, klev,1,klev,nvert, 32, - . "ave(X)", zsto,zout) +! TD CONVECTION + IF (iflag_con.GE.2) THEN + CALL histdef(nid_tra, "d_tr_cv_"//tname(iiq), & + "tendance convection"// ttext(iiq), "?", & + iim,jj_nb,nhori, klev,1,klev,nvert, 32, & + "ave(X)", zsto,zout) + ENDIF -c - CALL histend(nid_tra) - ndex2d = 0 - ndex3d = 0 - ndex = 0 -c$OMP END MASTER - END IF ! ecrit_tra>0. .AND. config_inca == 'none' +! TD COUCHE-LIMITE + CALL histdef(nid_tra, "d_tr_cl_"//tname(iiq), & + "tendance couche limite"// ttext(iiq), "?", & + iim,jj_nb,nhori, klev,1,klev,nvert, 32, & + "ave(X)", zsto,zout) + ENDDO +!--------------- +! +! VENT (niveau 1) + CALL histdef(nid_tra, "pyu1", "Vent niv 1", "-", & + iim,jj_nb,nhori, 1,1,1, -99, 32, & + "inst(X)", zout,zout) + CALL histdef(nid_tra, "pyv1", "Vent niv 1", "-", & + iim,jj_nb,nhori, 1,1,1, -99, 32, & + "inst(X)", zout,zout) + +! TEMPERATURE DU SOL + CALL histdef(nid_tra, "ftsol1", "temper sol", "-", & + iim,jj_nb,nhori, 1,1,1, -99, 32, & + "inst(X)", zout,zout) + CALL histdef(nid_tra, "ftsol2", "temper sol", "-", & + iim,jj_nb,nhori, 1,1,1, -99, 32, & + "inst(X)", zout,zout) + CALL histdef(nid_tra, "ftsol3", "temper sol", "-", & + iim,jj_nb,nhori, 1,1,1, -99, 32, & + "inst", zout,zout) + CALL histdef(nid_tra, "ftsol4", "temper sol", "-", & + iim,jj_nb,nhori, 1,1,1, -99, 32, & + "inst(X)", zout,zout) + +! NATURE DU SOL + CALL histdef(nid_tra, "psrf1", "nature sol", "-", & + iim,jj_nb,nhori, 1,1,1, -99, 32, & + "inst(X)", zout,zout) + CALL histdef(nid_tra, "psrf2", "nature sol", "-", & + iim,jj_nb,nhori, 1,1,1, -99, 32, & + "inst(X)", zout,zout) + CALL histdef(nid_tra, "psrf3", "nature sol", "-", & + iim,jj_nb,nhori, 1,1,1, -99, 32, & + "inst(X)", zout,zout) + CALL histdef(nid_tra, "psrf4", "nature sol", "-", & + iim,jj_nb,nhori, 1,1,1, -99, 32, & + "inst(X)", zout,zout) +! DIVERS + CALL histdef(nid_tra, "pplay", "flux u mont","-", & + iim,jj_nb,nhori, klev,1,klev,nvert, 32, & + "inst(X)", zout,zout) + CALL histdef(nid_tra, "t", "flux u mont","-", & + iim,jj_nb,nhori, klev,1,klev,nvert, 32, & + "inst(X)", zout,zout) + CALL histdef(nid_tra, "mfu", "flux u mont","-", & + iim,jj_nb,nhori, klev,1,klev,nvert, 32, & + "ave(X)", zsto,zout) + CALL histdef(nid_tra, "mfd", "flux u decen","-", & + iim,jj_nb,nhori, klev,1,klev,nvert, 32, & + "ave(X)", zsto,zout) + CALL histdef(nid_tra, "en_u", "flux u mont","-", & + iim,jj_nb,nhori, klev,1,klev,nvert, 32, & + "ave(X)", zsto,zout) + CALL histdef(nid_tra, "en_d", "flux u mont","-", & + iim,jj_nb,nhori, klev,1,klev,nvert, 32, & + "ave(X)", zsto,zout) + CALL histdef(nid_tra, "de_d", "flux u mont","-", & + iim,jj_nb,nhori, klev,1,klev,nvert, 32, & + "ave(X)", zsto,zout) + CALL histdef(nid_tra, "de_u", "flux u decen","-", & + iim,jj_nb,nhori, klev,1,klev,nvert, 32, & + "ave(X)", zsto,zout) + CALL histdef(nid_tra, "coefh", "turbulent coef","-", & + iim,jj_nb,nhori, klev,1,klev,nvert, 32, & + "ave(X)", zsto,zout) + + CALL histend(nid_tra) +!$OMP END MASTER + END IF ! ecrit_tra>0. .AND. config_inca == 'none' + +#endif + Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/init_be.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/init_be.F90 (revision 1191) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/init_be.F90 (revision 1191) @@ -0,0 +1,510 @@ +!$Id $ + +SUBROUTINE init_be(pctsrf,masktr,tautr,vdeptr,scavtr,srcbe) + + USE dimphy + USE comgeomphy + USE infotrac, ONLY : nbtr + + IMPLICIT NONE +!===================================================================== +! Objet : prescription d'une source de Beryllium 7 +! pour 19 niveaux verticaux +! (d'apres le diagramme de Lal and Peters, 1967) +! +! +! written by : O. Coindreau (CEA/LDG) 05/2005 +! last modified by : A. Jamelot (LMD/CEA) 04/03/2009 +!===================================================================== + + INCLUDE "YOMCST.h" + INCLUDE "YOECUMF.h" + INCLUDE "indicesol.h" + +! +! Input Arguments +! + REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf !Pourcentage de sol (f(nature du sol)) +! +! Output Arguments +! + REAL,DIMENSION(klon),INTENT(OUT) :: masktr ! Masque de l'echange avec la surface (possible => 1 ) + REAL,INTENT(OUT) :: tautr ! Constante de decroissance radioactive + REAL,INTENT(OUT) :: vdeptr ! Vitesse de depot sec dans la couche Brownienne + REAL,INTENT(OUT) :: scavtr ! Coefficient de lessivage + REAL,DIMENSION(klon,klev),INTENT(OUT) :: srcbe ! source volumique de 7Be +! +! Local Variables +! + REAL,DIMENSION(klon) :: rlatgeo ! latitudes geomagnetiques de la grille + REAL :: glt ! latitude du pole geomagnetique + REAL :: glg ! longitude du pole geomagnetique + REAL :: latgeo,qcos + INTEGER :: k,i + + WRITE(*,*)'PASSAGE init_be ...' + +! Source actuellement definie pour klev = 19 et klev >= 39 + IF (klev /= 19 .AND. klev<39) CALL abort_gcm("init_be","Source du be7 necessite klev=19 ou klev>=39",1) +! +! Definition des constantes +! ------------------------- + tautr = 6645000. + vdeptr = 1.E-3 + scavtr = 0.5 + + WRITE(*,*) '-------------- SOURCE DE BERYLLIUM ------------------- ' + WRITE(*,*)'Decroissance (s): ', tautr + WRITE(*,*)'Vitesse de depot sec: ',vdeptr + WRITE(*,*)'Facteur de lessivage: ',scavtr + + DO i = 1,klon + masktr(i) = 0. + IF ( NINT(pctsrf(i,1)) .EQ. 1 ) masktr(i) = 1. + END DO + +! Premiers niveaux: source nulle +! ------------------------------ + DO k = 1,6 + DO i = 1,klon + srcbe(i,k) = 0. + END DO + END DO +! +! Pour les autres niveaux: +! 1-passer des coordonnees geographiques a la latitude geomagnetique +! 2-prescrire la source de Be (en 10exp5 at/g/s) dans ce repere +! 3-mettre la source de Be ds la bonne unite (en at/kgA/s) +! + glt=78.5*rpi/180. + glg=69.0*rpi/180. + + DO i = 1,klon + qcos=sin(glt)*sin(rlatd(i)) + qcos=qcos+cos(glt)*cos(rlatd(i))*cos(rlond(i)+glg) + IF ( qcos .LT. -1.) qcos = -1. + IF ( qcos .GT. 1.) qcos = 1. + rlatgeo(i)=rpi/2.-acos(qcos) + ENDDO + +!=========================== +! Cas 19 niveaux verticaux +!=========================== + IF (klev.eq.19) then + DO k = 1,klev + DO i = 1,klon + latgeo=(180./rpi)*abs(rlatgeo(i)) + IF ( k .EQ. 1 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=0.1 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.09 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.07 + END IF + IF ( k .EQ. 2 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=0.12 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.1 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.09 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.07 + END IF + IF ( k .EQ. 3 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=0.14 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.12 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.1 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.09 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.07 + END IF + IF ( k .EQ. 4 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=0.175 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.16 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.14 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.12 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.1 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.09 + END IF + IF ( k .EQ. 5 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=0.28 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.26 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.23 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.175 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.14 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.12 + END IF + IF ( k .EQ. 6 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=0.56 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.49 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.42 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.28 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.26 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.245 + END IF + IF ( k .EQ. 7 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=1.05 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.875 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.7 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.52 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.44 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.385 + END IF + IF ( k .EQ. 8 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=2. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=1.8 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=1.5 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=1. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.8 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.75 + END IF + IF ( k .EQ. 9 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=4. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=3.5 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=3. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=2.5 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=1.8 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=1.4 + END IF + IF ( k .EQ. 10 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=8.5 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=8. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=7. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=4.5 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=3.5 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=3. + END IF + IF ( k .EQ. 11 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=17. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=15. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=11. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=8. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=5. + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=4. + END IF + IF ( k .EQ. 12 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=25. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=22. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=17. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=11. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=7.5 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7. + END IF + IF ( k .EQ. 13 ) THEN + IF (latgeo.GE.60.0) srcbe(i,k)=33. + IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=32. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=30. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=22. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=15. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=11. + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=8. + END IF + IF ( k .EQ. 14 ) THEN + IF (latgeo.GE.60.0) srcbe(i,k)=48. + IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=45. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=36. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=26. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=17.5 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=12.5 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10. + END IF + IF ( k .EQ. 15 ) THEN + IF (latgeo.GE.70.0) srcbe(i,k)=58. + IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=57. + IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=50. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=38. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=25. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=15. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=12.5 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10. + END IF + IF ( k .EQ. 16 ) THEN + IF (latgeo.GE.70.0) srcbe(i,k)=70. + IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=65. + IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=50. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=32. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=20. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=13. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=9. + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.5 + END IF + IF ( k .GE. 17 ) THEN + IF (latgeo.GE.70.0) srcbe(i,k)=80. + IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=70. + IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=45. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=27. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=17.5 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=12. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=8. + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7. + END IF + END DO + END DO + END IF ! fin de 19 niveaux verticaux + +!================================ +! Cas 39 niveaux verticaux +!================================ + IF (klev .ge. 39) then + DO k = 1,klev + DO i = 1,klon + latgeo=(180./rpi)*abs(rlatgeo(i)) + IF ( k .LE. 4 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=0.07 + END IF + IF ( k .EQ. 5 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=0.1 + IF (latgeo.GE.20.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.09 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.07 + END IF + IF ( k .EQ. 6 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=0.14 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.12 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.1 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.09 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.07 + END IF + IF ( k .EQ. 7 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=0.16 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.16 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.14 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.12 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.1 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.09 + END IF + IF ( k .EQ. 8 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=0.175 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.16 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.14 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.12 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.1 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.1 + END IF + IF ( k .EQ. 9 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=0.245 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.21 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.175 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.14 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.12 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.12 + END IF + IF ( k .EQ. 10 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=0.31 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.28 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.245 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.21 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.16 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.14 + END IF + IF ( k .EQ. 11 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=0.35 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.3 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.3 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.2 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.18 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.16 + END IF + IF ( k .EQ. 12 ) THEN + IF (latgeo.GE.40.0) srcbe(i,k)=0.5 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.4 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.35 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.3 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.25 + END IF + IF ( k .EQ. 13 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=0.8 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.7 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.6 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.5 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.4 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.35 + END IF + IF ( k .EQ. 14 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=1.2 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=1. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.75 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.6 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.5 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.4 + END IF + IF ( k .EQ. 15 ) THEN + IF (latgeo.GE.60.0) srcbe(i,k)=1.75 + IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=1.8 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=1.6 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=1.4 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.9 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.75 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.65 + END IF + IF ( k .EQ. 16 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=3. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=2.5 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=1.8 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=1.5 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=1.2 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.9 + END IF + IF ( k .EQ. 17 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=4. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=3. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=2.5 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=2. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=1.6 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=1.4 + END IF + IF ( k .EQ. 18 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=7. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=6. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=4.5 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=3.5 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=3. + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=2. + END IF + IF ( k .EQ. 19 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=8.5 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=8. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=7. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=4. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=3.5 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=3. + END IF + IF ( k .EQ. 20 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=12.5 + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=12. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=8. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=6. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=4. + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=3.5 + END IF + IF ( k .EQ. 21 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=16. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=13. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=10. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=7.5 + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=4.5 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=4. + END IF + IF ( k .EQ. 22 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=20. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=17.5 + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=12.5 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=9. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=6. + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=4.5 + END IF + IF ( k .EQ. 23 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=25. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=22. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=15. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=10. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=7.5 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=6. + END IF + IF ( k .EQ. 24 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=28. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=26. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=18. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=12. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=8.5 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7. + END IF + IF ( k .EQ. 25 ) THEN + IF (latgeo.GE.50.0) srcbe(i,k)=33. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=28. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=20. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=14. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=10. + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=8.5 + END IF + IF ( k .EQ. 26 ) THEN + IF (latgeo.GE.60.0) srcbe(i,k)=38. + IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=36. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=32. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=24. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=15. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=11.5 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10. + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=6. + END IF + IF ( k .EQ. 27 ) THEN + IF (latgeo.GE.60.0) srcbe(i,k)=46. + IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=44. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=35. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=25. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=16. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=12.5 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10. + END IF + IF ( k .EQ. 28 ) THEN + IF (latgeo.GE.60.0) srcbe(i,k)=53. + IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=48. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=37. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=25. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=16. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=12.5 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10. + END IF + IF ( k .EQ. 29 ) THEN + IF (latgeo.GE.70.0) srcbe(i,k)=58. + IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=56. + IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=50. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=36. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=24. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=15. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=11.5 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10. + END IF + IF ( k .EQ. 30 ) THEN + IF (latgeo.GE.70.0) srcbe(i,k)=65. + IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=60. + IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=50. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=35. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=22. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=14. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=10. + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=9. + END IF + IF ( k .EQ. 31 ) THEN + IF (latgeo.GE.70.0) srcbe(i,k)=70. + IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=62. + IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=48. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=32. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=21. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=13. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=9. + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.6 + END IF + IF ( k .EQ. 32 ) THEN + IF (latgeo.GE.70.0) srcbe(i,k)=80. + IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=60. + IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=46. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=30. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=17.5 + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=11. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=8. + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.4 + END IF + IF ( k .GE. 33 ) THEN + IF (latgeo.GE.70.0) srcbe(i,k)=80. + IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=70. + IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=45. + IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=27. + IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=15. + IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=10. + IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=7.6 + IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7. + END IF + END DO + END DO + END IF ! fin de 39 niveaux verticaux + + +!==================================== +! Conversion de la source en U/s/kgA +!==================================== + DO k = 1,klev + DO i = 1,klon + ! La source est at/min/m3 -> at/s/kgA + ! avec une masse volumique de l'air = 1.295 kg/m3 + ! 1/(60*1.295) = 0.01287 + srcbe(i,k)=srcbe(i,k)*0.01287 + ! La source est at/min/m3 -> at/s/m3 + ! srcbe(i,k)=srcbe(i,k)*0.0166667 + END DO + END DO + +END SUBROUTINE init_be Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/initrrnpb.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/initrrnpb.F (revision 1190) +++ (revision ) @@ -1,125 +1,0 @@ -! -! $Header$ -! - SUBROUTINE initrrnpb(ftsol,pctsrf,masktr,fshtr,hsoltr,tautr - . ,vdeptr,scavtr) - USE dimphy - USE infotrac, ONLY : nbtr - IMPLICIT none -c====================================================================== -c Auteur(s): AA + CG (LGGE/CNRS) Date 24-06-94 -c Objet: initialisation des constantes des traceurs -CAA Revison pour le controle avec la temperature du sol -cAA -CAA it = 1 radon ss controle de ts -cAA it = 2 plomb ss controle de ts -c====================================================================== -c Arguments: -c nbtr------input-I- nombre de vrais traceurs (sans l'eau) -c ftsol-------input-R- Temperature du sol (Kelvin) -c pctsrf-----input-R- Nature de sol (pourcentage de sol) -c masktr---output-R- Masque reservoir de sol traceur (1 = reservoir) -c fshtr----output-R- Flux surfacique de production dans le sol -c hsoltr---output-R- Epaisseur du reservoir de sol -c tautr----output-R- Constante de decroissance du traceur -c vdeptr---output-R- Vitesse de depot sec dans la couche Brownienne -c scavtr---output-R- Coefficient de lessivage -c====================================================================== -cym#include "dimensions.h" -cym#include "dimphy.h" -#include "indicesol.h" -c====================================================================== -C - INTEGER i, it - REAL pctsrf(klon,nbsrf) !Pourcentage de sol (f(nature du sol)) - REAL ftsol(klon,nbsrf) ! Temperature du sol pour le controle Rn -c ! le cas echeant - REAL masktr(klon,nbtr) ! Masque de l'echange avec la surface -c (possible => 1 ) - REAL fshtr(klon,nbtr) ! Flux surfacique dans le reservoir de sol - REAL hsoltr(nbtr) ! Epaisseur equivalente du reservoir de sol - REAL tautr(nbtr) ! Constante de decroissance radioactive - REAL vdeptr(nbtr) ! Vitesse de depot sec dans la couche Brownienne - REAL scavtr(nbtr) ! Coefficient de lessivage - REAL s -C - WRITE(*,'(''PASSAGE initrrnpb ...'',$)') - print*,'nbtr= ',nbtr - print*,'nbsrf= ',nbsrf - print*,'klon= ',klon -C -C Puis les initialisation specifiques a chaque traceur (pour le moment, Rn222) -C -C -C Radon it = 1 -c - IF ( nbtr .LE. 0 ) STOP 'initrrnpb pas glop pas glop' - it = 1 - s = 1.E4 ! Source: atome par m2 - hsoltr(it) = 0.1 ! Hauteur equivalente du reservoir : -c 1 m * porosite 0.1 - tautr(it) = 4.765E5 ! Decroissance du radon, secondes -cAA -c tautr(it) = 4.765E55 ! Decroissance du radon,infinie -cAA - vdeptr(it) = 0. ! Pas de depot sec pour le radon - scavtr(it) = 0. ! Pas de lessivage pour le radon - - print*, '-------------- SOURCE DU RADON ------------------------ ' - print*,'it = ',it - print*,'Source : ', s - print*,'Hauteur equivalente du reservoir de sol: ',hsoltr(it) - print*,'Decroissance (s): ', tautr(it) - print*,'Vitesse de depot sec: ',vdeptr(it) - print*,'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) - -cAA -cAA quelques tests -cAA POur l'instant le pctsrf(i,3) = 1.0 -cAA lorsqu'il ya de la terre mias ne prend aucune autre valeur -cAA il n'est donc pas necessaire de multiplier fshtr par pctsrf -cAA -c print*, '------------------------------------------ ' -c print*, 'masktr(',i,it,')= ',masktr(i,it) -c print*, 'fshtr(',i,it,')= ',fshtr(i,it) -c print*, 'pctsrf(',i,',1)= ',pctsrf(i,1) -c print*, 'pctsrf(',i,',2)= ',pctsrf(i,2) -c print*, 'pctsrf(',i,',3)= ',pctsrf(i,3) -c print*, 'pctsrf(',i,',4)= ',pctsrf(i,4) -c print*, 's = ',s -c print*, '------------------------------------------ ' - - END DO -C -C 210Pb it = 2 -C - IF ( nbtr .LE. 1 ) STOP 'initrrnpb pas glop pas glop' - it = 2 - s = 0. ! Pas de source !!! - hsoltr(it) = 10. ! Hauteur equivalente du reservoir -c a partir duquel le -c 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 - print*, '-------------- SOURCE DU PLOMB ------------------------ ' - print*,'it = ',it - print*,'Source : ', s - print*,'Hauteur equivalente du reservoir : ',hsoltr(it) - print*,'Decroissance (s): ', tautr(it) - print*,'Vitesse de depot sec: ',vdeptr(it) - print*,'Facteur de lessivage: ',scavtr(it) -c - WRITE(*,*) 'initialisation rnpb ok' -c - RETURN - END Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/initrrnpb.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/initrrnpb.F90 (revision 1191) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/initrrnpb.F90 (revision 1191) @@ -0,0 +1,92 @@ +! +! $Id $ +! +SUBROUTINE initrrnpb(ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr) + USE dimphy + USE infotrac, ONLY : nbtr + 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 +!====================================================================== +! Arguments: +! nbtr.............. nombre de vrais traceurs (sans l'eau) +! ftsol....input-R- Temperature du sol (Kelvin) +! pctsrf...input-R- Nature de sol (pourcentage de sol) +! masktr...output-R- Masque reservoir de sol traceur (1 = reservoir) +! fshtr....output-R- Flux surfacique de production dans le reservoir de sol +! hsoltr...output-R- Epaisseur equivalente du reservoir de sol +! tautr....output-R- Constante de decroissance radioactive du traceur +! vdeptr...output-R- Vitesse de depot sec dans la couche Brownienne +! scavtr...output-R- Coefficient de lessivage +!====================================================================== + INCLUDE "indicesol.h" +!====================================================================== + + REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf + REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol + REAL,DIMENSION(klon,nbtr),INTENT(OUT) :: masktr + REAL,DIMENSION(klon,nbtr),INTENT(OUT) :: fshtr + REAL,DIMENSION(nbtr),INTENT(OUT) :: hsoltr + REAL,DIMENSION(nbtr),INTENT(OUT) :: tautr + REAL,DIMENSION(nbtr),INTENT(OUT) :: vdeptr + REAL,DIMENSION(nbtr),INTENT(OUT) :: scavtr + INTEGER :: i, it + REAL :: s + + WRITE(*,*)'PASSAGE initrrnpb ...' +! +! Radon it = 1 +!---------------- + IF ( nbtr .LE. 0 ) STOP '**PHYTRAC:initrrnpb:** nbtr < 0; verifier RN dans traceur.def' + 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 ) STOP '**PHYTRAC**:initrrnpb:** nbtr <= 1; verifier PB dans traceur.def' + 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' + +END SUBROUTINE initrrnpb Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/minmaxqfi.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/minmaxqfi.F (revision 1190) +++ (revision ) @@ -1,34 +1,0 @@ -! -! $Header$ -! - SUBROUTINE minmaxqfi(zq,qmin,qmax,comment) - USE dimphy - IMPLICIT none - -cym#include "dimensions.h" -cym#include "dimphy.h" - - CHARACTER*(*) comment - real qmin,qmax - real zq(klon,klev) - - INTEGER jadrs(klon), jbad, k, i - - DO k = 1, klev - jbad = 0 - DO i = 1, klon - IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN - jbad = jbad + 1 - jadrs(jbad) = i - ENDIF - ENDDO - IF (jbad.GT.0) THEN - PRINT*, comment - DO i = 1, jbad - PRINT*, "i,k,q=", jadrs(i),k,zq(jadrs(i),k) - ENDDO - ENDIF - ENDDO - - return - end Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/minmaxqfi.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/minmaxqfi.F90 (revision 1191) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/minmaxqfi.F90 (revision 1191) @@ -0,0 +1,33 @@ +! +! $Id$ +! +SUBROUTINE minmaxqfi(zq,qmin,qmax,comment) + USE dimphy + IMPLICIT NONE + +! Entrees + REAL,DIMENSION(klon,klev), INTENT(IN) :: zq + REAL,INTENT(IN) :: qmin,qmax + CHARACTER(LEN=*),INTENT(IN) :: comment + +! Local + INTEGER,DIMENSION(klon) :: jadrs + INTEGER :: i, jbad, k + + DO k = 1, klev + jbad = 0 + DO i = 1, klon + IF (zq(i,k).GT.qmax .OR. zq(i,k).LT.qmin) THEN + jbad = jbad + 1 + jadrs(jbad) = i + ENDIF + ENDDO + IF (jbad.GT.0) THEN + WRITE(*,*)comment + DO i = 1, jbad + WRITE(*,*) "i,k,q=", jadrs(i),k,zq(jadrs(i),k) + ENDDO + ENDIF + ENDDO + +END SUBROUTINE minmaxqfi Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/nflxtr.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/nflxtr.F (revision 1190) +++ (revision ) @@ -1,166 +1,0 @@ -! -! $Header$ -! - SUBROUTINE nflxtr(pdtime,pmfu,pmfd,pen_u,pde_u,pen_d,pde_d, - . pplay,paprs,x,dx) - USE dimphy - IMPLICIT NONE -c===================================================================== -c Objet : Melange convectif de traceurs a partir des flux de masse -c Date : 13/12/1996 -- 13/01/97 -c Auteur: O. Boucher (LOA) sur inspiration de Z. X. Li (LMD), -c Brinkop et Sausen (1996) et Boucher et al. (1996). -c ATTENTION : meme si cette routine se veut la plus generale possible, -c elle a herite de certaines notations et conventions du -c schema de Tiedtke (1993). -c --En particulier, les couches sont numerotees de haut en bas !!! -c Ceci est valable pour les flux -c mais pas pour les entrees x, pplay, paprs !!!! -c --pmfu est positif, pmfd est negatif -c --Tous les flux d'entrainements et de detrainements sont positifs -c contrairement au schema de Tiedtke d'ou les changements de signe!!!! -c===================================================================== -c -cym#include "dimensions.h" -cym#include "dimphy.h" -#include "YOMCST.h" -#include "YOECUMF.h" -c - REAL pdtime -c--les flux sont definis au 1/2 niveaux -c--pmfu(klev+1) et pmfd(klev+1) sont implicitement nuls - REAL pmfu(klon,klev) ! flux de masse dans le panache montant - REAL pmfd(klon,klev) ! flux de masse dans le panache descendant - REAL pen_u(klon,klev) ! flux entraine dans le panache montant - REAL pde_u(klon,klev) ! flux detraine dans le panache montant - REAL pen_d(klon,klev) ! flux entraine dans le panache descendant - REAL pde_d(klon,klev) ! flux detraine dans le panache descendant - - REAL pplay(klon,klev) ! pression aux couches (bas en haut) - REAL paprs(klon,klev+1) ! pression aux 1/2 couches (bas en haut) - REAL x(klon,klev) ! q de traceur (bas en haut) - REAL dx(klon,klev) ! tendance de traceur (bas en haut) -c -c--flux convectifs mais en variables locales - REAL zmfu(klon,klev+1) - REAL zmfd(klon,klev+1) - REAL zen_u(klon,klev) - REAL zde_u(klon,klev) - REAL zen_d(klon,klev) - REAL zde_d(klon,klev) - real zmfe -c -c--variables locales -c--les flux de x sont definis aux 1/2 niveaux -c--xu et xd sont definis aux niveaux complets - REAL xu(klon,klev) ! q de traceurs dans le panache montant - REAL xd(klon,klev) ! q de traceurs dans le panache descendant - REAL zmfux(klon,klev+1) ! flux de x dans le panache montant - REAL zmfdx(klon,klev+1) ! flux de x dans le panache descendant - REAL zmfex(klon,klev+1) ! flux de x dans l'environnement - INTEGER i, k - REAL zmfmin - PARAMETER (zmfmin=1.E-10) - -c ========================================= -c -c -c Extension des flux UP et DN sur klev+1 niveaux -c ========================================= - do k=1,klev - do i=1,klon - zmfu(i,k)=pmfu(i,k) - zmfd(i,k)=pmfd(i,k) - enddo - enddo - do i=1,klon - zmfu(i,klev+1)=0. - zmfd(i,klev+1)=0. - enddo - -c--modif pour diagnostiquer les detrainements -c ========================================= -c on privilegie l'ajustement de l'entrainement dans l'ascendance. - - do k=1, klev - do i=1, klon - zen_d(i,k)=pen_d(i,k) - zde_u(i,k)=pde_u(i,k) - zde_d(i,k) =-zmfd(i,k+1)+zmfd(i,k)+zen_d(i,k) - zen_u(i,k) = zmfu(i,k+1)-zmfu(i,k)+zde_u(i,k) - enddo - enddo -c -c--calcul des flux dans le panache montant -c ========================================= -c -c Dans la premiere couche, on prend q comme valeur de qu -c - do i=1, klon - zmfux(i,1)=0.0 - enddo -c -c Autres couches - do k=1,klev - do i=1, klon - if ((zmfu(i,k+1)+zde_u(i,k)).lt.zmfmin) THEN - xu(i,k)=x(i,k) - else - xu(i,k)=(zmfux(i,k)+zen_u(i,k)*x(i,k)) - s /(zmfu(i,k+1)+zde_u(i,k)) - endif - zmfux(i,k+1)=zmfu(i,k+1)*xu(i,k) - enddo - enddo -c -c--calcul des flux dans le panache descendant -c ========================================= -c - do i=1, klon - zmfdx(i,klev+1)=0.0 - enddo -c - do k=klev,1,-1 - do i=1, klon - if ((zde_d(i,k)-zmfd(i,k)).lt.zmfmin) THEN - xd(i,k)=x(i,k) - else - xd(i,k)=(zmfdx(i,k+1)-zen_d(i,k)*x(i,k)) / - . (zmfd(i,k)-zde_d(i,k)) - endif - zmfdx(i,k)=zmfd(i,k)*xd(i,k) - enddo - enddo -c -c--introduction du flux de retour dans l'environnement -c ========================================= -c - do k=2, klev - do i=1, klon - zmfe=-zmfu(i,k)-zmfd(i,k) - if (zmfe.le.0.) then - zmfex(i,k)= zmfe*x(i,k) - else - zmfex(i,k)= zmfe*x(i,k-1) - endif - enddo - enddo - - do i=1, klon - zmfex(i,1)=0. - zmfex(i,klev+1)=0. - enddo -c -c--calcul final des tendances -c - do k=1, klev - do i=1, klon - dx(i,k)=RG/(paprs(i,k)-paprs(i,k+1))*pdtime* - . ( zmfux(i,k) - zmfux(i,k+1) + - . zmfdx(i,k) - zmfdx(i,k+1) + - . zmfex(i,k) - zmfex(i,k+1) ) - enddo - enddo -c - return - end Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/nflxtr.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/nflxtr.F90 (revision 1191) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/nflxtr.F90 (revision 1191) @@ -0,0 +1,159 @@ +! +! $Id $ +! +SUBROUTINE nflxtr(pdtime,pmfu,pmfd,pen_u,pde_u,pen_d,pde_d,pplay,paprs,x,dx) + USE dimphy + IMPLICIT NONE +!===================================================================== +! Objet : Melange convectif de traceurs a partir des flux de masse +! Date : 13/12/1996 -- 13/01/97 +! Auteur: O. Boucher (LOA) sur inspiration de Z. X. Li (LMD), +! Brinkop et Sausen (1996) et Boucher et al. (1996). +! ATTENTION : meme si cette routine se veut la plus generale possible, +! elle a herite de certaines notations et conventions du +! schema de Tiedtke (1993). +! 1. En particulier, les couches sont numerotees de haut en bas !!! +! Ceci est valable pour les flux +! mais pas pour les entrees x, pplay, paprs !!!! +! 2. pmfu est positif, pmfd est negatif +! 3. Tous les flux d'entrainements et de detrainements sont positifs +! contrairement au schema de Tiedtke d'ou les changements de signe!!!! +!===================================================================== +! + include "YOMCST.h" + include "YOECUMF.h" + + REAL,INTENT(IN) :: pdtime ! pdtphys +! +! les flux sont definis au 1/2 niveaux +! => pmfu(klev+1) et pmfd(klev+1) sont implicitement nuls +! + REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfu ! flux de masse dans le panache montant + REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfd ! flux de masse dans le panache descendant + REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_u ! flux entraine dans le panache montant + REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_u ! flux detraine dans le panache montant + REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_d ! flux entraine dans le panache descendant + REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_d ! flux detraine dans le panache descendant + + REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression aux couches (bas en haut) + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression aux 1/2 couches (bas en haut) + REAL,DIMENSION(klon,klev),INTENT(IN) :: x ! q de traceur (bas en haut) + REAL,DIMENSION(klon,klev),INTENT(INOUT) :: dx ! tendance de traceur (bas en haut) + +! flux convectifs mais en variables locales + REAL,DIMENSION(klon,klev+1) :: zmfu ! copie de pmfu avec klev+1 = 0 + REAL,DIMENSION(klon,klev+1) :: zmfd ! copie de pmfd avec klev+1 = 0 + REAL,DIMENSION(klon,klev) :: zen_u + REAL,DIMENSION(klon,klev) :: zde_u + REAL,DIMENSION(klon,klev) :: zen_d + REAL,DIMENSION(klon,klev) :: zde_d + REAL :: zmfe + +! variables locales +! les flux de x sont definis aux 1/2 niveaux +! xu et xd sont definis aux niveaux complets + REAL,DIMENSION(klon,klev) :: xu ! q de traceurs dans le panache montant + REAL,DIMENSION(klon,klev) :: xd ! q de traceurs dans le panache descendant + REAL,DIMENSION(klon,klev+1) :: zmfux ! flux de x dans le panache montant + REAL,DIMENSION(klon,klev+1) :: zmfdx ! flux de x dans le panache descendant + REAL,DIMENSION(klon,klev+1) :: zmfex ! flux de x dans l'environnement + INTEGER :: i, k + REAL,PARAMETER :: zmfmin=1.E-10 + +! ============================================== +! Extension des flux UP et DN sur klev+1 niveaux +! ============================================== + DO k=1,klev + DO i=1,klon + zmfu(i,k)=pmfu(i,k) + zmfd(i,k)=pmfd(i,k) + ENDDO + ENDDO + DO i=1,klon + zmfu(i,klev+1)=0. + zmfd(i,klev+1)=0. + ENDDO +! ========================================== +! modif pour diagnostiquer les detrainements +! ========================================== +! on privilegie l'ajustement de l'entrainement dans l'ascendance. + + DO k=1, klev + DO i=1, klon + zen_d(i,k)=pen_d(i,k) + zde_u(i,k)=pde_u(i,k) + zde_d(i,k) =-zmfd(i,k+1)+zmfd(i,k)+zen_d(i,k) + zen_u(i,k) = zmfu(i,k+1)-zmfu(i,k)+zde_u(i,k) + ENDDO + ENDDO +! ========================================= +! calcul des flux dans le panache montant +! ========================================= +! +! Dans la premiere couche, on prend q comme valeur de qu + + DO i=1, klon + zmfux(i,1)=0.0 + ENDDO + +! Autres couches + DO k=1,klev + DO i=1, klon + IF ((zmfu(i,k+1)+zde_u(i,k)).lt.zmfmin) THEN + xu(i,k)=x(i,k) + ELSE + xu(i,k)=(zmfux(i,k)+zen_u(i,k)*x(i,k))/(zmfu(i,k+1)+zde_u(i,k)) + ENDIF + zmfux(i,k+1)=zmfu(i,k+1)*xu(i,k) + ENDDO + ENDDO +! ========================================== +! calcul des flux dans le panache descendant +! ========================================== + + DO i=1, klon + zmfdx(i,klev+1)=0.0 + ENDDO + + DO k=klev,1,-1 + DO i=1, klon + IF ((zde_d(i,k)-zmfd(i,k)).lt.zmfmin) THEN + xd(i,k)=x(i,k) + ELSE + xd(i,k)=(zmfdx(i,k+1)-zen_d(i,k)*x(i,k))/(zmfd(i,k)-zde_d(i,k)) + ENDIF + zmfdx(i,k)=zmfd(i,k)*xd(i,k) + ENDDO + ENDDO +! =================================================== +! introduction du flux de retour dans l'environnement +! =================================================== + + DO k=2, klev + DO i=1, klon + zmfe=-zmfu(i,k)-zmfd(i,k) + IF (zmfe.le.0.) then + zmfex(i,k)= zmfe*x(i,k) + ELSE + zmfex(i,k)= zmfe*x(i,k-1) + ENDIF + ENDDO + ENDDO + + DO i=1, klon + zmfex(i,1)=0. + zmfex(i,klev+1)=0. + ENDDO +! ========================== +! calcul final des tendances +! ========================== + DO k=1, klev + DO i=1, klon + dx(i,k)=RG/(paprs(i,k)-paprs(i,k+1))*pdtime* & + ( zmfux(i,k) - zmfux(i,k+1) + & + zmfdx(i,k) - zmfdx(i,k+1) + & + zmfex(i,k) - zmfex(i,k+1) ) + ENDDO + ENDDO + +END SUBROUTINE nflxtr Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/phyetat0.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/phyetat0.F (revision 1190) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/phyetat0.F (revision 1191) @@ -19,4 +19,7 @@ USE iostart USE write_field_phy + USE infotrac + USE traclmdz_mod, ONLY : traclmdz_from_restart + IMPLICIT none c====================================================================== @@ -48,4 +51,5 @@ REAL run_off_lic_0(klon) REAL fractint(klon) + REAL trs(klon,nbtr) CHARACTER*6 ocean_in @@ -62,4 +66,5 @@ INTEGER length PARAMETER (length=100) + INTEGER it, iiq REAL tab_cntrl(length), tabcntr0(length) CHARACTER*7 str7 @@ -983,4 +988,27 @@ PRINT*,'(ecart-type) wake_fip:', xmin, xmax c +c Read and send field trs to traclmdz +c + IF (type_trac == 'lmdz') THEN + DO it=1,nbtr + iiq=niadv(it+2) + CALL get_field("trs_"//tname(iiq),trs(:,it),found) + IF (.NOT. found) THEN + PRINT*, + $ "phyetat0: Le champ est absent" + PRINT*, "Depart legerement fausse. Mais je continue" + trs(:,it) = 0. + ENDIF + xmin = 1.0E+20 + xmax = -1.0E+20 + xmin = MINval(trs(:,it)) + xmax = MAXval(trs(:,it)) + PRINT*,"(ecart-type) trs_"//tname(iiq)//" :", xmin, xmax + + END DO + + CALL traclmdz_from_restart(trs) + END IF + c on ferme le fichier @@ -1005,4 +1033,5 @@ CALL fonte_neige_init(run_off_lic_0) + RETURN END Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/phyredem.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/phyredem.F (revision 1190) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/phyredem.F (revision 1191) @@ -12,4 +12,6 @@ USE phys_state_var_mod USE iostart + USE traclmdz_mod, ONLY : traclmdz_to_restart + USE infotrac IMPLICIT none @@ -42,4 +44,5 @@ REAL agesno(klon,nbsrf) REAL run_off_lic_0(klon) + REAL trs(klon,nbtr) c INTEGER nid, nvarid, idim1, idim2, idim3 @@ -52,5 +55,6 @@ CHARACTER (len=7) :: str7 CHARACTER (len=2) :: str2 - + INTEGER :: it, iiq + c====================================================================== c @@ -311,4 +315,14 @@ CALL put_field("WAKE_FIP","",wake_fip) + +! trs from traclmdz_mod + IF (type_trac == 'lmdz') THEN + CALL traclmdz_to_restart(trs) + DO it=1,nbtr + iiq=niadv(it+2) + CALL put_field("trs_"//tname(iiq),"",trs(:,it)) + END DO + END IF + CALL close_restartphy !$OMP BARRIER Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F (revision 1190) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F (revision 1191) @@ -111,6 +111,4 @@ PARAMETER (ok_stratus=.FALSE.) c====================================================================== - LOGICAL, SAVE :: rnpb=.TRUE. -c$OMP THREADPRIVATE(rnpb) REAL amn, amx INTEGER igout @@ -3133,67 +3131,22 @@ c==================================================================== C - IF (config_inca /= 'none') rnpb=.FALSE. - - call phytrac ( rnpb, - I itap, - I julien, - I gmtime, - I debut, - I lafin, - I nlon, - I klev, - I dtime, - I u, - I v, - I t, - I paprs, - I pplay, - I pmfu, - I pmfd, - I pen_u, - I pde_u, - I pen_d, - I pde_d, - I cdragh, - I coefh, - I fm_therm, - I entr_therm, - I u1, - I v1, - I ftsol, - I pctsrf, - I rlat, - I frac_impa, - I frac_nucl, - I rlon, - I presnivs, - I pphis, - I pphi, - I albsol1, - I qx(1,1,1), - I rhcl, - I cldfra, - I rneb, - I diafra, - I cldliq, - I itop_con, - I ibas_con, - I pmflxr, - I pmflxs, - I prfl, - I psfl, - I da, - I phi, - I mp, - I upwd, - I dnwd, - I aerosol_couple, - I flxmass_w, - I tau_aero, - I piz_aero, - I cg_aero, - I ccm, - I rfname, - O tr_seri) + + call phytrac ( + I itap, julien, gmtime, debut, + I lafin, dtime, u, v, t, + I paprs, pplay, pmfu, pmfd, + I pen_u, pde_u, pen_d, pde_d, + I cdragh, coefh, fm_therm, entr_therm, + I u1, v1, ftsol, pctsrf, + I rlat, frac_impa, frac_nucl,rlon, + I presnivs, pphis, pphi, albsol1, + I qx(:,:,ivap),rhcl, cldfra, rneb, + I diafra, cldliq, itop_con, ibas_con, + I pmflxr, pmflxs, prfl, psfl, + I da, phi, mp, upwd, + I dnwd, aerosol_couple, flxmass_w, + I tau_aero, piz_aero, cg_aero, ccm, + I rfname, + O tr_seri) IF (offline) THEN Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F (revision 1190) +++ (revision ) @@ -1,833 +1,0 @@ -! -c -c - SUBROUTINE phytrac (rnpb, - I nstep, - I julien, - I gmtime, - I debutphy, - I lafin, - I nlon, - I nlev, - I pdtphys, - I u, - I v, - I t_seri, - I paprs, - I pplay, - I pmfu, - I pmfd, - I pen_u, - I pde_u, - I pen_d, - I pde_d, - I cdragh, - I coefh, - I fm_therm, - I entr_therm, - I yu1, - I yv1, - I ftsol, - I pctsrf, - I xlat, - I frac_impa, - I frac_nucl, - I xlon, - I presnivs, - I pphis, - I pphi, - I albsol, - I sh, - I rh, - I cldfra, - I rneb, - I diafra, - I cldliq, - I itop_con, - I ibas_con, - I pmflxr, - I pmflxs, - I prfl, - I psfl, - I da, - I phi, - I mp, - I upwd, - I dnwd, - I aerosol_couple, - I flxmass_w, - I tau_aero, - I piz_aero, - I cg_aero, - I ccm, - I rfname, - O tr_seri) - - USE ioipsl - USE dimphy - USE infotrac - USE mod_grid_phy_lmdz - USE mod_phys_lmdz_para - USE comgeomphy - USE iophy - USE vampir - - IMPLICIT none -c====================================================================== -c Auteur(s) FH -c Objet: Moniteur general des tendances traceurs -c -cAA Remarques en vrac: -cAA-------------------- -cAA 2/ Le choix du radon et du pb se fait juste avec un data -cAA (peu propre). Peut-etre pourrait-on prevoir dans l'avenir -cAA une variable "type de traceur" -c====================================================================== -#include "YOMCST.h" -#include "dimensions.h" -#include "indicesol.h" -#include "clesphys.h" -#include "temps.h" -#include "paramet.h" -#include "control.h" -#include "thermcell.h" -c====================================================================== - -c Arguments: -c -c EN ENTREE: -c ========== -c -c divers: -c ------- -c - integer nlon ! nombre de points horizontaux - integer nlev ! nombre de couches verticales - integer nstep ! appel physique - integer julien !jour julien - integer itop_con(nlon) - integer ibas_con(nlon) - real gmtime - real pdtphys ! pas d'integration pour la physique (seconde) - real t_seri(nlon,nlev) ! temperature - real tr_seri(nlon,nlev,nbtr) ! traceur - real u(klon,klev) - real v(klon,klev) - real sh(nlon,nlev) ! humidite specifique - real rh(nlon,nlev) ! humidite relative - real cldliq(nlon,nlev) ! eau liquide nuageuse - real cldfra(nlon,nlev) ! fraction nuageuse (tous les nuages) - real diafra(nlon,nlev) ! fraction nuageuse (convection ou stratus artificiels) - real rneb(nlon,nlev) ! fraction nuageuse (grande echelle) - real albsol(nlon) ! albedo surface - real paprs(nlon,nlev+1) ! pression pour chaque inter-couche (en Pa) - real ps(nlon) ! pression surface - real pplay(nlon,nlev) ! pression pour le mileu de chaque couche (en Pa) - real pphi(nlon,klev) ! geopotentiel - real pphis(klon) - REAL presnivs(klev) - logical debutphy ! le flag de l'initialisation de la physique - logical lafin ! le flag de la fin de la physique -c Olivia - integer nsplit - REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1) !--lessivage convection - REAL prfl(klon,klev+1), psfl(klon,klev+1) !--lessivage large-scale - LOGICAL aerosol_couple - - REAL flxmass_w(klon,klev) - CHARACTER(len=8) :: solsym(nbtr) - integer la - REAL :: tau_aero(klon,klev,9,2) - REAL :: piz_aero(klon,klev,9,2) - REAL :: cg_aero(klon,klev,9,2) - character*4 :: rfname(9) - REAL :: ccm(klon,klev,2) -c -c convection: -c ----------- -c - REAL pmfu(nlon,nlev) ! flux de masse dans le panache montant - REAL pmfd(nlon,nlev) ! flux de masse dans le panache descendant - REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant - -c -c thermiques: -c ----------- -c - real fm_therm(klon,klev+1),entr_therm(klon,klev) - real fm_therm1(klon,klev) -c - REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant - REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant - REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant -c KE - real da(nlon,nlev),phi(nlon,nlev,nlev),mp(nlon,nlev) - REAL upwd(nlon,nlev) ! saturated updraft mass flux - REAL dnwd(nlon,nlev) ! saturated downdraft mass flux - -c -c Couche limite: -c -------------- -c - REAL cdragh(nlon,nlev)! coeff drag pour T et Q - REAL coefh(nlon,nlev) ! coeff melange CL - REAL yu1(nlon) ! vents au premier niveau - REAL yv1(nlon) ! vents au premier niveau - REAL xlat(nlon) ! latitudes pour chaque point - REAL xlon(nlon) ! longitudes pour chaque point - -c -c Lessivage: -c ---------- -c -c pour le ON-LINE -c - REAL frac_impa(nlon,nlev) ! fraction d'aerosols impactes - REAL frac_nucl(nlon,nlev) ! fraction d'aerosols nuclees -c -cAA -cAA Arguments necessaires pour les sources et puits de traceur: -cAA ---------------- -cAA - real ftsol(nlon,nbsrf) ! Temperature du sol (surf)(Kelvin) - real pctsrf(nlon,nbsrf) ! Pourcentage de sol f(nature du sol) -c abder - real pftsol1(nlon),pftsol2(nlon),pftsol3(nlon),pftsol4(nlon) - real ppsrf1(nlon),ppsrf2(nlon),ppsrf3(nlon),ppsrf4(nlon) -c fin -cAA ---------------------------- -cAA VARIABLES LOCALES TRACEURS -cAA ---------------------------- -cAA -cAA Sources et puits des traceurs: -cAA ------------------------------ -cAA -cAA Pour l'instant seuls les cas du rn et du pb ont ete envisages. - - REAL source(klon,nbtr) ! a voir lorsque le flux est prescrit -cAA -cAA Pour la source de radon et son reservoir de sol -cAA ................................................ - - REAL,save,allocatable :: trs(:,:) ! Conc. radon ds le sol -c$OMP THREADPRIVATE(trs) - REAL :: trs_tmp(klon_glo) - - REAL,save,allocatable :: masktr(:,:) ! Masque reservoir de sol traceur -c Masque de l'echange avec la surface -c (1 = reservoir) ou (possible => 1 ) -c$OMP THREADPRIVATE(masktr) - REAL,save,allocatable :: fshtr(:,:) ! Flux surfacique dans le reservoir de sol -c$OMP THREADPRIVATE(fshtr) - REAL,save,allocatable :: hsoltr(:) ! Epaisseur equivalente du reservoir de sol -c$OMP THREADPRIVATE(hsoltr) - REAL,save,allocatable :: tautr(:) ! Constante de decroissance radioactive -c$OMP THREADPRIVATE(tautr) - REAL,save,allocatable :: vdeptr(:) ! Vitesse de depot sec dans la couche Brownienne -c$OMP THREADPRIVATE(vdeptr) - REAL,save,allocatable :: scavtr(:) ! Coefficient de lessivage -c$OMP THREADPRIVATE(scavtr) -cAA - CHARACTER*2 itn -C maf ioipsl - CHARACTER*2 str2 - INTEGER nhori, nvert - REAL zsto, zout, zjulian - INTEGER, SAVE :: nid_tra -c$OMP THREADPRIVATE(nid_tra) - INTEGER, SAVE :: nid_tra2,nid_tra3 -c$OMP THREADPRIVATE(nid_tra2,nid_tra3) - INTEGER ndex(1) - INTEGER ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev) - REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique - REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D - REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev) - REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1) -c - integer itau_w ! pas de temps ecriture = nstep + itau_phy -c - logical ok_sync - parameter (ok_sync = .true.) -C -C nature du traceur -c - logical,save,allocatable :: aerosol(:) ! Nature du traceur -c ! aerosol(it) = true => aerosol -c ! aerosol(it) = false => gaz -c ! nat_trac(it) = 1. aerosol - logical,save,allocatable :: clsol(:) ! clsol(it) = true => CL sol calculee - logical,save,allocatable :: radio(:) ! radio(it)=true => decroisssance radioactive -c$OMP THREADPRIVATE(aerosol,clsol,radio) -C -c====================================================================== -c -c Declaration des procedures appelees -c -c--modif convection tiedtke - INTEGER i, k, it - INTEGER iq, iiq - REAL delp(klon,klev) -c--end modif -c -c Variables liees a l'ecriture de la bande histoire physique -c -c Variables locales pour effectuer les appels en serie -c---------------------------------------------------- -c - REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs - REAL d_tr_cl(klon,klev,nbtr) ! tendance de traceurs couche limite - REAL d_tr_cv(klon,klev,nbtr) ! tendance de traceurs conv pour chq traceur - REAL d_tr_th(klon,klev,nbtr) ! la tendance des thermiques - REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance -c ! radioactive du rn - > pb - REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage -c ! par impaction - REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage -c ! par nucleation - REAL fluxrn(klon,klev) - REAL fluxpb(klon,klev) - REAL pbimpa(klon,klev) - REAL pbnucl(klon,klev) - REAL rn(klon,klev) - REAL pb(klon,klev) - REAL flestottr(klon,klev,nbtr) ! flux de lessivage -c ! dans chaque couche - real zmasse(klon,klev) - real ztra_th(klon,klev) - -C - character*20 modname - character*80 abort_message -c -c Controles -c------------- - logical first,couchelimite,convection,lessivage,sorties, - s rnpb,inirnpb - save first,couchelimite,convection,lessivage, - s sorties,inirnpb -c$OMP THREADPRIVATE(first,couchelimite,convection,lessivage, -c$OMP+ sorties,inirnpb) -c data first,couchelimite,convection,lessivage,sorties -c s /.true.,.true.,.false.,.true.,.true./ -c Olivia - data first,couchelimite,convection,lessivage, - s sorties - s /.true.,.true.,.true.,.true.,.true./ - -! Variables needed for configuration with INCA - INTEGER :: lastgas - INTEGER :: ncsec - - REAL, PARAMETER :: dry_mass = 28.966 - REAL, POINTER :: hbuf(:) - REAL, ALLOCATABLE :: obuf(:) - REAL :: calday - REAL :: pdel(klon,klev) -c -c====================================================================== - - modname='phytrac' - - ps(:)=paprs(:,1) - - if (debutphy) then - allocate( trs(klon,nbtr) ) - allocate( masktr(klon,nbtr)) - allocate( fshtr(klon,nbtr) ) - allocate( hsoltr(nbtr)) - allocate( tautr(nbtr)) - allocate( vdeptr(nbtr)) - allocate( scavtr(nbtr)) - allocate( aerosol(nbtr)) - allocate( clsol(nbtr)) - allocate( radio(nbtr)) - - -! FH 2008/05/09 correction de la frequence d'ecriture des traceurs -! ecrit_tra = FLOAT(NINT(86400./pdtphys *ecritphy)) - print*,'dans phytrac ',pdtphys,ecrit_tra - - inirnpb=rnpb - PRINT*, 'La frequence de sortie traceurs est ', ecrit_tra -C -c============================================================= -c Initialisation des sorties -c============================================================= - -#ifdef CPP_IOIPSL -#include "ini_histrac.h" -#endif - -c====================================================================== -c Initialisation de certaines variables pour le Rn et le Pb -c====================================================================== - -c Initialisation du traceur dans le sol (couche limite radonique) -c -c print*,'valeur de debut dans phytrac :',debutphy - trs(:,:) = 0. -c$OMP MASTER - if (is_mpi_root) then - trs_tmp(:)=0. - open (99,file='starttrac',status='old', - . err=999,form='formatted') - read(99,*) (trs_tmp(i),i=1,klon_glo) -999 close(99) - endif -c$OMP END MASTER - call Scatter(trs_tmp,trs(:,1)) - -c print*, 'apres starttrac' - -c Initialisation de la fraction d'aerosols lessivee -c - d_tr_lessi_impa(:,:,:) = 0. - d_tr_lessi_nucl(:,:,:) = 0. -c -c Initialisation de la nature des traceurs -c - DO it = 1, nbtr - aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut - radio(it) = .FALSE. ! Par defaut pas de passage par radiornpb - clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit - ENDDO -c - ENDIF ! fin debutphy -c Initialisation du traceur dans le sol (couche limite radonique) - if(inirnpb) THEN -c - radio(1)= .true. - radio(2)= .true. - clsol(1)= .true. - clsol(2)= .true. - aerosol(2) = .TRUE. ! le Pb est un aerosol -c - call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr - . ,vdeptr,scavtr) - inirnpb=.false. - endif - - - IF (config_inca == 'none') THEN - DO i=1,nlon - pftsol1(i) = ftsol(i,1) - pftsol2(i) = ftsol(i,2) - pftsol3(i) = ftsol(i,3) - pftsol4(i) = ftsol(i,4) - - ppsrf1(i) = pctsrf(i,1) - ppsrf2(i) = pctsrf(i,2) - ppsrf3(i) = pctsrf(i,3) - ppsrf4(i) = pctsrf(i,4) - - ENDDO - - ELSE ! config_inca /=none -#ifdef INCA - call VTe(VTphysiq) - call VTb(VTinca) - -!====================================================================== -! Chimie -!====================================================================== - - calday = FLOAT(julien) + gmtime - ncsec = NINT (86400.*gmtime) - - DO k = 1, nlev - pdel(:,k) = paprs(:,k) - paprs (:,k+1) - END DO - - IF (config_inca == 'aero') THEN - CALL aerosolmain (aerosol_couple, - $ tr_seri, - $ pdtphys, - $ pplay, - $ pdel, - $ prfl, - $ pmflxr, - $ psfl, - $ pmflxs, - $ pmfu, - $ itop_con, - $ ibas_con, - $ pphi, - $ airephy, ! paire, - $ nstep, - $ rneb, ! for chimiaq - $ t_seri, ! for chimiaq - $ rh, - $ tau_aero, - $ piz_aero, - $ cg_aero, - $ rfname, - $ ccm, - $ lafin) - END IF - - CALL chemmain (tr_seri, !mmr - $ nstep, !nstep - $ calday, !calday - $ julien, !ncdate - $ ncsec, !ncsec - $ 1, !lat - $ pdtphys, !delt - $ paprs(1,1), !ps - $ pplay, !pmid - $ pdel, !pdel - $ airephy, - $ pctsrf(1,1),!oro - $ ftsol, !tsurf - $ albsol, !albs - $ pphi, !zma - $ pphis, !phis - $ cldfra, !cldfr - $ rneb, !cldfr_st - $ diafra, !cldfr_cv - $ itop_con, !cldtop - $ ibas_con, !cldbot - $ cldliq, !cwat - $ prfl, !flxrst - $ pmflxr, !flxrcv - $ psfl, !flxsst - $ pmflxs, !flxscv - $ pmfu, !flxupd - $ flxmass_w, !flxmass_w - $ t_seri, !tfld - $ sh, !sh - $ rh, !rh - $ iip1, !nx - $ jjp1, !ny - $ source, - $ solsym) - - - call VTe(VTinca) - call VTb(VTphysiq) -#endif - END IF ! config_inca -c====================================================================== -c Calcul de l'effet de la convection -c====================================================================== -c print*,'Avant convection' - do it=1,nbtr - WRITE(itn,'(i2)') it -c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn) - enddo - - if (convection) then - -c print*,'Pas de temps dans phytrac : ',pdtphys - DO it=1, nbtr - - IF ( config_inca/='none') THEN - IF ( conv_flg(it) == 0 ) CYCLE - END IF - - if (iflag_con.lt.2) then - d_tr_cv=0. - else if (iflag_con.eq.2) then -c tiedke - CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, - . pplay, paprs, tr_seri(1,1,it), d_tr_cv(1,1,it)) - else -c KE - call cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(1,1,it), - . upwd,dnwd,d_tr_cv(1,1,it)) - endif - - DO k = 1, nlev - DO i = 1, klon - tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it) - ENDDO - ENDDO - - IF (config_inca == 'none') THEN - CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn) - ELSE - CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = ' - . //solsym(it)) - END IF - - ENDDO - - endif ! convection -c print*,'Apres convection' -c do it=1,nbtr -c WRITE(itn,'(i1)') it -c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn) -c enddo - - -c====================================================================== -c Calcul de l'effet des thermiques -c====================================================================== - - do k=1,klev - do i=1,klon - zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg - enddo - enddo - -c print*,'masse dans ph ',zmasse - do it=1,nbtr - do k=1,klev - do i=1,klon - d_tr_th(i,k,it)=0. - tr_seri(i,k,it)=max(tr_seri(i,k,it),0.) - tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10) - enddo - enddo - enddo - - if (iflag_thermals.gt.0) then -c print*,'calcul de leffet des thermiques' - nsplit=10 - DO it=1, nbtr -c WRITE(itn,'(i1)') it -c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'conv it='//itn) -c print*,'avant dqthermiquesretro' -c call dump2d(iim,jjm-1,tr_seri(2,1,1),'TR_SERI ') - - do isplit=1,nsplit -c Abderr 25 11 02 -C Thermiques -c print*,'Avant dans phytrac' - call dqthermcell(klon,klev,pdtphys/nsplit - . ,fm_therm,entr_therm,zmasse - . ,tr_seri(1:klon,1:klev,it),d_tr,ztra_th) - - do k=1,klev - do i=1,klon - d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit - d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k) - tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.) - enddo - enddo - enddo ! nsplit -c print*,'apres thermiques' -c call dump2d(iim,jjm-1,d_tr_th(1,1,1),'d_tr_th ') -c do k=1,klev -c print*,'d_tr_th(',k,')=',tr_seri(280,k,1) -c enddo - -c WRITE(itn,'(i1)') it -c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'therm it='//itn) - ENDDO ! it - endif ! Thermiques -c print*,'ATTENTION: sdans thermniques' - -c====================================================================== -c Calcul de l'effet de la couche limite -c====================================================================== -c print *,'Avant couchelimite' -c do it=1,nbtr -c WRITE(itn,'(i1)') it -c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn) -c enddo - - if (couchelimite) then - - DO k = 1, nlev - DO i = 1, klon - delp(i,k) = paprs(i,k)-paprs(i,k+1) - ENDDO - ENDDO - -C maf modif pour tenir compte du cas rnpb + traceur - DO it=1, nbtr - - IF ( config_inca/='none' ) THEN - IF( pbl_flg(it) == 0 ) CYCLE - END IF - -c print *,'it',it,clsol(it) - if (clsol(it)) then ! couche limite avec quantite dans le sol calculee - CALL cltracrn(it, pdtphys, yu1, yv1, - e cdragh, coefh,t_seri,ftsol,pctsrf, - e tr_seri(1,1,it),trs(1,it), - e paprs, pplay, delp, - e masktr(1,it),fshtr(1,it),hsoltr(it), - e tautr(it),vdeptr(it), - e xlat, - s d_tr_cl(1,1,it),d_trs) - DO k = 1, nlev - DO i = 1, klon - tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it) - ENDDO - ENDDO -c -c Traceur ds sol -c - DO i = 1, klon - trs(i,it) = trs(i,it) + d_trs(i) - END DO -C -C maf provisoire suppression des prints -C WRITE(itn,'(i1)') it -C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn) - else ! couche limite avec flux prescrit - - IF (config_inca == 'none') THEN -Cmaf provisoire source / traceur a creer - DO i=1, klon - source(i,it) = 0.0 ! pas de source, pour l'instant - ENDDO - END IF - - CALL cltrac(pdtphys, coefh,t_seri, - s tr_seri(1,1,it), source(:,it), - e paprs, pplay, delp, - s d_tr_cl(1,1,it)) - DO k = 1, nlev - DO i = 1, klon - tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it) - ENDDO - ENDDO -Cmaf WRITE(itn,'(i1)') it -cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn) - endif - ENDDO -c - endif ! couche limite - -c print*,'Apres couchelimite' -c do it=1,nbtr -c WRITE(itn,'(i1)') it -c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn) -c enddo - -c====================================================================== -c Calcul de l'effet du puits radioactif -c====================================================================== - -C MAF il faudrait faire une modification pour passer dans radiornpb -c si radio=true mais pour l'instant radiornpb propre au cas rnpb - if(rnpb) then -c print *, 'decroissance radiactive activee' - call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec) -C - DO it=1,nbtr - if(radio(it)) then - DO k = 1, nlev - DO i = 1, klon - tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it) - ENDDO - ENDDO - WRITE(itn,'(i1)') it - CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn) - endif - ENDDO -c - endif ! rnpb decroissance radioactive -C -c====================================================================== -c Calcul de l'effet de la precipitation -c====================================================================== - -c print*,'LESSIVAGE =',lessivage - IF (lessivage) THEN - -c print*,'avant lessivage' - - d_tr_lessi_nucl(:,:,:) = 0. - d_tr_lessi_impa(:,:,:) = 0. - flestottr(:,:,:) = 0. -c -c tendance des aerosols nuclees et impactes -c - DO it = 1, nbtr - IF (aerosol(it)) THEN - DO k = 1, nlev - DO i = 1, klon - d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) + - s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it) - d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) + - s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it) - ENDDO - ENDDO - ENDIF - ENDDO -c -c Mises a jour des traceurs + calcul des flux de lessivage -c Mise a jour due a l'impaction et a la nucleation -c -c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA') -c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL') -c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3') - DO it = 1, nbtr -c print*,'IT=',it,aerosol(it) - IF (aerosol(it)) THEN -c print*,'IT=',it,' On lessive' - DO k = 1, nlev - DO i = 1, klon - tr_seri(i,k,it)=tr_seri(i,k,it) - s *frac_impa(i,k)*frac_nucl(i,k) - ENDDO - ENDDO - ENDIF - ENDDO -c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B') -c -c Flux lessivage total -c - DO it = 1, nbtr - DO k = 1, nlev - DO i = 1, klon - flestottr(i,k,it) = flestottr(i,k,it) - - s ( d_tr_lessi_nucl(i,k,it) + - s d_tr_lessi_impa(i,k,it) ) * - s ( paprs(i,k)-paprs(i,k+1) ) / - s (RG * pdtphys) - ENDDO - ENDDO -c -Cmaf WRITE(itn,'(i1)') it -Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn) - ENDDO -c -c print*,'apres lessivage' - ENDIF -Cc - DO k = 1, nlev - DO i = 1, klon - fluxrn(i,k) = flestottr(i,k,1) - fluxpb(i,k) = flestottr(i,k,2) - rn(i,k) = tr_seri(i,k,1) - pb(i,k) = tr_seri(i,k,2) - pbnucl(i,k)=d_tr_lessi_nucl(i,k,2) - pbimpa(i,k)=d_tr_lessi_impa(i,k,2) - ENDDO - ENDDO - -c============================================================= -c Ecriture des sorties -c============================================================= - -#ifdef CPP_IOIPSL -#include "write_histrac.h" -#endif - -c============================================================= - - if (lafin) then - print*, 'c est la fin de la physique' - call Gather(trs(:,1),trs_tmp) -c$OMP MASTER - if (is_mpi_root) then - - open (99,file='restarttrac', form='formatted') - do i=1,klon_glo - write(99,*) trs_tmp(i) - enddo - PRINT*, 'Ecriture du fichier restarttrac' - close(99) - endif -c$OMP END MASTER - else -c print*, 'physique pas fini' - endif - - - RETURN - END Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F90 (revision 1191) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F90 (revision 1191) @@ -0,0 +1,410 @@ +!$Id $ + +SUBROUTINE phytrac( & + nstep, julien, gmtime, debutphy, & + lafin, pdtphys, u, v, t_seri, & + paprs, pplay, pmfu, pmfd, & + pen_u, pde_u, pen_d, pde_d, & + cdragh, coefh, fm_therm, entr_therm,& + yu1, yv1, ftsol, pctsrf, & + xlat, frac_impa,frac_nucl,xlon, & + presnivs, pphis, pphi, albsol, & + sh, rh, cldfra, rneb, & + diafra, cldliq, itop_con, ibas_con, & + pmflxr, pmflxs, prfl, psfl, & + da, phi, mp, upwd, & + dnwd, aerosol_couple, flxmass_w, & + tau_aero, piz_aero, cg_aero, ccm, & + rfname, & + tr_seri) +! +!====================================================================== +! Auteur(s) FH +! Objet: Moniteur general des tendances traceurs +!====================================================================== + + USE ioipsl + USE dimphy + USE infotrac + USE mod_grid_phy_lmdz + USE mod_phys_lmdz_para + USE comgeomphy + USE iophy + USE traclmdz_mod + USE tracinca_mod + + + IMPLICIT NONE + + INCLUDE "YOMCST.h" + INCLUDE "dimensions.h" + INCLUDE "indicesol.h" + INCLUDE "clesphys.h" + INCLUDE "temps.h" + INCLUDE "paramet.h" + INCLUDE "control.h" + INCLUDE "thermcell.h" +!========================================================================== +! -- ARGUMENT DESCRIPTION -- +!========================================================================== + +! Input arguments +!---------------- +!Configuration grille,temps: + INTEGER,INTENT(IN) :: nstep ! Appel physique + INTEGER,INTENT(IN) :: julien ! Jour julien + REAL,INTENT(IN) :: gmtime + REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde) + LOGICAL,INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique + LOGICAL,INTENT(IN) :: lafin ! le flag de la fin de la physique + + REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point + REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point +! +!Physique: +!-------- + REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature + REAL,DIMENSION(klon,klev),INTENT(IN) :: u ! + REAL,DIMENSION(klon,klev),INTENT(IN) :: v ! + REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique + REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) + REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) + REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentiel + REAL,DIMENSION(klon),INTENT(IN) :: pphis + REAL,DIMENSION(klev),INTENT(IN) :: presnivs + REAL,DIMENSION(klon,klev),INTENT(IN) :: cldliq ! eau liquide nuageuse + REAL,DIMENSION(klon,klev),INTENT(IN) :: cldfra ! fraction nuageuse (tous les nuages) + REAL,DIMENSION(klon,klev),INTENT(IN) :: diafra ! fraction nuageuse (convection ou stratus artificiels) + REAL,DIMENSION(klon,klev),INTENT(IN) :: rneb ! fraction nuageuse (grande echelle) + INTEGER,DIMENSION(klon),INTENT(IN) :: itop_con + INTEGER,DIMENSION(klon),INTENT(IN) :: ibas_con + REAL,DIMENSION(klon),INTENT(IN) :: albsol ! albedo surface +! +!Convection: +!---------- + REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfu ! flux de masse dans le panache montant + REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfd ! flux de masse dans le panache descendant + REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_u ! flux entraine dans le panache montant + REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_u ! flux detraine dans le panache montant + REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_d ! flux entraine dans le panache descendant + REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_d ! flux detraine dans le panache descendant + +!...Tiedke + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: pmflxr, pmflxs ! Flux precipitant de pluie, neige aux interfaces [convection] + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: prfl, psfl ! Flux precipitant de pluie, neige aux interfaces [large-scale] + + LOGICAL,INTENT(IN) :: aerosol_couple + REAL,DIMENSION(klon,klev),INTENT(IN) :: flxmass_w + REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: tau_aero + REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: piz_aero + REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: cg_aero + CHARACTER(len=4),DIMENSION(9),INTENT(IN) :: rfname + REAL,DIMENSION(klon,klev,2),INTENT(IN) :: ccm +!... K.Emanuel + REAL,DIMENSION(klon,klev),INTENT(IN) :: da + REAL,DIMENSION(klon,klev,klev),INTENT(IN):: phi + REAL,DIMENSION(klon,klev),INTENT(IN) :: mp + REAL,DIMENSION(klon,klev),INTENT(IN) :: upwd ! saturated updraft mass flux + REAL,DIMENSION(klon,klev),INTENT(IN) :: dnwd ! saturated downdraft mass flux +! +!Thermiques: +!---------- + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: fm_therm + REAL,DIMENSION(klon,klev),INTENT(IN) :: entr_therm +! +!Couche limite: +!-------------- +! + REAL,DIMENSION(klon,klev),INTENT(IN) :: cdragh ! coeff drag pour T et Q + REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh ! coeff melange CL (m**2/s) + REAL,DIMENSION(klon),INTENT(IN) :: yu1 ! vents au premier niveau + REAL,DIMENSION(klon),INTENT(IN) :: yv1 ! vents au premier niveau +! +!Lessivage: +!---------- +! +! pour le ON-LINE +! + REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_impa ! fraction d'aerosols non impactes + REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_nucl ! fraction d'aerosols non nuclees + +! Arguments necessaires pour les sources et puits de traceur: + REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol ! Temperature du sol (surf)(Kelvin) + REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol (nature du sol) + + +! Output argument +!---------------- + REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] + +!======================================================================================= +! -- LOCAL VARIABLES -- +!======================================================================================= + + INTEGER :: i, k, it + INTEGER :: nsplit + +!Sources et Reservoirs de traceurs (ex:Radon): +!-------------------------------------------- +! + REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: source ! a voir lorsque le flux de surface est prescrit +!$OMP THREADPRIVATE(source) + +! +!Entrees/Sorties: (cf ini_histrac.h et write_histrac.h) +!--------------- + INTEGER :: iiq, ierr + INTEGER :: nhori, nvert + REAL :: zsto, zout, zjulian + INTEGER,SAVE :: nid_tra ! pointe vers le fichier histrac.nc +!$OMP THREADPRIVATE(nid_tra) + REAL,DIMENSION(klon) :: zx_tmp_fi2d ! variable temporaire grille physique + INTEGER :: itau_w ! pas de temps ecriture = nstep + itau_phy + LOGICAL,PARAMETER :: ok_sync=.TRUE. + +! +! Nature du traceur +!------------------ + LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: aerosol ! aerosol(it) = true => aerosol => lessivage +!$OMP THREADPRIVATE(aerosol) ! aerosol(it) = false => gaz + REAL,DIMENSION(klon,klev) :: delp ! epaisseur de couche (Pa) +! +! Tendances de traceurs (Td): +!------------------------ +! + REAL,DIMENSION(klon,klev) :: d_tr ! Td dans l'atmosphere + REAL,DIMENSION(klon,klev,nbtr) :: d_tr_cl ! Td couche limite/traceur + REAL,DIMENSION(klon,klev,nbtr) :: d_tr_cv ! Td convection/traceur + REAL,DIMENSION(klon,klev,nbtr) :: d_tr_th ! Td thermique + REAL,DIMENSION(klon,klev,nbtr) :: d_tr_lessi_impa ! Td du lessivage par impaction + REAL,DIMENSION(klon,klev,nbtr) :: d_tr_lessi_nucl ! Td du lessivage par nucleation +! +! Physique +!---------- + REAL,DIMENSION(klon,klev,nbtr) :: flestottr ! flux de lessivage dans chaque couche + REAL,DIMENSION(klon,klev) :: zmasse ! densité atmosphérique Kg/m2 + REAL,DIMENSION(klon,klev) :: ztra_th + +!Controles: +!--------- + LOGICAL,SAVE :: couchelimite=.TRUE. + LOGICAL,SAVE :: convection=.TRUE. + LOGICAL,SAVE :: lessivage +!$OMP THREADPRIVATE(couchelimite,convection,lessivage) + + CHARACTER(len=8),DIMENSION(nbtr) :: solsym + + +!###################################################################### +! -- INITIALIZATION -- +!###################################################################### + IF (debutphy) THEN + WRITE(*,*) 'FIRST TIME IN PHYTRAC : pdtphys(sec) = ',pdtphys,'ecrit_tra (sec) = ',ecrit_tra + ALLOCATE( source(klon,nbtr), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 1',1) + + ALLOCATE( aerosol(nbtr), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 2',1) + + + ! Initialize module for specific tracers + SELECT CASE(type_trac) + CASE('lmdz') + CALL traclmdz_init(pctsrf, ftsol, aerosol, lessivage) + CASE('inca') + source(:,:)=0. + CALL tracinca_init(aerosol,lessivage) + END SELECT +! +! Initialize diagnostic output +! ---------------------------- + INCLUDE "ini_histrac.h" + + END IF +!############################################ END INITIALIZATION ####### + +!=============================================================================== +! -- Do specific treatment according to chemestry model or local LMDZ tracers +! +!=============================================================================== + SELECT CASE(type_trac) + CASE('lmdz') + ! -- Traitement des traceurs avec traclmdz + + CALL traclmdz(& + pdtphys, t_seri, & + paprs, pplay, cdragh, coefh, & + yu1, yv1, ftsol, pctsrf, & + xlat, couchelimite, & + tr_seri, source, solsym, d_tr_cl) + + CASE('inca') + ! -- CHIMIE INCA config_inca = aero or chem -- + + CALL tracinca(& + nstep, julien, gmtime, lafin, & + pdtphys, t_seri, paprs, pplay, & + pmfu, ftsol, pctsrf, pphis, & + pphi, albsol, sh, rh, & + cldfra, rneb, diafra, cldliq, & + itop_con, ibas_con, pmflxr, pmflxs, & + prfl, psfl, aerosol_couple, flxmass_w, & + tau_aero, piz_aero, cg_aero, ccm, & + rfname, & + tr_seri, source, solsym) + END SELECT + +!====================================================================== +! -- Calcul de l'effet de la convection -- +!====================================================================== + IF (convection) THEN + DO it=1, nbtr + IF ( conv_flg(it) == 0 ) CYCLE + + IF (iflag_con.LT.2) THEN + d_tr_cv(:,:,:)=0. + ELSE IF (iflag_con.EQ.2) THEN +!..Tiedke + CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, & + pplay, paprs, tr_seri(:,:,it), d_tr_cv(:,:,it)) + ELSE +!..K.Emanuel + CALL cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(:,:,it),& + upwd,dnwd,d_tr_cv(:,:,it)) + END IF + + DO k = 1, klev + DO i = 1, klon + tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it) + END DO + END DO + + CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'convection it = '//solsym(it)) + + END DO ! nbtr + END IF ! convection + +!====================================================================== +! -- Calcul de l'effet des thermiques -- +!====================================================================== + + DO k=1,klev + DO i=1,klon + zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg + END DO + END DO + + DO it=1,nbtr + DO k=1,klev + DO i=1,klon + d_tr_th(i,k,it)=0. + tr_seri(i,k,it)=MAX(tr_seri(i,k,it),0.) + tr_seri(i,k,it)=MIN(tr_seri(i,k,it),1.e10) + END DO + END DO + END DO + + IF (iflag_thermals.GT.0) THEN + nsplit=10 + DO it=1, nbtr + DO isplit=1,nsplit + + CALL dqthermcell(klon,klev,pdtphys/nsplit, & + fm_therm,entr_therm,zmasse, & + tr_seri(1:klon,1:klev,it),d_tr,ztra_th) + + DO k=1,klev + DO i=1,klon + d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit + d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k) + tr_seri(i,k,it)=MAX(tr_seri(i,k,it)+d_tr(i,k),0.) + END DO + END DO + END DO ! nsplit + END DO ! it + END IF ! Thermiques + +!====================================================================== +! -- Calcul de l'effet de la couche limite -- +!====================================================================== + + IF (couchelimite) THEN + + DO k = 1, klev + DO i = 1, klon + delp(i,k) = paprs(i,k)-paprs(i,k+1) + END DO + END DO + + DO it=1, nbtr + + IF( pbl_flg(it) /= 0 ) THEN + + CALL cltrac(pdtphys, coefh,t_seri, & + tr_seri(:,:,it), source(:,it), & + paprs, pplay, delp, & + d_tr_cl(:,:,it)) + + DO k = 1, klev + DO i = 1, klon + tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it) + END DO + END DO + END IF + + END DO + + END IF ! couche limite + + +!====================================================================== +! Calcul de l'effet de la precipitation +!====================================================================== + + IF (lessivage) THEN + + d_tr_lessi_nucl(:,:,:) = 0. + d_tr_lessi_impa(:,:,:) = 0. + flestottr(:,:,:) = 0. +!========================= +! LESSIVAGE LARGE SCALE : +!========================= + +! Tendance des aerosols nuclees et impactes +! ----------------------------------------- + DO it = 1, nbtr + IF (aerosol(it)) THEN + DO k = 1, klev + DO i = 1, klon + d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) + & + ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it) + d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) + & + ( 1 - frac_impa(i,k) )*tr_seri(i,k,it) + +! +! Flux lessivage total +! ------------------------------------------------------------ + flestottr(i,k,it) = flestottr(i,k,it) - & + ( d_tr_lessi_nucl(i,k,it) + & + d_tr_lessi_impa(i,k,it) ) * & + ( paprs(i,k)-paprs(i,k+1) ) / & + (RG * pdtphys) +! +! Mise a jour des traceurs due a l'impaction,nucleation +! ---------------------------------------------------------------------- + tr_seri(i,k,it)=tr_seri(i,k,it)*frac_impa(i,k)*frac_nucl(i,k) + END DO + END DO + END IF + END DO + + END IF ! lessivage + +!============================================================= +! Ecriture des sorties +!============================================================= + INCLUDE "write_histrac.h" + + +END SUBROUTINE phytrac Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/radio_decay.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/radio_decay.F90 (revision 1191) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/radio_decay.F90 (revision 1191) @@ -0,0 +1,61 @@ +! +! $Id $ +! +SUBROUTINE radio_decay(radio,rnpb,dtime,tautr,tr,d_tr) +! +! Caluclate radioactive decay for all tracers with radio(it)=true +! + USE dimphy + USE infotrac, ONLY : nbtr + 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. +!----------------------------------------------------------------------- +! +! Entrees +! + LOGICAL,DIMENSION(nbtr),INTENT(IN) :: radio ! .true. = traceur radioactif + LOGICAL,INTENT(IN) :: rnpb ! .true. = decroissance RN = source PB + REAL,INTENT(IN) :: dtime ! Pas de temps physique (secondes) + REAL,DIMENSION(nbtr),INTENT(IN) :: tautr ! Constante de decroissance radioactive + REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr ! Concentrations traceurs U/kgA +! +! Sortie +! + REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: d_tr ! Tendance de decroissance radioactive +! +! Locales +! + INTEGER :: i,k,it + + + DO it = 1,nbtr + IF ( radio(it) ) THEN + IF (tautr(it) .GT. 0.) THEN + DO k = 1,klev + DO i = 1,klon + d_tr(i,k,it) = - tr(i,k,it) * dtime / tautr(it) + 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] +!------------------------------------------------------- + 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) + ENDDO + ENDDO + ENDIF + +END SUBROUTINE radio_decay Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/radiornpb.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/radiornpb.F (revision 1190) +++ (revision ) @@ -1,55 +1,0 @@ -! -! $Header$ -! - SUBROUTINE radiornpb(tr,dtime,tautr,d_tr) - USE dimphy - USE infotrac, ONLY : nbtr - IMPLICIT none -c====================================================================== -c Auteur(s): AA + CG (LGGE/CNRS) Date 24-06-94 -c Objet: Decroissance radioactive d'un traceur dans l'atmosphere -CG240694 : Pour un traceur, le radon -CG161294 : Plus un 2eme traceur, le 210Pb. Le radon decroit en plomb. -c====================================================================== -c Arguments: -c====================================================================== -cym#include "dimensions.h" -cym#include "dimphy.h" -c====================================================================== -C - INTEGER i , k , it - REAL tr(klon,klev,nbtr) , d_tr(klon,klev,nbtr) - REAL dtime - REAL tautr(nbtr) -C -c WRITE(*,'(''PASSAGE radiornpb ... '',$)') -C Attention, pour un pas de temps beaucoup plus petit que la decroissance!!! - - DO it = 1,2 - IF ( tautr(it) .GT. 0. ) THEN - DO k = 1,klev - DO i = 1,klon - d_tr(i,k,it) = - tr(i,k,it) * dtime / tautr(it) - END DO - END DO - ELSE - DO k = 1,klev - DO i = 1,klon - d_tr(i,k,it) = 0. - END DO - END DO - END IF - END DO -C -CG161294 : Cas particulier radon 1 => plomb 2 -c - DO k = 1,klev - DO i = 1,klon - d_tr(i,k,2) = d_tr(i,k,2) - d_tr(i,k,1) - ENDDO - ENDDO -c -c WRITE(*,*) ' radiornpb OK' -c - RETURN - END Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/tracinca_mod.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/tracinca_mod.F90 (revision 1191) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/tracinca_mod.F90 (revision 1191) @@ -0,0 +1,191 @@ +!$Id $ +! +MODULE tracinca_mod +! +! This module prepares and calls the INCA main subroutines. +! + +CONTAINS + + SUBROUTINE tracinca_init(aerosol,lessivage) + ! This subroutine initialize some control varaibles. + + USE infotrac + IMPLICIT NONE + + ! Output variables + LOGICAL,DIMENSION(nbtr), INTENT(OUT) :: aerosol + LOGICAL,INTENT(OUT) :: lessivage + + + ! Initialization + lessivage =.FALSE. + aerosol(:) = .FALSE. + + END SUBROUTINE tracinca_init + + SUBROUTINE tracinca( & + nstep, julien, gmtime, lafin, & + pdtphys, t_seri, paprs, pplay, & + pmfu, ftsol, pctsrf, pphis, & + pphi, albsol, sh, rh, & + cldfra, rneb, diafra, cldliq, & + itop_con, ibas_con, pmflxr, pmflxs, & + prfl, psfl, aerosol_couple, flxmass_w, & + tau_aero, piz_aero, cg_aero, ccm, & + rfname, & + tr_seri, source, solsym) + +!======================================================== +! -- CHIMIE INCA -- +!======================================================== + + USE dimphy + USE infotrac + USE vampir + USE comgeomphy + + IMPLICIT NONE + + INCLUDE "indicesol.h" + INCLUDE "control.h" + INCLUDE "dimensions.h" + INCLUDE "paramet.h" + +!========================================================================== +! -- DESCRIPTION DES ARGUMENTS -- +!========================================================================== + + +! EN ENTREE ... +! +!Configuration grille,temps: + INTEGER,INTENT(IN) :: nstep ! Appel physique + INTEGER,INTENT(IN) :: julien ! Jour julien + REAL,INTENT(IN) :: gmtime + REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde) + LOGICAL,INTENT(IN) :: lafin ! le flag de la fin de la physique + + +!Physique: +!-------- + REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature + REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique + REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) + REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) + REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentiel + REAL,DIMENSION(klon),INTENT(IN) :: pphis + REAL,DIMENSION(klon,klev),INTENT(IN) :: cldliq ! eau liquide nuageuse + REAL,DIMENSION(klon,klev),INTENT(IN) :: cldfra ! fraction nuageuse (tous les nuages) + REAL,DIMENSION(klon,klev),INTENT(IN) :: diafra ! fraction nuageuse (convection ou stratus artificiels) + REAL,DIMENSION(klon,klev),INTENT(IN) :: rneb ! fraction nuageuse (grande echelle) + INTEGER,DIMENSION(klon),INTENT(IN) :: itop_con + INTEGER,DIMENSION(klon),INTENT(IN) :: ibas_con + REAL,DIMENSION(klon),INTENT(IN) :: albsol ! albedo surface +! +!Convection: +!---------- + REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfu ! flux de masse dans le panache montant + +!...Tiedke + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: pmflxr, pmflxs ! Flux precipitant de pluie, neige aux interfaces [convection] + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: prfl, psfl ! Flux precipitant de pluie, neige aux interfaces [large-scale] + + LOGICAL,INTENT(IN) :: aerosol_couple + REAL,DIMENSION(klon,klev),INTENT(IN) :: flxmass_w + REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: tau_aero + REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: piz_aero + REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: cg_aero + CHARACTER(len=4),DIMENSION(9),INTENT(IN) :: rfname + REAL,DIMENSION(klon,klev,2),INTENT(IN) :: ccm + +! Arguments necessaires pour les sources et puits de traceur: + REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol ! Temperature du sol (surf)(Kelvin) + REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol) + + + ! InOutput argument + REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] + + ! Output arguments + REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: source ! a voir lorsque le flux de surface est prescrit + CHARACTER(len=8),DIMENSION(nbtr), INTENT(OUT) :: solsym + +!======================================================================================= +! -- VARIABLES LOCALES TRACEURS -- +!======================================================================================= + + INTEGER :: k + REAL,DIMENSION(klon,klev) :: pdel + REAL :: calday + INTEGER :: ncsec + + CALL VTe(VTphysiq) + CALL VTb(VTinca) + + calday = FLOAT(julien) + gmtime + ncsec = NINT (86400.*gmtime) + + DO k = 1, klev + pdel(:,k) = paprs(:,k) - paprs (:,k+1) + END DO + + IF (config_inca == 'aero') THEN +#ifdef INCA + CALL aerosolmain( & + aerosol_couple,tr_seri,pdtphys, & + pplay,pdel,prfl,pmflxr,psfl, & + pmflxs,pmfu,itop_con,ibas_con, & + pphi,airephy,nstep,rneb,t_seri, & + rh,tau_aero,piz_aero,cg_aero, & + rfname,ccm,lafin) +#endif + END IF + +#ifdef INCA + CALL chemmain (tr_seri, & !mmr + nstep, & !nstep + calday, & !calday + julien, & !ncdate + ncsec, & !ncsec + 1, & !lat + pdtphys, & !delt + paprs(1,1), & !ps + pplay, & !pmid + pdel, & !pdel + airephy, & + pctsrf(1,1),& !oro + ftsol, & !tsurf + albsol, & !albs + pphi, & !zma + pphis, & !phis + cldfra, & !cldfr + rneb, & !cldfr_st + diafra, & !cldfr_cv + itop_con, & !cldtop + ibas_con, & !cldbot + cldliq, & !cwat + prfl, & !flxrst + pmflxr, & !flxrcv + psfl, & !flxsst + pmflxs, & !flxscv + pmfu, & !flxupd + flxmass_w, & !flxmass_w + t_seri, & !tfld + sh, & !sh + rh, & !rh + iip1, & !nx + jjp1, & !ny + source, & + solsym) +#endif + + CALL VTe(VTinca) + CALL VTb(VTphysiq) + + + END SUBROUTINE tracinca + + +END MODULE tracinca_mod Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/traclmdz_mod.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/traclmdz_mod.F90 (revision 1191) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/traclmdz_mod.F90 (revision 1191) @@ -0,0 +1,329 @@ +!$Id $ +! +MODULE traclmdz_mod +! +! In this module all tracers specific to LMDZ are treated. This module is used +! only if running without any other chemestry model as INCA or REPROBUS. +! + + REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: masktr ! Masque reservoir de sol traceur +!$OMP THREADPRIVATE(masktr) ! Masque de l'echange avec la surface (1 = reservoir) ou (possible >= 1 ) + REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: fshtr ! Flux surfacique dans le reservoir de sol +!$OMP THREADPRIVATE(fshtr) + REAL,DIMENSION(:),ALLOCATABLE,SAVE :: hsoltr ! Epaisseur equivalente du reservoir de sol +!$OMP THREADPRIVATE(hsoltr) +! +!Radioelements: +!-------------- +! + REAL,DIMENSION(:),ALLOCATABLE,SAVE :: tautr ! Constante de decroissance radioactive +!$OMP THREADPRIVATE(tautr) + REAL,DIMENSION(:),ALLOCATABLE,SAVE :: vdeptr ! Vitesse de depot sec dans la couche Brownienne +!$OMP THREADPRIVATE(vdeptr) + REAL,DIMENSION(:),ALLOCATABLE,SAVE :: scavtr ! Coefficient de lessivage +!$OMP THREADPRIVATE(scavtr) + REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: srcbe ! Production du beryllium7 dans l atmosphere (U/s/kgA) +!$OMP THREADPRIVATE(srcbe) + + LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: radio ! radio(it) = true => decroisssance radioactive +!$OMP THREADPRIVATE(radio) + + REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: trs ! Conc. radon ds le sol +!$OMP THREADPRIVATE(trs) + + INTEGER,SAVE :: id_be ! Activation et position du traceur Be7 [ id_be=0 -> desactive ] +!$OMP THREADPRIVATE(id_be) + + LOGICAL,SAVE :: rnpb=.TRUE. ! Presence du couple Rn222, Pb210 +!$OMP THREADPRIVATE(rnpb) + + +CONTAINS + + SUBROUTINE traclmdz_from_restart(trs_in) + ! This subroutine initialize the module saved variable trs with values from restart file (startphy.nc). + ! This subroutine is called from phyetat0 after the field trs_in has been read. + + USE dimphy + USE infotrac + IMPLICIT NONE + + ! Input argument + REAL,DIMENSION(klon,nbtr), INTENT(IN) :: trs_in + + ! Local variables + INTEGER :: ierr + + ! Allocate restart variables trs + ALLOCATE( trs(klon,nbtr), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('traclmdz_start', 'pb in allocation 1',1) + + ! Initialize trs with values read from restart file + trs(:,:) = trs_in(:,:) + + END SUBROUTINE traclmdz_from_restart + + + SUBROUTINE traclmdz_init(pctsrf, ftsol, aerosol, lessivage) + ! This subroutine allocates and initialize module variables and control variables. + USE dimphy + USE infotrac + + IMPLICIT NONE + + INCLUDE "indicesol.h" + +! Input variables + REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol) + REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol ! Temperature du sol (surf)(Kelvin) + +! Output variables + LOGICAL,DIMENSION(nbtr), INTENT(OUT) :: aerosol + LOGICAL,INTENT(OUT) :: lessivage + +! Local variables + INTEGER :: ierr, it, iiq + +! -------------------------------------------- +! Allocation +! -------------------------------------------- + + ALLOCATE( scavtr(nbtr), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 9',1) + scavtr(:)=1. + + ALLOCATE( radio(nbtr), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 11',1) + radio(:) = .false. ! Par defaut pas decroissance radioactive + + ALLOCATE( masktr(klon,nbtr), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 2',1) + + ALLOCATE( fshtr(klon,nbtr), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 3',1) + + ALLOCATE( hsoltr(nbtr), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 4',1) + + ALLOCATE( tautr(nbtr), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 5',1) + tautr(:) = 0. + + ALLOCATE( vdeptr(nbtr), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 6',1) + vdeptr(:) = 0. + + + lessivage = .TRUE. + aerosol(:) = .FALSE. ! Tous les traceurs sont des gaz par defaut + +! +! Recherche des traceurs connus : Be7, CO2,... +! -------------------------------------------- + DO it=1,nbtr + iiq=niadv(it+2) + ! + ! Recherche du Beryllium 7 + ! + IF ( tname(iiq) == "BE" .OR. tname(iiq) == "Be" .OR. & + tname(iiq) == "BE7" .OR. tname(iiq) == "Be7" ) THEN + id_be=it + ALLOCATE( srcbe(klon,klev) ) + radio(id_be) = .TRUE. + aerosol(id_be) = .TRUE. ! le Be est un aerosol + CALL init_be(pctsrf,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe) + WRITE(*,*) 'Initialisation srcBe: OK' + ELSE + id_be=0 + END IF + + ! Placer ici la recherche de nouveaux traceurs + ! ... + END DO +! +! 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 + + CALL initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr,vdeptr,scavtr) + END IF + + END SUBROUTINE traclmdz_init + + SUBROUTINE traclmdz( & + pdtphys, t_seri, & + paprs, pplay, cdragh, coefh, & + yu1, yv1, ftsol, pctsrf, & + xlat, couchelimite, & + tr_seri, source, solsym, d_tr_cl) + + USE dimphy + USE infotrac + + IMPLICIT NONE + + INCLUDE "YOMCST.h" + INCLUDE "indicesol.h" + +!========================================================================== +! -- DESCRIPTION DES ARGUMENTS -- +!========================================================================== + +! Input arguments +! +!Configuration grille,temps: + REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde) + REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point + +! +!Physique: +!-------- + REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) + REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) + + +!Couche limite: +!-------------- +! + REAL,DIMENSION(klon,klev),INTENT(IN) :: cdragh ! coeff drag pour T et Q + REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh ! coeff melange CL (m**2/s) + REAL,DIMENSION(klon),INTENT(IN) :: yu1 ! vents au premier niveau + REAL,DIMENSION(klon),INTENT(IN) :: yv1 ! vents au premier niveau + LOGICAL,INTENT(IN) :: couchelimite + +! Arguments necessaires pour les sources et puits de traceur: + REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol ! Temperature du sol (surf)(Kelvin) + REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol f(nature du sol) + +! InOutput argument + REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] + +! Output argument + CHARACTER(len=8),DIMENSION(nbtr), INTENT(OUT) :: solsym + REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: source ! a voir lorsque le flux de surface est prescrit + REAL,DIMENSION(klon,klev,nbtr), INTENT(OUT) :: d_tr_cl ! Td couche limite/traceur + + +!======================================================================================= +! -- VARIABLES LOCALES TRACEURS -- +!======================================================================================= + + INTEGER :: i, k, it + + + REAL,DIMENSION(klon) :: d_trs ! Td dans le reservoir + REAL,DIMENSION(klon,klev) :: delp ! epaisseur de couche (Pa) + + REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec ! Td radioactive + REAL :: zrho ! Masse Volumique de l'air KgA/m3 + +! +! +!================================================================= +! Ajout de la production en Be7 (Beryllium) srcbe U/s/kgA +!================================================================= +! + IF ( id_be /= 0 ) THEN + DO k = 1, klev + DO i = 1, klon + tr_seri(i,k,id_be) = tr_seri(i,k,id_be)+srcbe(i,k)*pdtphys + END DO + END DO + WRITE(*,*) 'Ajout srcBe dans tr_seri: OK' + END IF + + + DO it=1,nbtr + WRITE(solsym(it),'(i2)') it + END DO +!====================================================================== +! -- Calcul de l'effet de la couche limite -- +!====================================================================== + + IF (couchelimite) THEN + DO it=1, nbtr + IF ( it == id_be ) THEN + DO i=1, klon + zrho = pplay(i,1)/t_seri(i,1)/RD + source(i,it) = - vdeptr(it)*tr_seri(i,1,it)*zrho + END DO + ELSE + source(:,it) = 0.0 + END IF + END DO + END IF + + + DO k = 1, klev + DO i = 1, klon + delp(i,k) = paprs(i,k)-paprs(i,k+1) + END DO + END DO + + DO it=1, nbtr + IF (couchelimite .AND. pbl_flg(it) == 0 ) THEN ! couche limite avec quantite dans le sol calculee + CALL cltracrn(it, pdtphys, yu1, yv1, & + cdragh, coefh,t_seri,ftsol,pctsrf, & + tr_seri(:,:,it),trs(:,it), & + paprs, pplay, delp, & + masktr(:,it),fshtr(:,it),hsoltr(it),& + tautr(it),vdeptr(it), & + xlat,d_tr_cl(:,:,it),d_trs) + + DO k = 1, klev + DO i = 1, klon + tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it) + END DO + END DO + + ! Traceur dans le reservoir sol + DO i = 1, klon + trs(i,it) = trs(i,it) + d_trs(i) + END DO + END IF + END DO + +!====================================================================== +! Calcul de l'effet du puits radioactif +!====================================================================== + CALL radio_decay (radio,rnpb,pdtphys,tautr,tr_seri,d_tr_dec) + + DO it=1,nbtr + IF(radio(it)) then + DO k = 1, klev + DO i = 1, klon + tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it) + END DO + END DO + CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'puits rn it='//solsym(it)) + END IF + END DO + + + END SUBROUTINE traclmdz + + + SUBROUTINE traclmdz_to_restart(trs_out) + ! This subroutine is called from phyredem.F where the module + ! variable trs is written to restart file (restartphy.nc) + USE dimphy + USE infotrac + + IMPLICIT NONE + + REAL,DIMENSION(klon,nbtr), INTENT(OUT) :: trs_out + + trs_out(:,:) = trs(:,:) + + END SUBROUTINE traclmdz_to_restart + + +END MODULE traclmdz_mod Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/write_histrac.h =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/write_histrac.h (revision 1190) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/write_histrac.h (revision 1191) @@ -1,83 +1,85 @@ +!$Id $ +!*************************************** +! ECRITURE DU FICHIER : histrac.nc +!*************************************** +#ifdef CPP_IOIPSL + IF (ecrit_tra > 0. .AND. config_inca == 'none') THEN + + itau_w = itau_phy + nstep + + CALL histwrite_phy(nid_tra,"phis",itau_w,pphis) + CALL histwrite_phy(nid_tra,"aire",itau_w,airephy) + +!TRACEURS +!---------------- + DO it=1,nbtr + iiq=niadv(it+2) + +! CONCENTRATIONS + CALL histwrite_phy(nid_tra,tname(iiq),itau_w,tr_seri(:,:,it)) + +! TD LESSIVAGE + IF (lessivage .AND. aerosol(it)) THEN + CALL histwrite_phy(nid_tra,"fl"//tname(iiq),itau_w,flestottr(:,:,it)) + ENDIF + +! TD THERMIQUES + IF (iflag_thermals.gt.0) THEN + CALL histwrite_phy(nid_tra,"d_tr_th_"//tname(iiq),itau_w,d_tr_th(:,:,it)) + ENDIF + +! TD CONVECTION + IF (iflag_con.GE.2) THEN + CALL histwrite_phy(nid_tra,"d_tr_cv_"//tname(iiq),itau_w,d_tr_cv(:,:,it)) + ENDIF + +! TD COUCHE-LIMITE + CALL histwrite_phy(nid_tra,"d_tr_cl_"//tname(iiq),itau_w,d_tr_cl(:,:,it)) + ENDDO +!--------------- ! -! $Header$ ! +! VENT (niveau 1) + CALL histwrite_phy(nid_tra,"pyu1",itau_w,yu1) + CALL histwrite_phy(nid_tra,"pyv1",itau_w,yv1) +! +! TEMPERATURE DU SOL + zx_tmp_fi2d(:)=ftsol(:,1) + CALL histwrite_phy(nid_tra,"ftsol1",itau_w,zx_tmp_fi2d) + zx_tmp_fi2d(:)=ftsol(:,2) + CALL histwrite_phy(nid_tra,"ftsol2",itau_w,zx_tmp_fi2d) + zx_tmp_fi2d(:)=ftsol(:,3) + CALL histwrite_phy(nid_tra,"ftsol3",itau_w,zx_tmp_fi2d) + zx_tmp_fi2d(:)=ftsol(:,4) + CALL histwrite_phy(nid_tra,"ftsol4",itau_w,zx_tmp_fi2d) +! +! NATURE DU SOL + zx_tmp_fi2d(:)=pctsrf(:,1) + CALL histwrite_phy(nid_tra,"psrf1",itau_w,zx_tmp_fi2d) + zx_tmp_fi2d(:)=pctsrf(:,2) + CALL histwrite_phy(nid_tra,"psrf2",itau_w,zx_tmp_fi2d) + zx_tmp_fi2d(:)=pctsrf(:,3) + CALL histwrite_phy(nid_tra,"psrf3",itau_w,zx_tmp_fi2d) + zx_tmp_fi2d(:)=pctsrf(:,4) + CALL histwrite_phy(nid_tra,"psrf4",itau_w,zx_tmp_fi2d) + +! DIVERS + CALL histwrite_phy(nid_tra,"pplay",itau_w,pplay) + CALL histwrite_phy(nid_tra,"t",itau_w,t_seri) + CALL histwrite_phy(nid_tra,"mfu",itau_w,pmfu) + CALL histwrite_phy(nid_tra,"mfd",itau_w,pmfd) + CALL histwrite_phy(nid_tra,"en_u",itau_w,pen_u) + CALL histwrite_phy(nid_tra,"en_d",itau_w,pen_d) + CALL histwrite_phy(nid_tra,"de_d",itau_w,pde_d) + CALL histwrite_phy(nid_tra,"de_u",itau_w,pde_u) + CALL histwrite_phy(nid_tra,"coefh",itau_w,coefh) - IF (ecrit_tra>0. .AND. config_inca == 'none') THEN - ndex = 0 - ndex2d = 0 - ndex3d = 0 -c - itau_w = itau_phy + nstep + IF (ok_sync) THEN +!$OMP MASTER + CALL histsync(nid_tra) +!$OMP END MASTER + ENDIF - CALL histwrite_phy(nid_tra,"phis",itau_w,pphis) -C - CALL histwrite_phy(nid_tra,"aire",itau_w,airephy) + ENDIF !ecrit_tra>0. .AND. config_inca == 'none' - DO it=1,nbtr - iiq=niadv(it+2) - - CALL histwrite_phy(nid_tra,tname(iiq),itau_w,tr_seri(:,:,it)) - if (lessivage) THEN - CALL histwrite_phy(nid_tra,"fl"//tname(iiq),itau_w, - . flestottr(:,:,it)) - endif - -c----Olivia - CALL histwrite_phy(nid_tra,"d_tr_th_"//tname(iiq),itau_w, - . d_tr_th(:,:,it)) - - if(iflag_con.GE.2) then - CALL histwrite_phy(nid_tra,"d_tr_cv_"//tname(iiq),itau_w, - . d_tr_cv(:,:,it)) - endif !(iflag_con.GE.2) then - CALL histwrite_phy(nid_tra,"d_tr_cl_"//tname(iiq),itau_w, - . d_tr_cl(:,:,it)) -c---fin Olivia - - ENDDO - - -C abder - CALL histwrite_phy(nid_tra,"pyu1",itau_w,yu1) - - CALL histwrite_phy(nid_tra,"pyv1",itau_w,yv1) - - CALL histwrite_phy(nid_tra,"ftsol1",itau_w,pftsol1) - - CALL histwrite_phy(nid_tra,"ftsol2",itau_w,pftsol2) - - CALL histwrite_phy(nid_tra,"ftsol3",itau_w,pftsol3) - - CALL histwrite_phy(nid_tra,"ftsol4",itau_w,pftsol4) - - CALL histwrite_phy(nid_tra,"psrf1",itau_w,ppsrf1) - - CALL histwrite_phy(nid_tra,"psrf2",itau_w,ppsrf2) - - CALL histwrite_phy(nid_tra,"psrf3",itau_w,ppsrf3) - - CALL histwrite_phy(nid_tra,"psrf4",itau_w,ppsrf4) - CALL histwrite_phy(nid_tra,"pplay",itau_w,pplay) - - CALL histwrite_phy(nid_tra,"t",itau_w,t_seri) - CALL histwrite_phy(nid_tra,"mfu",itau_w,pmfu) - CALL histwrite_phy(nid_tra,"mfd",itau_w,pmfd) - CALL histwrite_phy(nid_tra,"en_u",itau_w,pen_u) - CALL histwrite_phy(nid_tra,"en_d",itau_w,pen_d) - CALL histwrite_phy(nid_tra,"de_d",itau_w,pde_d) - CALL histwrite_phy(nid_tra,"de_u",itau_w,pde_u) - CALL histwrite_phy(nid_tra,"coefh",itau_w,coefh) - - -c abder - - if (ok_sync) then -c$OMP MASTER - call histsync(nid_tra) -c$OMP END MASTER - endif - - END IF !ecrit_tra>0. .AND. config_inca == 'none' - - - +#endif