Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/infotrac.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/dyn3d/infotrac.F90 (revision 1178) +++ LMDZ4/branches/LMDZ4-dev/libf/dyn3d/infotrac.F90 (revision 1179) @@ -1,2 +1,4 @@ +! $Id$ +! MODULE infotrac @@ -299,4 +301,18 @@ END DO +! +! Test for advection schema. +! This version of LMDZ only garantees iadv=10 and iadv=14 (14 only for water vapour) . +! + DO iq=1,nqtot + IF (iadv(iq)/=10 .AND. iadv(iq)/=14) THEN + WRITE(lunout,*)'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ' + CALL abort_gcm('infotrac_init','In this version only iadv=10 and iadv=14 is tested!',1) + ELSE IF (iadv(iq)==14 .AND. iq/=1) THEN + WRITE(lunout,*)'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ' + CALL abort_gcm('infotrac_init','In this version iadv=14 is only permitted for water vapour!',1) + END IF + END DO + !----------------------------------------------------------------------- ! Finalize : Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/pres2lev.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/dyn3d/pres2lev.F (revision 1178) +++ LMDZ4/branches/LMDZ4-dev/libf/dyn3d/pres2lev.F (revision 1179) @@ -1,4 +1,3 @@ -! -! $Header$ +! $Id$ ! c****************************************************** @@ -21,19 +20,19 @@ c ARGUMENTS c """"""""" - LOGICAL ok_invertp - INTEGER lmo ! dimensions ancienne couches (input) - INTEGER lmn ! dimensions nouvelle couches (input) - INTEGER lmomx ! dimensions ancienne couches (input) - INTEGER lmnmx ! dimensions nouvelle couches (input) + LOGICAL, INTENT(IN) :: ok_invertp + INTEGER, INTENT(IN) :: lmo ! dimensions ancienne couches + INTEGER, INTENT(IN) :: lmn ! dimensions nouvelle couches + INTEGER lmomx ! dimensions ancienne couches + INTEGER lmnmx ! dimensions nouvelle couches parameter(lmomx=10000,lmnmx=10000) - real po(ni,nj,lmo)! niveau de pression ancienne grille - real pn(ni,nj,lmn) ! niveau de pression nouvelle grille + real, INTENT(IN) :: po(ni,nj,lmo) ! niveau de pression ancienne grille + real, INTENT(IN) :: pn(ni,nj,lmn) ! niveau de pression nouvelle grille - INTEGER i,j,Nhoriz,ni,nj ! nombre de point horizontale (input) + INTEGER, INTENT(IN) :: ni,nj ! nombre de point horizontale - REAL varo(ni,nj,lmo) ! var dans l'ancienne grille (input) - REAL varn(ni,nj,lmn) ! var dans la nouvelle grille (output) + REAL, INTENT(IN) :: varo(ni,nj,lmo) ! var dans l'ancienne grille + REAL, INTENT(OUT) :: varn(ni,nj,lmn) ! var dans la nouvelle grille real zvaro(lmomx),zpo(lmomx) @@ -41,5 +40,5 @@ c Autres variables c """""""""""""""" - INTEGER n, ln ,lo + INTEGER n, ln ,lo, i, j, Nhoriz REAL coef Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/infotrac.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/infotrac.F90 (revision 1178) +++ LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/infotrac.F90 (revision 1179) @@ -1,31 +1,27 @@ +! $Id$ +! MODULE infotrac ! nqtot : total number of tracers and higher order of moment, water vapor and liquid included INTEGER, SAVE :: nqtot -!!$OMP THREADPRIVATE(nqtot) ! nbtr : number of tracers not including higher order of moment or water vapor or liquid ! number of tracers used in the physics INTEGER, SAVE :: nbtr -!!$OMP THREADPRIVATE(nbtr) ! Name variables CHARACTER(len=20), ALLOCATABLE, DIMENSION(:), SAVE :: tname ! tracer short name for restart and diagnostics CHARACTER(len=23), ALLOCATABLE, DIMENSION(:), SAVE :: ttext ! tracer long name for diagnostics -!!$OMP THREADPRIVATE(tname,ttext) ! iadv : index of trasport schema for each tracer INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: iadv -!!$OMP THREADPRIVATE(iadv) ! niadv : vector keeping the coorspondance between all tracers(nqtot) treated in the ! dynamic part of the code and the tracers (nbtr+2) used in the physics part of the code. INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: niadv ! equivalent dyn / physique -!!$OMP THREADPRIVATE(niadv) ! Variables for INCA INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: conv_flg INTEGER, ALLOCATABLE, DIMENSION(:), SAVE :: pbl_flg -!!$OMP THREADPRIVATE(conv_flg, pbl_flg) CONTAINS @@ -305,4 +301,18 @@ END DO +! +! Test for advection schema. +! This version of LMDZ only garantees iadv=10 and iadv=14 (14 only for water vapour) . +! + DO iq=1,nqtot + IF (iadv(iq)/=10 .AND. iadv(iq)/=14) THEN + WRITE(lunout,*)'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ' + CALL abort_gcm('infotrac_init','In this version only iadv=10 and iadv=14 is tested!',1) + ELSE IF (iadv(iq)==14 .AND. iq/=1) THEN + WRITE(lunout,*)'STOP : The option iadv=',iadv(iq),' is not tested in this version of LMDZ' + CALL abort_gcm('infotrac_init','In this version iadv=14 is only permitted for water vapour!',1) + END IF + END DO + !----------------------------------------------------------------------- ! Finalize : Index: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/pres2lev.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/pres2lev.F (revision 1178) +++ LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/pres2lev.F (revision 1179) @@ -1,4 +1,3 @@ -! -! $Header$ +! $Id$ ! c****************************************************** @@ -21,19 +20,19 @@ c ARGUMENTS c """"""""" - LOGICAL ok_invertp - INTEGER lmo ! dimensions ancienne couches (input) - INTEGER lmn ! dimensions nouvelle couches (input) - INTEGER lmomx ! dimensions ancienne couches (input) - INTEGER lmnmx ! dimensions nouvelle couches (input) + LOGICAL, INTENT(IN) :: ok_invertp + INTEGER, INTENT(IN) :: lmo ! dimensions ancienne couches + INTEGER, INTENT(IN) :: lmn ! dimensions nouvelle couches + INTEGER lmomx ! dimensions ancienne couches + INTEGER lmnmx ! dimensions nouvelle couches parameter(lmomx=10000,lmnmx=10000) - real po(ni,nj,lmo)! niveau de pression ancienne grille - real pn(ni,nj,lmn) ! niveau de pression nouvelle grille + real, INTENT(IN) :: po(ni,nj,lmo) ! niveau de pression ancienne grille + real, INTENT(IN) :: pn(ni,nj,lmn) ! niveau de pression nouvelle grille - INTEGER i,j,Nhoriz,ni,nj ! nombre de point horizontale (input) + INTEGER, INTENT(IN) :: ni,nj ! nombre de point horizontale - REAL varo(ni,nj,lmo) ! var dans l'ancienne grille (input) - REAL varn(ni,nj,lmn) ! var dans la nouvelle grille (output) + REAL, INTENT(IN) :: varo(ni,nj,lmo) ! var dans l'ancienne grille + REAL, INTENT(OUT) :: varn(ni,nj,lmn) ! var dans la nouvelle grille real zvaro(lmomx),zpo(lmomx) @@ -41,5 +40,5 @@ c Autres variables c """""""""""""""" - INTEGER n, ln ,lo + INTEGER n, ln ,lo, i, j, Nhoriz REAL coef Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/aerosol_optic.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/aerosol_optic.F90 (revision 1178) +++ (revision ) @@ -1,145 +1,0 @@ -SUBROUTINE aerosol_optic(debut, new_aod, flag_aerosol, rjourvrai, pdtphys, & - pplay, paprs, t_seri, rhcl, & - maerosol, maerosol_pi, & - tau_aero, piz_aero, cg_aero ) - - USE dimphy - IMPLICIT NONE - -! Input arguments - LOGICAL, INTENT(IN) :: debut - LOGICAL, INTENT(IN) :: new_aod - INTEGER, INTENT(IN) :: flag_aerosol - REAL, INTENT(IN) :: rjourvrai - REAL, INTENT(IN) :: pdtphys - REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay - REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs - REAL, DIMENSION(klon,klev), INTENT(IN) :: t_seri - REAL, DIMENSION(klon,klev), INTENT(IN) :: rhcl ! humidite relative ciel clair - -! Output arguments - REAL, DIMENSION(klon,klev), INTENT(OUT) :: maerosol - REAL, DIMENSION(klon,klev), INTENT(OUT) :: maerosol_pi - REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: tau_aero - REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: piz_aero - REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: cg_aero - -! Local variables - REAL, DIMENSION(klon) :: aerindex ! POLDER aerosol index - REAL, DIMENSION(klon,10) :: fractnat_allaer - REAL, DIMENSION(klon,klev) :: sulfate ! SO4 aerosol concentration [ug/m3] - REAL, DIMENSION(klon,klev) :: bcsol ! BC soluble concentration [ug/m3] - REAL, DIMENSION(klon,klev) :: bcins ! BC insoluble concentration [ug/m3] - REAL, DIMENSION(klon,klev) :: pomsol ! POM soluble concentration [ug/m3] - REAL, DIMENSION(klon,klev) :: pomins ! POM insoluble concentration [ug/m3] - REAL, DIMENSION(klon,klev) :: sulfate_pi - REAL, DIMENSION(klon,klev) :: bcsol_pi - REAL, DIMENSION(klon,klev) :: bcins_pi - REAL, DIMENSION(klon,klev) :: pomsol_pi - REAL, DIMENSION(klon,klev) :: pomins_pi - REAL, DIMENSION(klon,klev) :: pdel - REAL, DIMENSION(klon,klev,8) :: m_allaer - - INTEGER :: k - -! -! Read aersosols from file -! - maerosol(:,:) = 0. - maerosol_pi(:,:) = 0. - - bcsol(:,:) = 0. - bcins(:,:) = 0. - pomsol(:,:) = 0. - pomins(:,:) = 0. - sulfate(:,:) = 0. - -! Read sulfate - IF ( flag_aerosol .EQ. 1 .OR. & - flag_aerosol .EQ. 4 .OR. & - flag_aerosol .EQ. 6 ) THEN - - CALL readaerosol(5,rjourvrai, debut, sulfate) - CALL readaerosol_preind(5,rjourvrai, & - debut, sulfate_pi) - - maerosol(:,:) = maerosol(:,:) + sulfate(:,:) - maerosol_pi(:,:) = maerosol_pi(:,:) + sulfate_pi(:,:) - ENDIF - - -! Read bcscol and bcins - IF ( flag_aerosol .EQ. 2 .OR. & - flag_aerosol .EQ. 4 .OR. & - flag_aerosol .EQ. 5 ) THEN - - ! Get bc aerosol distribution - CALL readaerosol(3,rjourvrai, debut,bcsol ) - CALL readaerosol_preind(3,rjourvrai, debut, bcsol_pi) - - CALL readaerosol(7,rjourvrai, debut,bcins ) - CALL readaerosol_preind(7,rjourvrai, debut, bcins_pi) - - maerosol(:,:) = maerosol(:,:) + bcsol(:,:) - maerosol_pi(:,:) = maerosol_pi(:,:) + bcsol_pi(:,:) - ENDIF - - -! Read pomsol and pomins - IF ( flag_aerosol .EQ. 3 .OR. & - flag_aerosol .EQ. 4 .OR. & - flag_aerosol .EQ. 5 .OR. & - flag_aerosol .EQ. 6 ) THEN - - CALL readaerosol(4,rjourvrai, debut,pomsol) - CALL readaerosol_preind(4,rjourvrai,debut,pomsol_pi ) - - CALL readaerosol(8,rjourvrai, debut,pomins) - CALL readaerosol_preind(8,rjourvrai,debut,pomins_pi ) - - maerosol(:,:) = maerosol(:,:) + pomsol - maerosol_pi(:,:) = maerosol_pi(:,:) + pomsol_pi - ENDIF - - -! Save all aerosols in one variable -! -! ACo pour couplage aerosol offline 07/04/2009 -! Tableau contenant les masses pour tous les aerosols -! les valeurs a zero seront a remplacer par les bons -! tableaux lorsque les routines de lectures seront -! ajoutees. - m_allaer(:,:,1) = 0. ! SSSSM || CSSSM - m_allaer(:,:,2) = 0. ! ASSSM - m_allaer(:,:,3) = bcsol(:,:) ! ASBCM - m_allaer(:,:,4) = pomsol(:,:) ! ASPOMM - m_allaer(:,:,5) = sulfate(:,:) ! ASSO4M || CSSO4M - m_allaer(:,:,6) = 0. ! CIDUSTM - m_allaer(:,:,7) = bcins(:,:) ! AIBCM - m_allaer(:,:,8) = pomins(:,:) ! AIPOMM - -! -! Calculate the optical properties for the aerosols -! - DO k = 1, klev - pdel(:,k) = paprs(:,k) - paprs (:,k+1) - END DO - - IF (.NOT. new_aod) THEN - CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl, & - tau_aero(:,:,5,:), piz_aero(:,:,5,:), cg_aero(:,:,5,:), aerindex) - ELSE - fractnat_allaer(:,:) = 0. - CALL aeropt_2bands( & - pdel, m_allaer, pdtphys, rhcl, & - tau_aero, piz_aero, cg_aero, & - fractnat_allaer, flag_aerosol, & - pplay, t_seri) - - CALL aeropt_5wv( & - pdel, m_allaer, & - pdtphys, rhcl, aerindex, & - flag_aerosol, pplay, t_seri) - ENDIF - -END SUBROUTINE aerosol_optic Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/newmicro.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/newmicro.F (revision 1178) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/newmicro.F (revision 1179) @@ -1,5 +1,4 @@ -! -! $Header$ -! +! $Id$ +! SUBROUTINE newmicro (paprs, pplay,ok_newmicro, . t, pqlwp, pclc, pcltau, pclemi, @@ -7,5 +6,5 @@ s xflwp, xfiwp, xflwc, xfiwc, e ok_aie, - e sulfate, sulfate_pi, + e mass_ins_aero, mass_ins_aero_pi, e bl95_b0, bl95_b1, s cldtaupi, re, fl) @@ -22,6 +21,6 @@ c c ok_aie--input-L-apply aerosol indirect effect or not -c sulfate-input-R-sulfate aerosol mass concentration [um/m^3] -c sulfate_pi-input-R-dito, pre-industrial value +c mass_ins_aero-----input-R-total mass concentration for all indissoluble aerosols[ug/m^3] +c mass_ins_aero_pi--input-R-dito, pre-industrial value c bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land) c bl95_b1-input-R-a parameter, may be varied for tests ( -"- ) @@ -94,8 +93,8 @@ LOGICAL ok_a1lwpdep ! a1 LWP dependent? - REAL sulfate(klon, klev) ! sulfate aerosol mass concentration [ug m-3] + REAL mass_ins_aero(klon, klev) ! total mass concentration for all indissoluble aerosols [ug m-3] + REAL mass_ins_aero_pi(klon, klev) ! - " - (pre-industrial value) REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3] REAL re(klon, klev) ! cloud droplet effective radius [um] - REAL sulfate_pi(klon, klev) ! sulfate aerosol mass concentration [ug m-3] (pre-industrial value) REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value) REAL re_pi(klon, klev) ! cloud droplet effective radius [um] (pi value) @@ -157,5 +156,5 @@ ! cdnc(i,k) = 10.**(bl95_b0+bl95_b1* - & log(MAX(sulfate(i,k),1.e-4))/log(10.))*1.e6 !-m-3 + & log(MAX(mass_ins_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3 ! Cloud droplet number concentration (CDNC) is restricted ! to be within [20, 1000 cm^3] @@ -165,5 +164,5 @@ ! cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1* - & log(MAX(sulfate_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3 + & log(MAX(mass_ins_aero_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3 cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k))) ENDDO Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/nuage.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/nuage.F (revision 1178) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/nuage.F (revision 1179) @@ -1,4 +1,3 @@ -! -! $Header$ +! $Id$ ! SUBROUTINE nuage (paprs, pplay, @@ -6,5 +5,5 @@ . pch, pcl, pcm, pct, pctlwp, e ok_aie, - e sulfate, sulfate_pi, + e mass_ins_aero, mass_ins_aero_pi, e bl95_b0, bl95_b1, s cldtaupi, re, fl) @@ -20,6 +19,6 @@ c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1) c ok_aie--input-L-apply aerosol indirect effect or not -c sulfate-input-R-sulfate aerosol mass concentration [um/m^3] -c sulfate_pi-input-R-dito, pre-industrial value +c mass_ins_aero-----input-R-total mass concentration for all indissoluble aerosols[ug/m^3] +c mass_ins_aero_pi--input-R-dito, pre-industrial value c bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land) c bl95_b1-input-R-a parameter, may be varied for tests ( -"- ) @@ -74,8 +73,8 @@ LOGICAL ok_aie ! Apply AIE or not? - REAL sulfate(klon, klev) ! sulfate aerosol mass concentration [ug m-3] + REAL mass_ins_aero(klon, klev) ! total mass concentration for all indissoluble aerosols[ug m-3] + REAL mass_ins_aero_pi(klon, klev) ! - " - pre-industrial value REAL cdnc(klon, klev) ! cloud droplet number concentration [m-3] REAL re(klon, klev) ! cloud droplet effective radius [um] - REAL sulfate_pi(klon, klev) ! sulfate aerosol mass concentration [ug m-3] (pre-industrial value) REAL cdnc_pi(klon, klev) ! cloud droplet number concentration [m-3] (pi value) REAL re_pi(klon, klev) ! cloud droplet effective radius [um] (pi value) @@ -108,5 +107,5 @@ ! cdnc(i,k) = 10.**(bl95_b0+bl95_b1* - . log(MAX(sulfate(i,k),1.e-4))/log(10.))*1.e6 !-m-3 + . log(MAX(mass_ins_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3 ! Cloud droplet number concentration (CDNC) is restricted ! to be within [20, 1000 cm^3] @@ -114,5 +113,5 @@ cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k))) cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1* - . log(MAX(sulfate_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3 + . log(MAX(mass_ins_aero_pi(i,k),1.e-4))/log(10.))*1.e6 !-m-3 cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k))) ! Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_local_var_mod.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_local_var_mod.F90 (revision 1178) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_local_var_mod.F90 (revision 1179) @@ -69,4 +69,16 @@ REAL, SAVE, ALLOCATABLE :: d_ts(:,:), d_tr(:,:,:) !$OMP THREADPRIVATE(d_ts, d_tr) + +! diagnostique pour le rayonnement + REAL, SAVE, ALLOCATABLE :: topswad_aero(:), solswad_aero(:) ! diag + !$OMP THREADPRIVATE(topswad_aero,solswad_aero) + REAL, SAVE, ALLOCATABLE :: topswai_aero(:), solswai_aero(:) ! diag + !$OMP THREADPRIVATE(topswai_aero,solswai_aero) + REAL, SAVE, ALLOCATABLE :: topswad0_aero(:), solswad0_aero(:) ! pas utilise, eventuellment pour diag + !$OMP THREADPRIVATE(topswad0_aero,solswad0_aero) + REAL, SAVE, ALLOCATABLE :: topsw_aero(:,:), solsw_aero(:,:) ! pas utilise, eventuellment pour diag + !$OMP THREADPRIVATE(topsw_aero,solsw_aero) + REAL, SAVE, ALLOCATABLE :: topsw0_aero(:,:), solsw0_aero(:,:) ! pas utilise, eventuellment pour diag + !$OMP THREADPRIVATE(topsw0_aero,solsw0_aero) CONTAINS @@ -100,4 +112,9 @@ allocate(d_u_lif(klon,klev),d_v_lif(klon,klev)) allocate(d_ts(klon,klev), d_tr(klon,klev,nbtr)) + allocate(topswad_aero(klon), solswad_aero(klon)) + allocate(topswai_aero(klon), solswai_aero(klon)) + allocate(topswad0_aero(klon), solswad0_aero(klon)) + allocate(topsw_aero(klon,9), solsw_aero(klon,9)) + allocate(topsw0_aero(klon,9), solsw0_aero(klon,9)) allocate(d_u_hin(klon,klev),d_v_hin(klon,klev),d_t_hin(klon,klev)) @@ -132,4 +149,9 @@ deallocate(d_u_lif,d_v_lif) deallocate(d_ts, d_tr) + deallocate(topswad_aero,solswad_aero) + deallocate(topswai_aero,solswai_aero) + deallocate(topswad0_aero,solswad0_aero) + deallocate(topsw_aero,solsw_aero) + deallocate(topsw0_aero,solsw0_aero) deallocate(d_u_hin,d_v_hin,d_t_hin) Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_state_var_mod.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_state_var_mod.F90 (revision 1178) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_state_var_mod.F90 (revision 1179) @@ -274,4 +274,6 @@ !$OMP THREADPRIVATE(topswai,solswai) + REAL,SAVE,ALLOCATABLE :: tau_aero(:,:,:,:), piz_aero(:,:,:,:), cg_aero(:,:,:,:) +!$OMP THREADPRIVATE(tau_aero, piz_aero, cg_aero) REAL,SAVE,ALLOCATABLE :: ccm(:,:,:) !$OMP THREADPRIVATE(ccm) @@ -379,5 +381,5 @@ ALLOCATE(topswad(klon), solswad(klon)) ALLOCATE(topswai(klon), solswai(klon)) - + ALLOCATE(tau_aero(klon,klev,9,2),piz_aero(klon,klev,9,2),cg_aero(klon,klev,9,2)) ALLOCATE(ccm(klon,klev,2)) @@ -466,5 +468,5 @@ deallocate(topswad, solswad) deallocate(topswai, solswai) - + deallocate(tau_aero,piz_aero,cg_aero) deallocate(ccm) Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F (revision 1178) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F (revision 1179) @@ -1,3 +1,2 @@ -! ! $Id$ ! @@ -537,11 +536,4 @@ c c Variables propres a la physique -c -c INTEGER radpas -c SAVE radpas ! frequence d'appel rayonnement -ccccccccc$OMP THREADPRIVATE(radpas) -c -cc INTEGER iflag_con -c INTEGER itap SAVE itap ! compteur pour la physique @@ -1074,13 +1066,6 @@ CHARACTER*4, DIMENSION(9) :: rfname REAL, DIMENSION(klon) :: aerindex ! POLDER aerosol index - REAL, DIMENSION(klon,klev) :: maerosol ! aerosol concentration [ug/m3] - REAL, DIMENSION(klon,klev) :: maerosol_pi ! aerosol concentration [ug/m3] (pre-industrial value) - REAL, DIMENSION(klon,klev,9,2) :: tau_aero, piz_aero, cg_aero - REAL, DIMENSION(klon) :: topswad_aero, solswad_aero ! diag - REAL, DIMENSION(klon) :: topswai_aero, solswai_aero ! diag - REAL, DIMENSION(klon) :: topswad0_aero, solswad0_aero ! pas utilise, eventuellment pour diag - REAL, DIMENSION(klon,9) :: topsw_aero, solsw_aero ! pas utilise - REAL, DIMENSION(klon,9) :: topsw0_aero, solsw0_aero ! pas utilise - + REAL, DIMENSION(klon,klev) :: mass_ins_aero! total mass concentration for all indissoluble aerosols[ug/m3] + REAL, DIMENSION(klon,klev) :: mass_ins_aero_pi ! - " - (pre-industrial value) ! Parameters @@ -1228,5 +1213,11 @@ tau_overturning_th(:)=0. - IF (config_inca /= 'none') ccm(:,:,:) = 0. + IF (config_inca /= 'none') THEN + ! jg : initialisation jusqu'au ces variables sont dans restart + ccm(:,:,:) = 0. + tau_aero(:,:,:,:) = 0. + piz_aero(:,:,:,:) = 0. + cg_aero(:,:,:,:) = 0. + END IF rnebcon0(:,:) = 0.0 @@ -2644,8 +2635,8 @@ IF (ok_ade.OR.ok_aie) THEN IF (.NOT. aerosol_couple) - & CALL aerosol_optic( + & CALL readaerosol_optic( & debut, new_aod, flag_aerosol, rjourvrai, pdtphys, & pplay, paprs, t_seri, rhcl, - & maerosol, maerosol_pi, + & mass_ins_aero, mass_ins_aero_pi, & tau_aero, piz_aero, cg_aero ) ELSE @@ -2829,6 +2820,6 @@ IF (aerosol_couple) THEN - maerosol(:,:) = ccm(:,:,1) - maerosol_pi(:,:) = ccm(:,:,2) + mass_ins_aero(:,:) = ccm(:,:,1) + mass_ins_aero_pi(:,:) = ccm(:,:,2) END IF @@ -2839,5 +2830,5 @@ . flwp, fiwp, flwc, fiwc, e ok_aie, - e maerosol, maerosol_pi, + e mass_ins_aero, mass_ins_aero_pi, e bl95_b0, bl95_b1, s cldtaupi, re, fl) @@ -2847,5 +2838,5 @@ . cldh, cldl, cldm, cldt, cldq, e ok_aie, - e maerosol, maerosol_pi, + e mass_ins_aero, mass_ins_aero_pi, e bl95_b0, bl95_b1, s cldtaupi, re, fl) @@ -2922,7 +2913,7 @@ - ENDIF + ENDIF ! aerosol_couple itaprad = 0 - ENDIF + ENDIF ! MOD(itaprad,radpas) itaprad = itaprad + 1 Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90 (revision 1178) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90 (revision 1179) @@ -1,700 +1,519 @@ +! $Id$ ! -! $Header$ +MODULE readaerosol_mod + + REAL, SAVE :: not_valid=-333. + +CONTAINS + +SUBROUTINE readaerosol(name_aero, type, iyr_in, klev_src, pt_ap, pt_b, pt_out, psurf, load) + +!**************************************************************************************** +! This routine will read the aersosol from file. ! -SUBROUTINE readaerosol (id_aero, r_day, first, massvar_p) +! Read a year data with get_aero_fromfile depending on aer_type : +! - actuel : read year 1980 +! - preind : read natural data +! - scenario : read one or two years and do eventually linare time interpolation +! +! Return pointer, pt_out, to the year read or result from interpolation +!**************************************************************************************** + USE dimphy + + IMPLICIT NONE + + INCLUDE "iniprint.h" + + ! Input arguments + CHARACTER(len=7), INTENT(IN) :: name_aero + CHARACTER(len=*), INTENT(IN) :: type ! correspond to aer_type in clesphys.h + INTEGER, INTENT(IN) :: iyr_in + + ! Output + INTEGER, INTENT(OUT) :: klev_src + REAL, POINTER, DIMENSION(:) :: pt_ap ! Pointer for describing the vertical levels + REAL, POINTER, DIMENSION(:) :: pt_b ! Pointer for describing the vertical levels + REAL, POINTER, DIMENSION(:,:,:) :: pt_out ! The massvar distributions, DIMENSION(klon, klev_src, 12) + REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf ! Surface pression for 12 months + REAL, DIMENSION(klon,12), INTENT(OUT) :: load ! Aerosol mass load in each column for 12 months + + ! Local variables + CHARACTER(len=4) :: cyear + REAL, POINTER, DIMENSION(:,:,:) :: pt_2 + REAL, DIMENSION(klon,12) :: psurf2, load2 + REAL :: p0 ! Reference pressure + INTEGER :: iyr1, iyr2, klev_src2 + INTEGER :: it, k, i + LOGICAL, PARAMETER :: lonlyone=.FALSE. + +!**************************************************************************************** +! Read data depending on aer_type +! +!**************************************************************************************** + + IF (type == 'actuel') THEN +! Read and return data for year 1980 +!**************************************************************************************** + cyear='1980' + ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month + ! pt_out has dimensions (klon, klev_src, 12) + CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load) + + + ELSE IF (type == 'preind') THEN +! Read and return data from file with suffix .nat +!**************************************************************************************** + cyear='.nat' + ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month + ! pt_out has dimensions (klon, klev_src, 12) + CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load) + + ELSE IF (type == 'scenario') THEN +! Read data depending on actual year and interpolate if necessary +!**************************************************************************************** + IF (iyr_in .LT. 1850) THEN + cyear='.nat' + WRITE(lunout,*) 'get_aero 1 iyr_in=', iyr_in,' ',cyear + ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month + ! pt_out has dimensions (klon, klev_src, 12) + CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load) + + ELSE IF (iyr_in .GE. 2100) THEN + cyear='2100' + WRITE(lunout,*) 'get_aero 2 iyr_in=', iyr_in,' ',cyear + ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month + ! pt_out has dimensions (klon, klev_src, 12) + CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load) + + ELSE + ! Read data from 2 decades and interpolate to actual year + ! a) from actual 10-yr-period + IF (iyr_in.LT.1900) THEN + iyr1 = 1850 + iyr2 = 1900 + ELSE IF (iyr_in.GE.1900.AND.iyr_in.LT.1920) THEN + iyr1 = 1900 + iyr2 = 1920 + ELSE + iyr1 = INT(iyr_in/10)*10 + iyr2 = INT(1+iyr_in/10)*10 + ENDIF + + WRITE(cyear,'(I4)') iyr1 + WRITE(lunout,*) 'get_aero 3 iyr_in=', iyr_in,' ',cyear + ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month + ! pt_out has dimensions (klon, klev_src, 12) + CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load) + + ! If to read two decades: + IF (.NOT.lonlyone) THEN + + ! b) from the next following one + WRITE(cyear,'(I4)') iyr2 + WRITE(lunout,*) 'get_aero 4 iyr_in=', iyr_in,' ',cyear + + NULLIFY(pt_2) + ! get_aero_fromfile returns pt_2 allocated and initialized with data for 12 month + ! pt_2 has dimensions (klon, klev_src, 12) + CALL get_aero_fromfile(name_aero, cyear, klev_src2, pt_ap, pt_b, p0, pt_2, psurf2, load2) + ! Test for same number of vertical levels + IF (klev_src /= klev_src2) THEN + WRITE(lunout,*) 'Two aerosols files with different number of vertical levels is not allowded' + CALL abort_gcm('readaersosol','Error in number of vertical levels',1) + END IF + + ! Linare interpolate to the actual year: + DO it=1,12 + DO k=1,klev_src + DO i = 1, klon + pt_out(i,k,it) = & + pt_out(i,k,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * & + (pt_out(i,k,it) - pt_2(i,k,it)) + END DO + END DO + + DO i = 1, klon + psurf(i,it) = & + psurf(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * & + (psurf(i,it) - psurf2(i,it)) + + load(i,it) = & + load(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * & + (load(i,it) - load2(i,it)) + END DO + END DO + + ! Deallocate pt_2 no more needed + DEALLOCATE(pt_2) + + END IF ! lonlyone + END IF ! iyr_in .LT. 1850 + + ELSE + WRITE(lunout,*)'This option is not implemented : aer_type = ', type + CALL abort_gcm('readaerosol','Error : aer_type parameter not accepted',1) + END IF ! type + + +END SUBROUTINE readaerosol + + + SUBROUTINE get_aero_fromfile(varname, cyr, klev_src, pt_ap, pt_b, p0, pt_year, psurf_out, load_out) +!**************************************************************************************** +! Read 12 month aerosol from file and distribute to local process on physical grid. +! Vertical levels, klev_src, may differ from model levels if new file format. +! +! For mpi_root and master thread : +! 1) Open file +! 2) Find vertical dimension klev_src +! 3) Read field month by month +! 4) Close file +! 5) Transform the global field from 2D(iim, jjp+1) to 1D(klon_glo) +! - Also the levels and the latitudes have to be inversed +! +! For all processes and threads : +! 6) Scatter global field(klon_glo) to local process domain(klon) +! 7) Test for negative values +!**************************************************************************************** + + USE netcdf + USE dimphy + USE mod_grid_phy_lmdz + USE mod_phys_lmdz_para + + IMPLICIT NONE + + INCLUDE "dimensions.h" + INCLUDE "iniprint.h" + +! Input argumets + CHARACTER(len=7), INTENT(IN) :: varname + CHARACTER(len=4), INTENT(IN) :: cyr + +! Output arguments + INTEGER, INTENT(OUT) :: klev_src ! Number of vertical levels in file + REAL, POINTER, DIMENSION(:) :: pt_ap ! Pointer for describing the vertical levels + REAL, POINTER, DIMENSION(:) :: pt_b ! Pointer for describing the vertical levels + REAL :: p0 ! Reference pressure value + REAL, POINTER, DIMENSION(:,:,:) :: pt_year ! Pointer-variabale from file, 12 month, grid : klon,klev_src + REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf_out ! Surface pression for 12 months + REAL, DIMENSION(klon,12), INTENT(OUT) :: load_out ! Aerosol mass load in each column + +! Local variables + CHARACTER(len=30) :: fname + CHARACTER(len=30) :: cvar + INTEGER :: ncid, dimid, varid + INTEGER :: imth, i, j, k, ierr + REAL :: npole, spole + REAL, ALLOCATABLE, DIMENSION(:,:,:) :: varmth + REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear ! Global variable read from file, 12 month + REAL, ALLOCATABLE, DIMENSION(:,:,:) :: varyear_glo1D !(klon_glo, klev_src, 12) + REAL, ALLOCATABLE, DIMENSION(:) :: varktmp + + REAL, DIMENSION(iim,jjm+1,12) :: psurf_glo2D ! Surface pression for 12 months on dynamics global grid + REAL, DIMENSION(klon_glo,12) :: psurf_glo1D ! -"- on physical global grid + REAL, DIMENSION(iim,jjm+1,12) :: load_glo2D ! Load for 12 months on dynamics global grid + REAL, DIMENSION(klon_glo,12) :: load_glo1D ! -"- on physical global grid + REAL, DIMENSION(iim,jjm+1) :: vartmp + LOGICAL :: new_file + + + ! Deallocate pointers + IF (ASSOCIATED(pt_ap)) DEALLOCATE(pt_ap) + IF (ASSOCIATED(pt_b)) DEALLOCATE(pt_b) + +!$OMP MASTER + IF (is_mpi_root) THEN + +! 1) Open file +!**************************************************************************************** + fname = TRIM(varname)//'.run'//cyr//'.nc' - USE dimphy, ONLY : klev - USE mod_grid_phy_lmdz, klon=>klon_glo - USE mod_phys_lmdz_para + WRITE(lunout,*) 'reading ', TRIM(fname) + CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid) ) + +! 2) Check if old or new file is avalabale. +! New type of file should contain the dimension 'lev' +! Old type of file should contain the dimension 'PRESNIVS' +!**************************************************************************************** + ierr = nf90_inq_dimid(ncid, 'lev', dimid) + IF (ierr /= NF90_NOERR) THEN + ! Coordinate axe lev not found. Check for presnivs. + ierr = nf90_inq_dimid(ncid, 'PRESNIVS', dimid) + IF (ierr /= NF90_NOERR) THEN + ! Dimension PRESNIVS not found either + CALL abort_gcm('get_aero_fromfile', 'dimension lev or presnivs not in file',1) + ELSE + ! Old file found + new_file=.FALSE. + WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will not be done' + END IF + ELSE + ! New file found + new_file=.TRUE. + WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done' + END IF + +! 2) Find vertical dimension klev_src +!**************************************************************************************** + CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src) ) + + ! Allocate variables depending on the number of vertical levels + ALLOCATE(varmth(iim, jjm+1, klev_src), varyear(iim, jjm+1, klev_src, 12), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 1',1) + + ALLOCATE(pt_ap(klev_src), pt_b(klev_src), varktmp(klev_src), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 2',1) + +! 3) Read all variables from file +! There is 2 options for the file structure : +! new_file=TRUE : read varyear, ps, pt_ap and pt_b +! new_file=FALSE : read varyear month by month +!**************************************************************************************** + + IF (new_file) THEN + +! ++) Read the aerosol concentration month by month and concatenate to total variable varyear +!**************************************************************************************** + ! Get variable id + CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid) ) + + ! Get the variable + CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)) ) + +! ++) Read surface pression, 12 month in one variable +!**************************************************************************************** + ! Get variable id + CALL check_err( nf90_inq_varid(ncid, "ps", varid) ) + ! Get the variable + CALL check_err( nf90_get_var(ncid, varid, psurf_glo2D) ) + +! ++) Read mass load, 12 month in one variable +!**************************************************************************************** + ! Get variable id + CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) ) + ! Get the variable + CALL check_err( nf90_get_var(ncid, varid, load_glo2D) ) + +! ++) Read ap +!**************************************************************************************** + ! Get variable id + CALL check_err( nf90_inq_varid(ncid, "ap", varid) ) + ! Get the variable + CALL check_err( nf90_get_var(ncid, varid, pt_ap) ) + +! ++) Read b +!**************************************************************************************** + ! Get variable id + CALL check_err( nf90_inq_varid(ncid, "b", varid) ) + ! Get the variable + CALL check_err( nf90_get_var(ncid, varid, pt_b) ) + +! ++) Read p0 : reference pressure +!**************************************************************************************** + ! Get variable id + CALL check_err( nf90_inq_varid(ncid, "p0", varid) ) + ! Get the variable + CALL check_err( nf90_get_var(ncid, varid, p0) ) + + + ELSE ! old file + +! ++) Read the aerosol concentration month by month and concatenate to total variable varyear +!**************************************************************************************** + DO imth=1, 12 + IF (imth.EQ.1) THEN + cvar=TRIM(varname)//'JAN' + ELSE IF (imth.EQ.2) THEN + cvar=TRIM(varname)//'FEB' + ELSE IF (imth.EQ.3) THEN + cvar=TRIM(varname)//'MAR' + ELSE IF (imth.EQ.4) THEN + cvar=TRIM(varname)//'APR' + ELSE IF (imth.EQ.5) THEN + cvar=TRIM(varname)//'MAY' + ELSE IF (imth.EQ.6) THEN + cvar=TRIM(varname)//'JUN' + ELSE IF (imth.EQ.7) THEN + cvar=TRIM(varname)//'JUL' + ELSE IF (imth.EQ.8) THEN + cvar=TRIM(varname)//'AUG' + ELSE IF (imth.EQ.9) THEN + cvar=TRIM(varname)//'SEP' + ELSE IF (imth.EQ.10) THEN + cvar=TRIM(varname)//'OCT' + ELSE IF (imth.EQ.11) THEN + cvar=TRIM(varname)//'NOV' + ELSE IF (imth.EQ.12) THEN + cvar=TRIM(varname)//'DEC' + END IF + + ! Get variable id + CALL check_err( nf90_inq_varid(ncid, TRIM(cvar), varid) ) + + ! Get the variable + CALL check_err( nf90_get_var(ncid, varid, varmth) ) + + ! Store in variable for the whole year + varyear(:,:,:,imth)=varmth(:,:,:) + + END DO + + ! Putting dummy + psurf_glo2D(:,:,:) = not_valid + load_glo2D(:,:,:) = not_valid + pt_ap(:) = not_valid + pt_b(:) = not_valid + + END IF + +! 4) Close file +!**************************************************************************************** + CALL check_err( nf90_close(ncid) ) + + +! 5) Transform the global field from 2D(iim, jjp+1) to 1D(klon_glo) +!**************************************************************************************** +! Test if vertical levels have to be inversed + + IF ((pt_b(1) < pt_b(klev_src)) .OR. .NOT. new_file) THEN + WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' needs to be inversed' + WRITE(lunout,*) 'before pt_ap = ', pt_ap + WRITE(lunout,*) 'before pt_b = ', pt_b + + ! Inverse vertical levels for varyear + DO imth=1, 12 + varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly + DO k=1, klev_src + DO j=1, jjm+1 + DO i=1,iim + varyear(i,j,k,imth) = varmth(i,j,klev_src+1-k) + END DO + END DO + END DO + END DO + + ! Inverte vertical axes for pt_ap and pt_b + varktmp(:) = pt_ap(:) + DO k=1, klev_src + pt_ap(k) = varktmp(klev_src+1-k) + END DO + + varktmp(:) = pt_b(:) + DO k=1, klev_src + pt_b(k) = varktmp(klev_src+1-k) + END DO + WRITE(lunout,*) 'after pt_ap = ', pt_ap + WRITE(lunout,*) 'after pt_b = ', pt_b + + ELSE + WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done' + WRITE(lunout,*) 'pt_ap = ', pt_ap + WRITE(lunout,*) 'pt_b = ', pt_b + END IF + +! - Also the latitudes have to be inversed + DO imth=1, 12 + ! Inverse latitudes + varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly + DO k=1,klev_src + DO j=1,jjm+1 + DO i=1,iim + varyear(i,j,k,imth) = varmth(i,jjm+1+1-j,k) + END DO + END DO + END DO + + ! Inverse latitudes for surface pressure + vartmp(:,:) = psurf_glo2D(:,:,imth) + DO j=1, jjm+1 + DO i=1,iim + psurf_glo2D(i,j,imth)= vartmp(i,jjm+1+1-j) + END DO + END DO + + ! Inverse latitudes for the load + vartmp(:,:) = load_glo2D(:,:,imth) + DO j=1, jjm+1 + DO i=1,iim + load_glo2D(i,j,imth)= vartmp(i,jjm+1+1-j) + END DO + END DO + + ! Do zonal mead at poles and distribut at whole first and last latitude + DO k=1, klev_src + npole=0. ! North pole, j=1 + spole=0. ! South pole, j=jjm+1 + DO i=1,iim + npole = npole + varyear(i,1,k,imth) + spole = spole + varyear(i,jjm+1,k,imth) + END DO + npole = npole/FLOAT(iim) + spole = spole/FLOAT(iim) + varyear(:,1, k,imth) = npole + varyear(:,jjm+1,k,imth) = spole + END DO + END DO ! imth + + ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 3',1) + + ! Transform from 2D to 1D field + CALL grid2Dto1D_glo(varyear,varyear_glo1D) + CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D) + CALL grid2Dto1D_glo(load_glo2D,load_glo1D) + + END IF ! is_mpi_root +!$OMP END MASTER BARRIER - IMPLICIT NONE - -! Content: -! -------- -! This routine reads in monthly mean values of massvar aerosols and -! returns a linearly interpolated dayly-mean field. -! -! -! Author: -! ------- -! Johannes Quaas (quaas@lmd.jussieu.fr) -! Celine Deandreis & Anne Cozic -! 19/02/09 -! - -! -! Include-files: -! -------------- - INCLUDE "YOMCST.h" - INCLUDE "chem.h" - INCLUDE "dimensions.h" - INCLUDE "temps.h" - INCLUDE "clesphys.h" - INCLUDE "iniprint.h" -! -! Input: -! ------ - INTEGER, INTENT(in) :: id_aero - REAL, INTENT(in) :: r_day ! Day of integration - LOGICAL, INTENT(in) :: first ! First timestep - ! (and therefore initialization necessary) -! -! Output: -! ------- - REAL, INTENT(out) :: massvar_p(klon_omp,klev) ! Mass of massvar (monthly mean data,from file) [ug AIBCM/m3] - -! -! Local Variables: -! ---------------- - INTEGER :: i, ig, k, it - INTEGER :: j, iday, iyr, iyr1, iyr2 - INTEGER :: im, day1, day2, im2 - INTEGER, PARAMETER :: ny=jjm+1 - CHARACTER(len=4) :: cyear - CHARACTER(len=7),DIMENSION(8) :: name_aero - REAL :: var_1(iim, jjm+1, klev, 12) - REAL, DIMENSION(klon,klev) :: massvar - REAL, DIMENSION(iim, jjm+1, klev, 12) :: var_2 ! The massvar distributions - REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: var ! VAR in right dimension - REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: var_out -!$OMP THREADPRIVATE(var,var_out) - - LOGICAL :: lnewday - LOGICAL, PARAMETER :: lonlyone=.FALSE. - LOGICAL,SAVE :: first2=.TRUE. -!$OMP THREADPRIVATE(first2) - -! Variable pour definir a partir d'un numero d'aerosol, son nom - name_aero(1) = "SSSSM " - name_aero(2) = "ASSSM " - name_aero(3) = "ASBCM " - name_aero(4) = "ASPOMM " - name_aero(5) = "SO4 " - name_aero(6) = "CIDUSTM" - name_aero(7) = "AIBCM " - name_aero(8) = "AIPOMM " - -!$OMP MASTER - IF (first2) THEN - ALLOCATE( var(klon, klev, 12,8) ) - ALLOCATE( var_out(klon, klev,8)) - first2=.FALSE. - - IF (aer_type /= 'actuel ' .AND. aer_type /= 'preind ' .AND. & - aer_type /= 'scenario') THEN - WRITE(lunout,*)' *** Warning ***' - WRITE(lunout,*)'Option aer_type non prevu pour les aerosols = ', & - aer_type - CALL abort_gcm('readaerosol','Error : aer_type parameter not accepted',1) - ENDIF - ENDIF - - IF (is_mpi_root) THEN - - iday = INT(r_day) - ! Get the year of the run - iyr = iday/360 - ! Get the day of the actual year: - iday = iday-iyr*360 - ! 0.02 is about 0.5/24, namly less than half an hour - lnewday = (r_day-FLOAT(iday).LT.0.02) - - ! --------------------------------------------- - ! All has to be done only, if a new day begins! - ! --------------------------------------------- - - IF (lnewday.OR.first) THEN - - im = iday/30 +1 ! the actual month - ! annee_ref is the initial year (defined in temps.h) - iyr = iyr + annee_ref - - ! Do I have to read new data? (Is this the first day of a year?) - IF (first.OR.iday.EQ.1.) THEN - ! Initialize values - DO it=1,12 - DO k=1,klev - DO i=1,klon - var(i,k,it,id_aero)=0. - ENDDO - ENDDO - ENDDO - - - IF (aer_type == 'actuel ') then - - cyear='1980' - CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero)) - - ELSE IF (aer_type == 'preind ') THEN - - cyear='.nat' - CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero)) - - ELSE ! aer_type=scenario - - IF (iyr .LT. 1850) THEN - cyear='.nat' - WRITE(lunout,*) 'get_aero iyr=', iyr,' ',cyear - CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero)) - - ELSE IF (iyr .GE. 2100) THEN - cyear='2100' - WRITE(lunout,*) 'get_aero iyr=', iyr,' ',cyear - CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero)) - - ELSE ! Read data from 2 decades - ! a) from actual 10-yr-period - IF (iyr.LT.1900) THEN - iyr1 = 1850 - iyr2 = 1900 - ELSE IF (iyr.GE.1900.AND.iyr.LT.1920) THEN - iyr1 = 1900 - iyr2 = 1920 - ELSE - iyr1 = INT(iyr/10)*10 - iyr2 = INT(1+iyr/10)*10 - ENDIF - - WRITE(cyear,'(I4)') iyr1 - WRITE(lunout,*) 'get_aero iyr=', iyr,' ',cyear - CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero)) - - ! If to read two decades: - IF (.NOT.lonlyone) THEN - - ! b) from the next following one - WRITE(cyear,'(I4)') iyr2 - WRITE(lunout,*) 'get_aero iyr=', iyr,' ',cyear - - CALL get_aero_fromfile(cyear, var_2, name_aero(id_aero)) - - ! Interpolate linarily to the actual year: - DO it=1,12 - DO k=1,klev - DO j=1,jjm - DO i=1,iim - var_1(i,j,k,it) = & - var_1(i,j,k,it) - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1) * & - (var_1(i,j,k,it) - var_2(i,j,k,it)) - ENDDO - ENDDO - ENDDO - ENDDO - - ENDIF ! lonlyone - ENDIF ! iyr .LT. 1850 - ENDIF ! aer_type - - ! Transform the horizontal 2D-field into the physics-field - ! (Also the levels and the latitudes have to be inversed) - - DO it=1,12 - DO k=1, klev - ! a) at the poles, use the zonal mean: - DO i=1,iim - ! North pole - var(1,k,it,id_aero) = & - var(1,k,it,id_aero)+ var_1(i,jjm+1,klev+1-k,it) - ! South pole - var(klon,k,it,id_aero)= & - var(klon,k,it,id_aero)+ var_1(i,1,klev+1-k,it) - ENDDO - var(1,k,it,id_aero) = var(1,k,it,id_aero)/FLOAT(iim) - var(klon,k,it,id_aero)= var(klon,k,it,id_aero)/FLOAT(iim) - - ! b) the values between the poles: - ig=1 - DO j=2,jjm - DO i=1,iim - ig=ig+1 - - IF (ig.GT.klon) THEN - WRITE(lunout,*) 'Attention dans readaerosol', & - name_aero(id_aero), 'ig > klon', ig, klon - ENDIF - - var(ig,k,it,id_aero) = var_1(i,jjm+1+1-j,klev+1-k,it) - - ENDDO - ENDDO - IF (ig.NE.klon-1) THEN - PRINT *,'Error in readaerosol', name_aero(id_aero), 'conversion' - CALL abort_gcm('readaerosol','Error in readaerosol 1',1) - ENDIF - ENDDO ! Loop over k (vertical) - ENDDO ! Loop over it (months) - - ENDIF ! Had to read new data? - - - ! Interpolate to actual day: - IF (iday.LT.im*30-15) THEN - ! in the first half of the month use month before and actual month - im2=im-1 - day2 = im2*30-15 - day1 = im2*30+15 - IF (im2.LE.0) THEN - ! the month is january, thus the month before december - im2=12 - ENDIF - DO k=1,klev - DO i=1,klon - massvar(i,k) = & - var(i,k,im2,id_aero)- FLOAT(iday-day2)/FLOAT(day1-day2) * & - (var(i,k,im2,id_aero) - var(i,k,im,id_aero)) - - IF (massvar(i,k).LT.0.) THEN - IF (iday-day2.LT.0.) WRITE(lunout,*) 'iday-day2',iday-day2 - IF (var(i,k,im2,id_aero) - var(i,k,im,id_aero).LT.0.) THEN - PRINT*, name_aero(id_aero),'(i,k,im2)- ', & - name_aero(id_aero),'(i,k,im)', & - var(i,k,im2,id_aero) - var(i,k,im,id_aero) - ENDIF - - IF (day1-day2.LT.0.) WRITE(lunout,*) 'day1-day2',day1-day2 - WRITE(lunout,*) 'stop',name_aero(id_aero) - CALL abort_gcm('readaerosol','Error in interpolation 2',1) - ENDIF - ENDDO - ENDDO - ELSE - ! the second half of the month - im2=im+1 - IF (im2.GT.12) THEN - ! the month is december, the following thus january - im2=1 - ENDIF - day2 = im*30-15 - day1 = im*30+15 - DO k=1,klev - DO i=1,klon - massvar(i,k) = & - var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * & - (var(i,k,im2,id_aero) - var(i,k,im,id_aero)) - - IF (massvar(i,k).LT.0.) THEN - IF (iday-day2.LT.0.) & - WRITE(lunout,*) 'iday-day2',iday-day2 - IF (var(i,k,im2,id_aero) - var(i,k,im,id_aero).LT.0.) THEN - WRITE(lunout,*) name_aero(id_aero),'(i,k,im2)-', & - name_aero(id_aero),'(i,k,im)', & - var(i,k,im2,id_aero) - var(i,k,im,id_aero) - ENDIF - IF (day1-day2.LT.0.) & - WRITE(lunout,*) 'day1-day2',day1-day2 - WRITE(lunout,*) 'stop',name_aero(id_aero) - CALL abort_gcm('readaerosol','Error in interpolation 3',1) - ENDIF - ENDDO - ENDDO - ENDIF - - - ! The massvar concentration [molec cm-3] is read in. - ! Convert it into mass [ug VAR/m3] - ! masse_var in [g/mol], n_avogadro in [molec/mol] - ! The sulfate mass [ug VAR/m3] is read in. - DO k=1,klev - DO i=1,klon - var_out(i,k,id_aero) = massvar(i,k) - IF (var_out(i,k,id_aero).LT.0) THEN - PRINT*, 'Attention il y a un probleme dans readaerosol', & - name_aero(id_aero), 'la masse est negative' - CALL abort_gcm('readaerosol','Error in readaerosol 4',1) - ENDIF - ENDDO - ENDDO - ELSE ! if no new day, use old data: - DO k=1,klev - DO i=1,klon - massvar(i,k) = var_out(i,k,id_aero) - IF (var_out(i,k,id_aero).LT.0) THEN - PRINT*, 'Attention il y a un probleme dans readaerosol', & - name_aero(id_aero), 'la masse est negative' - CALL abort_gcm('readaerosol','Error in readaerosol 5',1) - ENDIF - ENDDO - ENDDO - - - ENDIF ! Did I have to do anything (was it a new day?) - - ENDIF ! phy_rank==0 - -!$OMP END MASTER - CALL Scatter(massvar,massvar_p) - -END SUBROUTINE readaerosol - - - - - -!----------------------------------------------------------------------------- -! Read in /calculate pre-industrial values of sulfate -!----------------------------------------------------------------------------- - -SUBROUTINE readaerosol_preind (id_aero, r_day, first, pi_massvar_p) - - USE dimphy, ONLY : klev - USE mod_grid_phy_lmdz, klon=>klon_glo - USE mod_phys_lmdz_para - IMPLICIT NONE - - ! Content: - ! -------- - ! This routine reads in monthly mean values of massvar aerosols and - ! returns a linearly interpolated dayly-mean field. - ! - ! It does so for the preindustriel values of the massvar, to a large part - ! analogous to the routine readmassvar above. - ! - ! Only Pb: Variables must be saved and don t have to be overwritten! - ! - ! Author: - ! ------- - ! Celine Deandreis & Anne Cozic (LSCE) 2009 - ! Johannes Quaas (quaas@lmd.jussieu.fr) - ! 26/06/01 - ! - ! Modifications: - ! -------------- - ! see above - - - ! Include-files: - ! -------------- - - INCLUDE "YOMCST.h" - INCLUDE "chem.h" - INCLUDE "dimensions.h" - INCLUDE "temps.h" - INCLUDE "iniprint.h" - - ! Input: - ! ------ - INTEGER, INTENT(in) :: id_aero - REAL, INTENT(in) :: r_day ! Day of integration - LOGICAL, INTENT(in) :: first ! First timestep (and therefore initialization necessary) - - - ! - ! Output: - ! ------- - REAL, DIMENSION(klon_omp,klev), INTENT(out) :: pi_massvar_p ! Number conc. massvar (monthly mean data, from file) - - - ! - ! Local Variables: - ! ---------------- - INTEGER :: i, ig, k, it - INTEGER :: j, iday, iyr - INTEGER, PARAMETER :: ny=jjm+1 - INTEGER :: im, day1, day2, im2 - - REAL, DIMENSION(iim, jjm+1, klev, 12) :: pi_var_1 - REAL, DIMENSION(klon,klev) :: pi_massvar - REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: pi_var ! VAR in right dimension - REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: pi_var_out -!$OMP THREADPRIVATE(pi_var,pi_var_out) - - CHARACTER(len=4) :: cyear - CHARACTER(len=7),DIMENSION(8) :: name_aero - LOGICAL :: lnewday - LOGICAL,SAVE :: first2=.TRUE. -!$OMP THREADPRIVATE(first2) - - -! Variable pour definir a partir d'un numero d'aerosol, son nom - name_aero(1) = "SSSSM " - name_aero(2) = "ASSSM " - name_aero(3) = "ASBCM " - name_aero(4) = "ASPOMM " - name_aero(5) = "SO4 " - name_aero(6) = "CIDUSTM" - name_aero(7) = "AIBCM " - name_aero(8) = "AIPOMM " - - -!$OMP MASTER - IF (first2) THEN - ALLOCATE( pi_var(klon, klev, 12,8) ) - ALLOCATE( pi_var_out(klon, klev,8)) - first2=.FALSE. - ENDIF - - IF (is_mpi_root) THEN - iday = INT(r_day) - ! Get the year of the run - iyr = iday/360 - ! Get the day of the actual year: - iday = iday-iyr*360 - ! 0.02 is about 0.5/24, namly less than half an hour - lnewday = (r_day-FLOAT(iday).LT.0.02) - - ! --------------------------------------------- - ! All has to be done only, if a new day begins! - ! --------------------------------------------- - - IF (lnewday.OR.first) THEN - - im = iday/30 +1 ! the actual month - ! annee_ref is the initial year (defined in temps.h) - iyr = iyr + annee_ref - - IF (first) THEN - cyear='.nat' - CALL get_aero_fromfile(cyear,pi_var_1, name_aero(id_aero)) - - ! Transform the horizontal 2D-field into the physics-field - ! (Also the levels and the latitudes have to be inversed) - - ! Initialize field - DO it=1,12 - DO k=1,klev - DO i=1,klon - pi_var(i,k,it,id_aero)=0. - ENDDO - ENDDO - ENDDO - - WRITE(lunout,*) 'preind: finished reading', FLOAT(iim) - DO it=1,12 - DO k=1, klev - ! a) at the poles, use the zonal mean: - DO i=1,iim - ! North pole - pi_var(1,k,it,id_aero) = & - pi_var(1,k,it,id_aero) + pi_var_1(i,jjm+1,klev+1-k,it) - ! South pole - pi_var(klon,k,it,id_aero) = & - pi_var(klon,k,it,id_aero) + pi_var_1(i,1,klev+1-k,it) - ENDDO - pi_var(1,k,it,id_aero) = pi_var(1,k,it,id_aero)/FLOAT(iim) - pi_var(klon,k,it,id_aero) = pi_var(klon,k,it,id_aero)/FLOAT(iim) - - ! b) the values between the poles: - ig=1 - DO j=2,jjm - DO i=1,iim - ig=ig+1 - IF (ig.GT.klon) THEN - WRITE(lunout,*) 'Attention dans readaerosol_preind', & - name_aero(id_aero), 'ig > klon', ig, klon - ENDIF - pi_var(ig,k,it,id_aero) = & - pi_var_1(i,jjm+1+1-j,klev+1-k,it) - ENDDO - ENDDO - IF (ig.NE.klon-1) THEN - WRITE(lunout,*) 'Error in readaerosol_preind (', name_aero(id_aero), ' conversion)' - CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 1',1) - ENDIF - ENDDO ! Loop over k (vertical) - ENDDO ! Loop over it (months) - - ENDIF ! Had to read new data? - - - ! Interpolate to actual day: - IF (iday.LT.im*30-15) THEN - ! in the first half of the month use month before and actual month - im2=im-1 - day1 = im2*30+15 - day2 = im2*30-15 - IF (im2.LE.0) THEN - ! the month is january, thus the month before december - im2=12 - ENDIF - DO k=1,klev - DO i=1,klon - pi_massvar(i,k) = & - pi_var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * & - (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero)) - - IF (pi_massvar(i,k).LT.0.) THEN - IF (iday-day2.LT.0.) WRITE(lunout,*) 'iday-day2',iday-day2 - IF (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero).LT.0.) THEN - WRITE(lunout,*) 'pi_',name_aero(id_aero),'(i,k,im2) - pi_', & - name_aero(id_aero),'(i,k,im)', & - pi_var(i,k,im2,id_aero) -pi_var(i,k,im,id_aero) - ENDIF - IF (day1-day2.LT.0.)WRITE(lunout,*) 'day1-day2',day1-day2 - PRINT *, 'pi_',name_aero(id_aero) - CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 2',1) - ENDIF - ENDDO - ENDDO - ELSE - - ! the second half of the month - im2=im+1 - day1 = im*30+15 - IF (im2.GT.12) THEN - ! the month is december, the following thus january - im2=1 - ENDIF - day2 = im*30-15 - - DO k=1,klev - DO i=1,klon - pi_massvar(i,k) = & - pi_var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * & - (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero)) - IF (pi_massvar(i,k).LT.0.) THEN - IF (iday-day2.LT.0.) WRITE(lunout,*) 'iday-day2',iday-day2 - IF (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero).LT.0.) THEN - WRITE(lunout,*) 'pi_', name_aero(id_aero),'(i,k,im2) - pi_',& - name_aero(id_aero), '(i,k,im)',& - pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero) - ENDIF - IF (day1-day2.LT.0.) & - WRITE(lunout,*) 'day1-day2',day1-day2 - PRINT *, 'pi_',name_aero(id_aero) - CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 3',1) - ENDIF - ENDDO - ENDDO - ENDIF - - - ! The massvar concentration [molec cm-3] is read in. - ! Convert it into mass [ug VAR/m3] - ! masse_var in [g/mol], n_avogadro in [molec/mol] - DO k=1,klev - DO i=1,klon - pi_var_out(i,k,id_aero) = pi_massvar(i,k) - ENDDO - ENDDO - - ELSE ! If no new day, use old data: - DO k=1,klev - DO i=1,klon - pi_massvar(i,k) = pi_var_out(i,k,id_aero) - ENDDO - ENDDO - ENDIF ! Was this the beginning of a new day? - ENDIF ! is_mpi_root - -!$OMP END MASTER - CALL Scatter(pi_massvar,pi_massvar_p) - - END SUBROUTINE readaerosol_preind - - - -!----------------------------------------------------------------------------- -! Routine for reading VAR data from files -!----------------------------------------------------------------------------- - - -SUBROUTINE get_aero_fromfile (cyr, var, varname) - USE dimphy - IMPLICIT NONE - - INCLUDE "netcdf.inc" - INCLUDE "dimensions.h" - INCLUDE "iniprint.h" - - CHARACTER(len=4), INTENT(in) :: cyr - REAL, DIMENSION(iim, jjm+1, klev, 12), INTENT(out) :: var - CHARACTER(len=*), INTENT(in) :: varname - - - CHARACTER(len=30) :: fname - CHARACTER(len=30) :: cvar - INTEGER, DIMENSION(3) :: START, COUNT - INTEGER :: STATUS, NCID, VARID - INTEGER :: imth, i, j, k - INTEGER, PARAMETER :: ny=jjm+1 - REAL, DIMENSION(iim, ny, klev) :: varmth -! -! -! - fname = TRIM(varname)//'.run'//cyr//'.cdf' - - WRITE(lunout,*) 'reading ', TRIM(fname) - STATUS = NF_OPEN (TRIM(fname), NF_NOWRITE, NCID) - IF (STATUS .NE. NF_NOERR) THEN - WRITE(lunout,*) 'err in open ',status - CALL abort_gcm('get_aero_fromfile','Error in opening file',1) - ENDIF - DO imth=1, 12 - IF (imth.EQ.1) THEN - cvar=TRIM(varname)//'JAN' - ELSEIF (imth.EQ.2) THEN - cvar=TRIM(varname)//'FEB' - ELSEIF (imth.EQ.3) THEN - cvar=TRIM(varname)//'MAR' - ELSEIF (imth.EQ.4) THEN - cvar=TRIM(varname)//'APR' - ELSEIF (imth.EQ.5) THEN - cvar=TRIM(varname)//'MAY' - ELSEIF (imth.EQ.6) THEN - cvar=TRIM(varname)//'JUN' - ELSEIF (imth.EQ.7) THEN - cvar=TRIM(varname)//'JUL' - ELSEIF (imth.EQ.8) THEN - cvar=TRIM(varname)//'AUG' - ELSEIF (imth.EQ.9) THEN - cvar=TRIM(varname)//'SEP' - ELSEIF (imth.EQ.10) THEN - cvar=TRIM(varname)//'OCT' - ELSEIF (imth.EQ.11) THEN - cvar=TRIM(varname)//'NOV' - ELSEIF (imth.EQ.12) THEN - cvar=TRIM(varname)//'DEC' - ENDIF - start(1)=1 - start(2)=1 - start(3)=1 - COUNT(1)=iim - COUNT(2)=ny - COUNT(3)=klev - STATUS = NF_INQ_VARID (NCID, TRIM(cvar), VARID) - WRITE(lunout,*) ncid,imth,TRIM(cvar), varid - - IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in read ',status - -#ifdef NC_DOUBLE - status = NF_GET_VARA_DOUBLE(NCID, VARID, START, COUNT,varmth) -#else - status = NF_GET_VARA_REAL(NCID, VARID, START, COUNT, varmth) -#endif - IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in read data',status - - DO k=1,klev - DO j=1,jjm+1 - DO i=1,iim - IF (varmth(i,j,k).LT.0.) THEN - WRITE(lunout,*) 'this is shit' - WRITE(lunout,*) varname,'(',i,j,k,') =',varmth(i,j,k) - ENDIF - var(i,j,k,imth)=varmth(i,j,k) - ENDDO - ENDDO - ENDDO - ENDDO - - STATUS = NF_CLOSE(NCID) - IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in closing file',status - - -END SUBROUTINE get_aero_fromfile - - - - - - - - - - - - - - - - +! 6) Distribute to all processes +! Scatter global field(klon_glo) to local process domain(klon) +! and distribute klev_src to all processes +!**************************************************************************************** + + ! Distribute klev_src + CALL bcast(klev_src) + + ! Allocate and distribute pt_ap and pt_b + IF (.NOT. ASSOCIATED(pt_ap)) THEN ! if pt_ap is allocated also pt_b is allocated + ALLOCATE(pt_ap(klev_src), pt_b(klev_src), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 4',1) + END IF + CALL bcast(pt_ap) + CALL bcast(pt_b) + + ! Allocate space for output pointer variable at local process + IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year) + ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 5',1) + + ! Scatter global field to local domain at local process + CALL scatter(varyear_glo1D, pt_year) + CALL scatter(psurf_glo1D, psurf_out) + CALL scatter(load_glo1D, load_out) + +! 7) Test for negative values +!**************************************************************************************** + IF (MINVAL(pt_year) < 0.) THEN + WRITE(lunout,*) 'Warning! Negative values read from file :', fname + END IF + + END SUBROUTINE get_aero_fromfile + + + SUBROUTINE check_err(status) + USE netcdf + IMPLICIT NONE + + INCLUDE "iniprint.h" + INTEGER, INTENT (IN) :: status + + IF (status /= NF90_NOERR) THEN + WRITE(lunout,*) 'Error in get_aero_fromfile ',status + CALL abort_gcm('get_aero_fromfile',trim(nf90_strerror(status)),1) + END IF + + END SUBROUTINE check_err + + +END MODULE readaerosol_mod Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol_interp.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol_interp.F90 (revision 1179) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol_interp.F90 (revision 1179) @@ -0,0 +1,453 @@ +! $Id$ +! +SUBROUTINE readaerosol_interp(id_aero, r_day, first, pplay, paprs, mass_out, pi_mass_out) +! +! This routine will return the mass concentration at actual day(mass_out) and +! the pre-industrial values(pi_mass_out) for aerosol corresponding to "id_aero". +! The mass concentrations for all aerosols are saved in this routine but each +! call to this routine only treats the aerosol "id_aero". +! +! 1) Read in data for the whole year, only at first time step +! 2) Interpolate to the actual day, only at new day +! 3) Interpolate to the model vertical grid (target grid), only at new day +! 4) Test for negative mass values + + USE dimphy, ONLY : klev,klon + USE mod_phys_lmdz_para, ONLY : mpi_rank + USE readaerosol_mod + USE write_field_phy + + IMPLICIT NONE + + INCLUDE "YOMCST.h" + INCLUDE "chem.h" + INCLUDE "temps.h" + INCLUDE "clesphys.h" + INCLUDE "iniprint.h" + INCLUDE "dimensions.h" + INCLUDE "comvert.h" +! +! Input: +!**************************************************************************************** + INTEGER, INTENT(IN) :: id_aero! Identity number for the aerosol to treat + REAL, INTENT(IN) :: r_day ! Day of integration + LOGICAL, INTENT(IN) :: first ! First model timestep + REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! pression at model mid-layers + REAL, DIMENSION(klon,klev+1),INTENT(IN):: paprs ! pression between model layers +! +! Output: +!**************************************************************************************** + REAL, INTENT(OUT) :: mass_out(klon,klev) ! Mass of aerosol (monthly mean data,from file) [ug AIBCM/m3] + REAL, INTENT(OUT) :: pi_mass_out(klon,klev) ! Mass of preindustrial aerosol (monthly mean data,from file) [ug AIBCM/m3] +! +! Local Variables: +!**************************************************************************************** + INTEGER :: i, k, ierr + INTEGER :: iday, iyr + INTEGER :: im, day1, day2, im2 + INTEGER :: pi_klev_src ! Only for testing purpose + INTEGER, SAVE :: klev_src ! Number of vertical levles in source field +!$OMP THREADPRIVATE(klev_src) + INTEGER, PARAMETER :: nb_aero=8 + + CHARACTER(len=7),DIMENSION(nb_aero) :: name_aero + + REAL, DIMENSION(klon) :: psurf_day, pi_psurf_day + REAL, DIMENSION(klon) :: load_src, pi_load_src ! Mass load at source grid + REAL, DIMENSION(klon) :: load_tgt, load_tgt_test + REAL, DIMENSION(klon,klev) :: delp ! pressure difference in each model layer + + REAL, ALLOCATABLE, DIMENSION(:,:) :: pplay_src ! pression mid-layer at source levels + REAL, ALLOCATABLE, DIMENSION(:,:) :: tmp1, tmp2 ! Temporary variables + REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: var_year ! VAR in right dimension for the total year + REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE :: pi_var_year ! pre-industrial VAR, -"- +!$OMP THREADPRIVATE(var_year,pi_var_year) + REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: var_day ! VAR interpolated to the actual day and model grid + REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE :: pi_var_day ! pre-industrial VAR, -"- +!$OMP THREADPRIVATE(var_day,pi_var_day) + REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: psurf_year, pi_psurf_year ! surface pressure for the total year +!$OMP THREADPRIVATE(psurf_year, pi_psurf_year) + REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: load_year, pi_load_year ! load in the column for the total year +!$OMP THREADPRIVATE(load_year, pi_load_year) + + REAL, DIMENSION(:,:,:), POINTER :: pt_tmp ! Pointer allocated in readaerosol + REAL, POINTER, DIMENSION(:), SAVE :: pt_ap, pt_b ! Pointer for describing the vertical levels +!$OMP THREADPRIVATE(pt_ap, pt_b) + + LOGICAL :: lnewday ! Indicates if first time step at a new day + LOGICAL,SAVE :: vert_interp ! Indicates if vertical interpolation will be done + LOGICAL,SAVE :: debug=.TRUE. ! Debugging in this subroutine +!$OMP THREADPRIVATE(vert_interp, debug) + + +!**************************************************************************************** +! Initialization +! +!**************************************************************************************** +! Variable containing aerosols name + name_aero(1) = "SSSSM " + name_aero(2) = "ASSSM " + name_aero(3) = "ASBCM " + name_aero(4) = "ASPOMM " + name_aero(5) = "SO4 " + name_aero(6) = "CIDUSTM" + name_aero(7) = "AIBCM " + name_aero(8) = "AIPOMM " + +! Calculation to find if it is a new day + iday = INT(r_day) + iyr = iday/360 + iday = iday-iyr*360 ! day of the actual year + iyr = iyr + annee_ref ! year of the run + im = iday/30 +1 ! the actual month + + ! 0.02 is about 0.5/24, namly less than half an hour + lnewday = (r_day-FLOAT(iday) < 0.02) + + IF (.NOT. ALLOCATED(var_day)) THEN + ALLOCATE( var_day(klon, klev, nb_aero), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 1',1) + ALLOCATE( pi_var_day(klon, klev, nb_aero), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 2',1) + + ALLOCATE( psurf_year(klon, 12, nb_aero), pi_psurf_year(klon, 12, nb_aero), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 3',1) + + ALLOCATE( load_year(klon, 12, nb_aero), pi_load_year(klon, 12, nb_aero), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 4',1) + + lnewday=.TRUE. + + NULLIFY(pt_ap) + NULLIFY(pt_b) + END IF + +!**************************************************************************************** +! 1) Read in data : corresponding to the actual year and preindustrial data. +! Only for the first day of the year. +! +!**************************************************************************************** + IF ( (first .OR. iday==1.) .AND. lnewday ) THEN + NULLIFY(pt_tmp) + + ! Reading values corresponding to the closest year taking into count the choice of aer_type. + ! For aer_type=scenario interpolation between 2 data sets is done in readaerosol. + CALL readaerosol(name_aero(id_aero), aer_type, iyr, klev_src, pt_ap, pt_b, pt_tmp, & + psurf_year(:,:,id_aero), load_year(:,:,id_aero)) + IF (.NOT. ALLOCATED(var_year)) THEN + ALLOCATE(var_year(klon, klev_src, 12, nb_aero), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 5',1) + END IF + var_year(:,:,:,id_aero) = pt_tmp(:,:,:) + + ! Reading values corresponding to the preindustrial concentrations. + CALL readaerosol(name_aero(id_aero), 'preind', iyr, pi_klev_src, pt_ap, pt_b, pt_tmp, & + pi_psurf_year(:,:,id_aero), pi_load_year(:,:,id_aero)) + + ! klev_src must be the same in both files. + ! Also supposing pt_ap and pt_b to be the same in the 2 files without testing. + IF (pi_klev_src /= klev_src) THEN + WRITE(lunout,*) 'Error! All forcing files for the same aerosol must have the same vertical dimension' + WRITE(lunout,*) 'Aerosol : ', name_aero(id_aero) + CALL abort_gcm('readaerosol_interp','Differnt vertical axes in aerosol forcing files',1) + END IF + + IF (.NOT. ALLOCATED(pi_var_year)) THEN + ALLOCATE(pi_var_year(klon, klev_src, 12, nb_aero), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 6',1) + END IF + pi_var_year(:,:,:,id_aero) = pt_tmp(:,:,:) + + IF (debug) THEN + CALL writefield_phy('psurf_src',psurf_year(:,:,id_aero),1) + CALL writefield_phy('pi_psurf_src',pi_psurf_year(:,:,id_aero),1) + CALL writefield_phy('load_year_src',load_year(:,:,id_aero),1) + CALL writefield_phy('pi_load_year_src',pi_load_year(:,:,id_aero),1) + END IF + + ! Pointer no more useful, deallocate. + DEALLOCATE(pt_tmp) + + ! Test if vertical interpolation will be needed. + IF (psurf_year(1,1,id_aero)==not_valid .OR. pi_psurf_year(1,1,id_aero)==not_valid ) THEN + ! Pressure=not_valid indicates old file format, see module readaerosol + vert_interp = .FALSE. + + ! If old file format, both psurf_year and pi_psurf_year must be not_valid + IF ( psurf_year(1,1,id_aero) /= pi_psurf_year(1,1,id_aero) ) THEN + WRITE(lunout,*) 'Warning! All forcing files for the same aerosol must have the same structure' + CALL abort_gcm('readaerosol_interp', 'The aerosol files have not the same format',1) + END IF + + IF (klev /= klev_src) THEN + WRITE(lunout,*) 'Old format of aerosol file do not allowed vertical interpolation' + CALL abort_gcm('readaerosol_interp', 'Old aerosol file not possible',1) + END IF + + ELSE + vert_interp = .TRUE. + END IF + + END IF ! IF ( (first .OR. iday==1.) .AND. lnewday ) THEN + +!**************************************************************************************** +! - 2) Interpolate to the actual day. +! - 3) Interpolate to the model vertical grid. +! +!**************************************************************************************** + + IF (lnewday) THEN ! only if new day +!**************************************************************************************** +! 2) Interpolate to the actual day +! +!**************************************************************************************** + ! Find which months and days to use for time interpolation + IF (iday < im*30-15) THEN + ! in the first half of the month use month before and actual month + im2=im-1 + day2 = im2*30-15 + day1 = im2*30+15 + IF (im2 <= 0) THEN + ! the month is january, thus the month before december + im2=12 + END IF + ELSE + ! the second half of the month + im2=im+1 + IF (im2 > 12) THEN + ! the month is december, the following thus january + im2=1 + ENDIF + day2 = im*30-15 + day1 = im*30+15 + END IF + + ! Time interpolation, still on vertical source grid + ALLOCATE(tmp1(klon,klev_src), tmp2(klon,klev_src),stat=ierr) + IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 7',1) + + ALLOCATE(pplay_src(klon,klev_src), stat=ierr) + IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 8',1) + + + DO k=1,klev_src + DO i=1,klon + tmp1(i,k) = & + var_year(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * & + (var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero)) + + tmp2(i,k) = & + pi_var_year(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * & + (pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero)) + END DO + END DO + + ! Time interpolation for pressure at surface, still on vertical source grid + DO i=1,klon + psurf_day(i) = & + psurf_year(i,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * & + (psurf_year(i,im2,id_aero) - psurf_year(i,im,id_aero)) + + pi_psurf_day(i) = & + pi_psurf_year(i,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * & + (pi_psurf_year(i,im2,id_aero) - pi_psurf_year(i,im,id_aero)) + END DO + + ! Time interpolation for the load, still on vertical source grid + DO i=1,klon + load_src(i) = & + load_year(i,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * & + (load_year(i,im2,id_aero) - load_year(i,im,id_aero)) + + pi_load_src(i) = & + pi_load_year(i,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * & + (pi_load_year(i,im2,id_aero) - pi_load_year(i,im,id_aero)) + END DO + +!**************************************************************************************** +! 3) Interpolate to the model vertical grid (target grid) +! +!**************************************************************************************** + + IF (vert_interp) THEN + + ! - Interpolate variable tmp1 (on source grid) to var_day (on target grid) + !******************************************************************************** + ! a) calculate pression at vertical levels for the source grid using the + ! hybrid-sigma coordinates ap and b and the surface pressure, variables from file. + DO k = 1, klev_src + DO i = 1, klon + pplay_src(i,k)= pt_ap(k) + pt_b(k)*psurf_day(i) + END DO + END DO + + IF (debug) THEN + CALL writefield_phy('psurf_day_src',psurf_day(:),1) + CALL writefield_phy('pplay_src',pplay_src(:,:),klev_src) + CALL writefield_phy('pplay',pplay(:,:),klev) + CALL writefield_phy('day_src',tmp1,klev_src) + CALL writefield_phy('pi_day_src',tmp2,klev_src) + END IF + + ! b) vertical interpolation on pressure leveles + CALL pres2lev(tmp1(:,:), var_day(:,:,id_aero), klev_src, klev, pplay_src, pplay, & + 1, klon, .FALSE.) + + IF (debug) CALL writefield_phy('day_tgt',var_day(:,:,id_aero),klev) + + ! c) adjust to conserve total aerosol mass load in the vertical pillar + ! Calculate the load in the actual pillar and compare with the load + ! read from aerosol file. + + ! Find the pressure difference in each model layer + DO k = 1, klev + DO i = 1, klon + delp(i,k) = paprs(i,k) - paprs (i,k+1) + END DO + END DO + + ! Find the mass load in the actual pillar, on target grid + load_tgt(:) = 0. + DO k= 1, klev + DO i = 1, klon + load_tgt(i) = load_tgt(i) + 1/RG * var_day(i,k,id_aero)*delp(i,k) + END DO + END DO + + ! Adjust, uniform + DO k = 1, klev + DO i = 1, klon + var_day(i,k,id_aero) = var_day(i,k,id_aero)*load_src(i)/load_tgt(i) + END DO + END DO + + IF (debug) THEN + load_tgt_test(:) = 0. + DO k= 1, klev + DO i = 1, klon + load_tgt_test(i) = load_tgt_test(i) + 1/RG * var_day(i,k,id_aero)*delp(i,k) + END DO + END DO + + CALL writefield_phy('day_tgt2',var_day(:,:,id_aero),klev) + CALL writefield_phy('load_tgt',load_tgt(:),1) + CALL writefield_phy('load_tgt_test',load_tgt_test(:),1) + CALL writefield_phy('load_src',load_src(:),1) + END IF + + ! - Interpolate variable tmp2 (source grid) to pi_var_day (target grid) + !******************************************************************************** + ! a) calculate pression at vertical levels at source grid + DO k = 1, klev_src + DO i = 1, klon + pplay_src(i,k)= pt_ap(k) + pt_b(k)*pi_psurf_day(i) + END DO + END DO + + IF (debug) THEN + CALL writefield_phy('pi_psurf_day_src',pi_psurf_day(:),1) + CALL writefield_phy('pi_pplay_src',pplay_src(:,:),klev_src) + END IF + + ! b) vertical interpolation on pressure leveles + CALL pres2lev(tmp2(:,:), pi_var_day(:,:,id_aero), klev_src, klev, pplay_src, pplay, & + 1, klon, .FALSE.) + + IF (debug) CALL writefield_phy('pi_day_tgt',pi_var_day(:,:,id_aero),klev) + + ! c) adjust to conserve total aerosol mass load in the vertical pillar + ! Calculate the load in the actual pillar and compare with the load + ! read from aerosol file. + + ! Find the load in the actual pillar, on target grid + load_tgt(:) = 0. + DO k = 1, klev + DO i = 1, klon + load_tgt(i) = load_tgt(i) + 1/RG * pi_var_day(i,k,id_aero)*delp(i,k) + END DO + END DO + + DO k = 1, klev + DO i = 1, klon + pi_var_day(i,k,id_aero) = pi_var_day(i,k,id_aero)*pi_load_src(i)/load_tgt(i) + END DO + END DO + + IF (debug) THEN + load_tgt_test(:) = 0. + DO k = 1, klev + DO i = 1, klon + load_tgt_test(i) = load_tgt_test(i) + 1/RG * pi_var_day(i,k,id_aero)*delp(i,k) + END DO + END DO + CALL writefield_phy('pi_day_tgt2',pi_var_day(:,:,id_aero),klev) + CALL writefield_phy('pi_load_tgt',load_tgt(:),1) + CALL writefield_phy('pi_load_tgt_test',load_tgt_test(:),1) + CALL writefield_phy('pi_load_src',pi_load_src(:),1) + END IF + + + ELSE ! No vertical interpolation done + + var_day(:,:,id_aero) = tmp1(:,:) + pi_var_day(:,:,id_aero) = tmp2(:,:) + + END IF ! vert_interp + + + ! Deallocation + DEALLOCATE(tmp1, tmp2, pplay_src, stat=ierr) + +!**************************************************************************************** +! 4) Test for negative mass values +! +!**************************************************************************************** + IF (MINVAL(var_day(:,:,id_aero)) < 0.) THEN + DO k=1,klev + DO i=1,klon + ! Test for var_day + IF (var_day(i,k,id_aero) < 0.) THEN + IF (iday-day2 < 0.) WRITE(lunout,*) 'iday-day2=',iday-day2 + IF (var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero) < 0.) THEN + WRITE(lunout,*) trim(name_aero(id_aero)),'(i,k,im2)-', & + trim(name_aero(id_aero)),'(i,k,im)=', & + var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero) + END IF + + WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero) + CALL abort_gcm('readaerosol_interp','Error in interpolation 1',1) + END IF + END DO + END DO + END IF + + IF (MINVAL(pi_var_day(:,:,id_aero)) < 0. ) THEN + DO k=1, klev + DO i=1,klon + ! Test for pi_var_day + IF (pi_var_day(i,k,id_aero) < 0.) THEN + IF (iday-day2 < 0.) WRITE(lunout,*) 'iday-day2=',iday-day2 + IF (pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero) < 0.) THEN + WRITE(lunout,*) trim(name_aero(id_aero)),'(i,k,im2)-', & + trim(name_aero(id_aero)),'(i,k,im)=', & + pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero) + END IF + + WRITE(lunout,*) 'stop for aerosol : ',name_aero(id_aero) + CALL abort_gcm('readaerosol_interp','Error in interpolation 2',1) + END IF + END DO + END DO + END IF + + END IF ! lnewday + +!**************************************************************************************** +! Copy output from saved variables +! +!**************************************************************************************** + + mass_out(:,:) = var_day(:,:,id_aero) + pi_mass_out(:,:) = pi_var_day(:,:,id_aero) + +END SUBROUTINE readaerosol_interp Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol_optic.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol_optic.F90 (revision 1179) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol_optic.F90 (revision 1179) @@ -0,0 +1,152 @@ +! $Id$ +! +SUBROUTINE readaerosol_optic(debut, new_aod, flag_aerosol, rjourvrai, pdtphys, & + pplay, paprs, t_seri, rhcl, & + mass_ins_aero, mass_ins_aero_pi, & + tau_aero, piz_aero, cg_aero ) + +! This routine will : +! 1) recevie the aerosols(already read and interpolated) corresponding to flag_aerosol +! 2) calculate the optical properties for the aerosols +! + + USE dimphy + IMPLICIT NONE + +! Input arguments +!**************************************************************************************** + LOGICAL, INTENT(IN) :: debut + LOGICAL, INTENT(IN) :: new_aod + INTEGER, INTENT(IN) :: flag_aerosol + REAL, INTENT(IN) :: rjourvrai + REAL, INTENT(IN) :: pdtphys + REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay + REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs + REAL, DIMENSION(klon,klev), INTENT(IN) :: t_seri + REAL, DIMENSION(klon,klev), INTENT(IN) :: rhcl ! humidite relative ciel clair + +! Output arguments +!**************************************************************************************** + REAL, DIMENSION(klon,klev), INTENT(OUT) :: mass_ins_aero ! Total mass for all indissoluble aerosols + REAL, DIMENSION(klon,klev), INTENT(OUT) :: mass_ins_aero_pi ! -"- preindustrial values + REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: tau_aero ! Aerosol optical thickness + REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: piz_aero ! Single scattering albedo aerosol + REAL, DIMENSION(klon,klev,9,2), INTENT(OUT) :: cg_aero ! asymmetry parameter aerosol + +! Local variables +!**************************************************************************************** + REAL, DIMENSION(klon) :: aerindex ! POLDER aerosol index + REAL, DIMENSION(klon,10) :: fractnat_allaer + REAL, DIMENSION(klon,klev) :: sulfate ! SO4 aerosol concentration [ug/m3] + REAL, DIMENSION(klon,klev) :: bcsol ! BC soluble concentration [ug/m3] + REAL, DIMENSION(klon,klev) :: bcins ! BC insoluble concentration [ug/m3] + REAL, DIMENSION(klon,klev) :: pomsol ! POM soluble concentration [ug/m3] + REAL, DIMENSION(klon,klev) :: pomins ! POM insoluble concentration [ug/m3] + REAL, DIMENSION(klon,klev) :: sulfate_pi + REAL, DIMENSION(klon,klev) :: bcsol_pi + REAL, DIMENSION(klon,klev) :: bcins_pi + REAL, DIMENSION(klon,klev) :: pomsol_pi + REAL, DIMENSION(klon,klev) :: pomins_pi + REAL, DIMENSION(klon,klev) :: pdel + REAL, DIMENSION(klon,klev,8) :: m_allaer + + INTEGER :: k, i + +!**************************************************************************************** +! 1) Get aerosol mass +! +!**************************************************************************************** +! Read and interpolate sulfate + IF ( flag_aerosol .EQ. 1 .OR. & + flag_aerosol .EQ. 4 .OR. & + flag_aerosol .EQ. 6 ) THEN + + CALL readaerosol_interp(5, rjourvrai, debut, pplay, paprs, sulfate, sulfate_pi) + ELSE + sulfate(:,:) = 0. ; sulfate_pi(:,:) = 0. + END IF + +! Read and interpolate bcsol and bcins + IF ( flag_aerosol .EQ. 2 .OR. & + flag_aerosol .EQ. 4 .OR. & + flag_aerosol .EQ. 5 ) THEN + + ! Get bc aerosol distribution + CALL readaerosol_interp(3, rjourvrai, debut, pplay, paprs, bcsol, bcsol_pi ) + CALL readaerosol_interp(7, rjourvrai, debut, pplay, paprs, bcins, bcins_pi ) + ELSE + bcsol(:,:) = 0. ; bcsol_pi(:,:) = 0. + bcins(:,:) = 0. ; bcins_pi(:,:) = 0. + END IF + + +! Read and interpolate pomsol and pomins + IF ( flag_aerosol .EQ. 3 .OR. & + flag_aerosol .EQ. 4 .OR. & + flag_aerosol .EQ. 5 .OR. & + flag_aerosol .EQ. 6 ) THEN + + CALL readaerosol_interp(4, rjourvrai, debut, pplay, paprs, pomsol, pomsol_pi) + CALL readaerosol_interp(8, rjourvrai, debut, pplay, paprs, pomins, pomins_pi) + ELSE + pomsol(:,:) = 0. ; pomsol_pi(:,:) = 0. + pomins(:,:) = 0. ; pomins_pi(:,:) = 0. + END IF + + +! Store all aerosols in one variable +! +! ACo pour couplage aerosol offline 07/04/2009 +! Tableau contenant les masses pour tous les aerosols +! les valeurs a zero seront a remplacer par les bons +! tableaux lorsque les routines de lectures seront +! ajoutees. + m_allaer(:,:,1) = 0. ! SSSSM || CSSSM ! Coarse Soluble Sea Salt Mass + m_allaer(:,:,2) = 0. ! ASSSM + m_allaer(:,:,3) = bcsol(:,:) ! ASBCM + m_allaer(:,:,4) = pomsol(:,:) ! ASPOMM + m_allaer(:,:,5) = sulfate(:,:) ! ASSO4M || CSSO4M + m_allaer(:,:,6) = 0. ! CIDUSTM ! Coarse Insoluble DUST Mass + m_allaer(:,:,7) = bcins(:,:) ! AIBCM + m_allaer(:,:,8) = pomins(:,:) ! AIPOMM + +! +! Calculate the total mass of all indissoluble aersosols +! + mass_ins_aero(:,:) = sulfate(:,:) + bcsol(:,:) + pomsol(:,:) + mass_ins_aero_pi(:,:) = sulfate_pi(:,:) + bcsol_pi(:,:) + pomsol_pi(:,:) + +!**************************************************************************************** +! 2) Calculate optical properties for the aerosols +! +!**************************************************************************************** + DO k = 1, klev + DO i = 1, klon + pdel(i,k) = paprs(i,k) - paprs (i,k+1) + END DO + END DO + + IF (new_aod) THEN + + fractnat_allaer(:,:) = 0. + CALL aeropt_2bands( & + pdel, m_allaer, pdtphys, rhcl, & + tau_aero, piz_aero, cg_aero, & + fractnat_allaer, flag_aerosol, & + pplay, t_seri) + + ! aeropt_5wv only for validation and diagnostics. + ! In this version no diagnostics are set. + ! jg : may be desactivated if no diagnostics added. + CALL aeropt_5wv( & + pdel, m_allaer, & + pdtphys, rhcl, aerindex, & + flag_aerosol, pplay, t_seri) + ELSE + + CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl, & + tau_aero(:,:,5,:), piz_aero(:,:,5,:), cg_aero(:,:,5,:), aerindex) + + END IF + +END SUBROUTINE readaerosol_optic