Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt.F (revision 1149) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt.F (revision 1150) @@ -16,12 +16,12 @@ c Arguments: c - REAL paprs(klon,klev+1) - REAL pplay(klon,klev), t_seri(klon,klev) - REAL msulfate(klon,klev) ! masse sulfate ug SO4/m3 [ug/m^3] - REAL RHcl(klon,klev) ! humidite relative ciel clair - REAL tau_ae(klon,klev,2) ! epaisseur optique aerosol - REAL piz_ae(klon,klev,2) ! single scattering albedo aerosol - REAL cg_ae(klon,klev,2) ! asymmetry parameter aerosol - REAL ai(klon) ! POLDER aerosol index + REAL, INTENT(in) :: paprs(klon,klev+1) + REAL, INTENT(in) :: pplay(klon,klev), t_seri(klon,klev) + REAL, INTENT(in) :: msulfate(klon,klev) ! masse sulfate ug SO4/m3 [ug/m^3] + REAL, INTENT(in) :: RHcl(klon,klev) ! humidite relative ciel clair + REAL, INTENT(out) :: tau_ae(klon,klev,2) ! epaisseur optique aerosol + REAL, INTENT(out) :: piz_ae(klon,klev,2) ! single scattering albedo aerosol + REAL, INTENT(out) :: cg_ae(klon,klev,2) ! asymmetry parameter aerosol + REAL, INTENT(out) :: ai(klon) ! POLDER aerosol index c c Local Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_2bands.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_2bands.F90 (revision 1150) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_2bands.F90 (revision 1150) @@ -0,0 +1,507 @@ +SUBROUTINE AEROPT_2BANDS( & + pdel, m_allaer, delt, RHcl, & + tau_allaer, piz_allaer, & + cg_allaer, fractnat_allaer, & + flag_aerosol, pplay, t_seri) + + USE dimphy + + + ! Yves Balkanski le 12 avril 2006 + ! Celine Deandreis + ! Anne Cozic Avril 2009 + ! a partir d'une sous-routine de Johannes Quaas pour les sulfates + ! + IMPLICIT NONE + + INCLUDE "YOMCST.h" + INCLUDE "iniprint.h" + + ! + ! Input arguments: + ! + REAL, INTENT(in) :: pdel(KLON,KLEV) + REAL, INTENT(in) :: delt + REAL, DIMENSION(klon,klev,8), INTENT(in) :: m_allaer + REAL, INTENT(in) :: RHcl(KLON,KLEV) ! humidite relative ciel clair + REAL, DIMENSION(klon,10), INTENT(in) :: fractnat_allaer + INTEGER, INTENT(in) :: flag_aerosol + REAL, INTENT(in) :: pplay(klon,klev) + REAL, INTENT(in) :: t_seri(klon,klev) + + ! + ! Output arguments: + ! + REAL, DIMENSION(klon,klev,9,2), INTENT(out) :: tau_allaer ! epaisseur optique aerosol + REAL, DIMENSION(klon,klev,9,2), INTENT(out) :: piz_allaer ! single scattering albedo aerosol + REAL, DIMENSION(klon,klev,9,2), INTENT(out) :: cg_allaer ! asymmetry parameter aerosol + + ! + ! Local + ! + REAL, DIMENSION(klon,klev,10,2) :: tau_ae + REAL, DIMENSION(klon,klev,10,2) :: piz_ae + REAL, DIMENSION(klon,klev,10,2) :: cg_ae + LOGICAL :: soluble + INTEGER :: i, k, inu, m, mrfspecies + INTEGER :: spsol, spinsol + INTEGER :: RH_num, nbre_RH, nbsol_compaer, nbinsol_compaer + + PARAMETER (nbre_RH=12) + PARAMETER (nbsol_compaer=5) ! 1- Seasalt AS: 2- Sesalt CS; 3- BC soluble; 4- POM soluble; 5- SO4. + PARAMETER (nbinsol_compaer=3) ! 1- Dust; 2- BC insoluble; 3- POM insoluble + REAL:: RH_tab(nbre_RH) + REAL:: RH_MAX, DELTA, rh + REAL:: tau_ae2b_int(KLON,KLEV,2) ! Intermediate computation of epaisseur optique aerosol + REAL:: piz_ae2b_int(KLON,KLEV,2) ! Intermediate computation of Single scattering albedo + REAL:: cg_ae2b_int(KLON,KLEV,2) ! Intermediate computation of Assymetry parameter + PARAMETER (RH_MAX=95.) + DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./ + REAL :: zrho + REAL :: fac + REAL :: zdp1(klon,klev) + REAL, PARAMETER :: gravit = 9.80616 ! m2/s + INTEGER, ALLOCATABLE, DIMENSION(:) :: aerosol_name + INTEGER, PARAMETER :: id_SSSSM = 1 + INTEGER, PARAMETER :: id_CSSSM = 2 + INTEGER, PARAMETER :: id_ASSSM = 3 + INTEGER, PARAMETER :: id_ASBCM = 4 + INTEGER, PARAMETER :: id_ASPOMM = 5 + INTEGER, PARAMETER :: id_ASSO4M = 6 + INTEGER, PARAMETER :: id_CSSO4M = 7 + INTEGER, PARAMETER :: id_CIDUSTM = 8 + INTEGER, PARAMETER :: id_AIBCM = 9 + INTEGER, PARAMETER :: id_AIPOMM = 10 + INTEGER :: nb_aer + REAL, DIMENSION(klon,klev,8) :: mass_temp + + ! + ! Proprietes optiques + ! + REAL:: alpha_aers_2bands(nbre_RH,2,nbsol_compaer) !--unit m2/g SO4 + REAL:: alpha_aeri_2bands(2,nbinsol_compaer) + REAL:: cg_aers_2bands(nbre_RH,2,nbsol_compaer) !--unit + REAL:: cg_aeri_2bands(2,nbinsol_compaer) + REAL:: piz_aers_2bands(nbre_RH,2,nbsol_compaer) !-- unit + REAL:: piz_aeri_2bands(2,nbinsol_compaer) !-- unit + + + spsol = 0 + spinsol = 0 + + + DATA alpha_aers_2bands/ & + ! seasalt soluble Coarse Soluble (CS) + 0.5090,0.6554,0.7129,0.7767,0.8529,1.2728, & + 1.3820,1.5792,1.9173,2.2002,2.7173,4.1487, & + 0.5167,0.6613,0.7221,0.7868,0.8622,1.3027, & + 1.4227,1.6317,1.9887,2.2883,2.8356,4.3453, & + ! seasalt soluble Accumulation Soluble (AS) + 4.125, 4.674, 5.005, 5.434, 5.985, 10.006, & + 11.175,13.376,17.264,20.540,26.604, 42.349,& + 4.187, 3.939, 3.919, 3.937, 3.995, 5.078, & + 5.511, 6.434, 8.317,10.152,14.024, 26.537, & + ! bc soluble + 7.675,7.675,7.675,7.675,7.675,7.675, & + 7.675,7.675,10.433,11.984,13.767,15.567,& + 4.720,4.720,4.720,4.720,4.720,4.720, & + 4.720,4.720,6.081,6.793,7.567,9.344, & + ! pom soluble + 5.503,5.503,5.503,5.503,5.588,5.957, & + 6.404,7.340,8.545,10.319,13.595,20.398, & + 1.402,1.402,1.402,1.402,1.431,1.562, & + 1.715,2.032,2.425,2.991,4.193,7.133, & + ! sulfate + 4.681,5.062,5.460,5.798,6.224,6.733, & + 7.556,8.613,10.687,12.265,16.32,21.692, & + 1.107,1.239,1.381,1.490,1.635,1.8030, & + 2.071,2.407,3.126,3.940,5.539,7.921/ + + + DATA alpha_aeri_2bands/ & + ! dust insoluble + 0.7661,0.7123,& + ! bc insoluble + 10.360,4.437, & + ! pom insoluble + 3.741,0.606/ + + + DATA cg_aers_2bands/ & + ! seasalt Coarse soluble (CS) + 0.727, 0.747, 0.755, 0.761, 0.770, 0.788, & + 0.792, 0.799, 0.805, 0.809, 0.815, 0.826, & + 0.717, 0.738, 0.745, 0.752, 0.761, 0.779, & + 0.781, 0.786, 0.793, 0.797, 0.803, 0.813, & + ! Sesalt Accumulation Soluble (AS) + 0.727, 0.741, 0.748, 0.754, 0.761, 0.782, & + 0.787, 0.792, 0.797, 0.799, 0.801, 0.799, & + 0.606, 0.645, 0.658, 0.669, 0.681, 0.726, & + 0.734, 0.746, 0.761, 0.770, 0.782, 0.798, & + ! bc soluble + .612, .612, .612, .612, .612, .612, & + .612, .612, .702, .734, .760, .796, & + .433, .433, .433, .433, .433, .433, & + .433, .433, .534, .575, .613, .669, & + ! pom soluble + .663, .663, .663, .663, .666, .674, & + .685, .702, .718, .737, .757, .777, & + .544, .544, .544, .544, .547, .554, & + .565, .583, .604, .631, .661, .698, & + ! sulfate + .658, .669, .680, .688, .698, .707, & + .719, .733, .752, .760, .773, .786, & + .544, .555, .565, .573, .583, .593, & + .610, .628, .655, .666, .692, .719/ + + + DATA cg_aeri_2bands/ & + ! dust insoluble + .701, .670, & + ! bc insoluble + .471, .297, & + ! pom insoluble + .568, .365/ + + DATA piz_aers_2bands/& + ! seasalt Coarse soluble (CS) + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 0.992,0.989,0.987,0.986,0.986,0.980, & + 0.980,0.978,0.976,0.976,0.974,0.971, & + ! seasalt Accumulation Soluble (AS) + 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, & + 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, & + 0.970, 0.975, 0.976, 0.977, 0.978, 0.982, & + 0.982, 0.983, 0.984, 0.984, 0.985, 0.985, & + ! bc soluble + .445, .445, .445, .445, .445, .445, & + .445, .445, .461, .480, .505, .528, & + .362, .362, .362, .362, .362, .362, & + .362, .362, .381, .405, .437, .483, & + ! pom soluble + .972, .972, .972, .972, .972, .974, & + .976, .979, .982, .986, .989, .992, & + .924, .924, .924, .924, .925, .927, & + .932, .938, .945, .952, .961, .970, & + ! sulfate + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + .992, .988, .988, .987, .986, .985, & + .985, .985, .984, .984, .984, .984 / + + + DATA piz_aeri_2bands/ & + ! dust insoluble + .963, .987, & + ! bc insoluble + .395, .264, & + ! pom insoluble + .966, .859/ + + + DO k=1, klev + DO i=1, klon + IF (t_seri(i,k).EQ.0.) THEN + WRITE(lunout,*) 't_seri(i,k)=0 for i=',i,'k=',k + CALL abort_gcm('aeropt_2bands','t_seri=0',1) + END IF + IF (pplay(i,k).EQ.0.) THEN + WRITE(lunout,*) 'pplay(i,k)=0 for i=',i,'k=',k + CALL abort_gcm('aeropt_2bands','pplay=0',1) + END IF + + zrho=pplay(i,k)/t_seri(i,k)/RD ! kg/m3 + mass_temp(i,k,:) = m_allaer(i,k,:) / zrho / 1.e+9 + + ENDDO + ENDDO + + IF (flag_aerosol .EQ. 1) THEN + nb_aer = 1 + ALLOCATE (aerosol_name(nb_aer)) + aerosol_name(1) = id_ASSO4M + ELSEIF (flag_aerosol .EQ. 2) THEN + nb_aer = 2 + ALLOCATE (aerosol_name(nb_aer)) + aerosol_name(1) = id_ASBCM + aerosol_name(2) = id_AIBCM + ELSEIF (flag_aerosol .EQ. 3) THEN + nb_aer = 2 + ALLOCATE (aerosol_name(nb_aer)) + aerosol_name(1) = id_ASPOMM + aerosol_name(2) = id_AIPOMM + ELSEIF (flag_aerosol .EQ. 4) THEN + nb_aer = 5 + ALLOCATE (aerosol_name(nb_aer)) + aerosol_name(1) = id_ASSO4M + aerosol_name(2) = id_ASBCM + aerosol_name(3) = id_AIBCM + aerosol_name(4) = id_ASPOMM + aerosol_name(5) = id_AIPOMM + ELSEIF (flag_aerosol .EQ. 5) THEN + nb_aer = 4 + ALLOCATE (aerosol_name(nb_aer)) + aerosol_name(1) = id_ASBCM + aerosol_name(2) = id_AIBCM + aerosol_name(3) = id_ASPOMM + aerosol_name(4) = id_AIPOMM + ELSEIF (flag_aerosol .EQ. 6) THEN + nb_aer = 3 + ALLOCATE (aerosol_name(nb_aer)) + aerosol_name(1) = id_ASSO4M + aerosol_name(2) = id_ASPOMM + aerosol_name(3) = id_AIPOMM + ENDIF + + + ! + ! loop over modes, use of precalculated nmd and corresponding sigma + ! loop over wavelengths + ! for each mass species in mode + ! interpolate from Sext to retrieve Sext_at_gridpoint_per_species + ! compute optical_thickness_at_gridpoint_per_species + + tau_ae(:,:,:,:)=0. + piz_ae(:,:,:,:)=0. + cg_ae(:,:,:,:)=0. + tau_allaer(:,:,:,:)=0. + piz_allaer(:,:,:,:)=0. + cg_allaer(:,:,:,:)=0. + + ! + ! Calculations that need to be done since we are not in the subroutines INCA + ! + ! air mass auxiliary variable --> zdp1 [kg/(m^2 *s)] + zdp1(:,:)=pdel(:,:)/(gravit*delt) + + + DO m=1,nb_aer ! tau is only computed for each mass + + fac=1.0 + IF (aerosol_name(m).EQ.id_SSSSM) THEN ! for now + soluble=.TRUE. + spsol=1 + ELSEIF (aerosol_name(m).EQ.id_CSSSM) THEN + soluble=.TRUE. + spsol=1 + ELSEIF (aerosol_name(m).EQ.id_ASSSM) THEN + soluble=.TRUE. + spsol=2 + ELSEIF (aerosol_name(m).EQ.id_ASBCM) THEN + soluble=.TRUE. + spsol=3 + ELSEIF (aerosol_name(m).EQ.id_ASPOMM) THEN + soluble=.TRUE. + spsol=4 + ELSEIF ((aerosol_name(m).EQ.id_ASSO4M) .OR. (aerosol_name(m).EQ.id_CSSO4M)) THEN + soluble=.TRUE. + spsol=5 + fac=1.375 ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD + ELSEIF (aerosol_name(m).EQ.id_CIDUSTM) THEN + soluble=.FALSE. + spinsol=1 + ELSEIF (aerosol_name(m).EQ.id_AIBCM) THEN + soluble=.FALSE. + spinsol=2 + ELSEIF (aerosol_name(m).EQ.id_AIPOMM) THEN + soluble=.FALSE. + spinsol=3 + ELSE + CYCLE + ENDIF + + + tau_ae2b_int(:,:,:)=0. + piz_ae2b_int(:,:,:)=0. + cg_ae2b_int(:,:,:)=0. + + DO k=1, KLEV + DO i=1, KLON + + rh=MIN(RHcl(i,k)*100.,RH_MAX) + RH_num = INT( rh/10. + 1.) + + IF (rh.GT.85.) RH_num=10 + IF (rh.GT.90.) RH_num=11 + DELTA=(rh-RH_tab(RH_num))/(RH_tab(RH_num+1)-RH_tab(RH_num)) + + DO inu=1,2 + IF (soluble) THEN + tau_ae2b_int(i,k,inu)= & + alpha_aers_2bands(RH_num,inu,spsol)+ & + DELTA* (alpha_aers_2bands(RH_num+1,inu,spsol) - & + alpha_aers_2bands(RH_num,inu,spsol)) + + piz_ae2b_int(i,k,inu)= & + piz_aers_2bands(RH_num,inu,spsol) + & + DELTA* (piz_aers_2bands(RH_num+1,inu,spsol) - & + piz_aers_2bands(RH_num,inu,spsol)) + + cg_ae2b_int(i,k,inu)= & + cg_aers_2bands(RH_num,inu,spsol) + & + DELTA* (cg_aers_2bands(RH_num+1,inu,spsol) - & + cg_aers_2bands(RH_num,inu,spsol)) + + tau_ae(i,k,aerosol_name(m),inu) = mass_temp(i,k,spsol)*1000.*zdp1(i,k)*delt*tau_ae2b_int(i,k,inu)*fac + + ELSE + tau_ae2b_int(i,k,inu) = alpha_aeri_2bands(inu,spinsol) + piz_ae2b_int(i,k,inu) = piz_aeri_2bands(inu,spinsol) + cg_ae2b_int(i,k,inu) = cg_aeri_2bands(inu,spinsol) + + tau_ae(i,k,aerosol_name(m),inu) = mass_temp(i,k,5+ spinsol)*1000.*zdp1(i,k)*delt*tau_ae2b_int(i,k,inu)*fac + ENDIF + + piz_ae(i,k,aerosol_name(m),inu) = piz_ae2b_int(i,k,inu) + + cg_ae(i,k,aerosol_name(m),inu)= cg_ae2b_int(i,k,inu) + + + ENDDO ! boucle sur les bandes spectrale + ENDDO ! Boucle sur les points géographiques (grille horizontale) + ENDDO ! Boucle sur les niveaux verticaux + ENDDO ! Boucle sur les masses de traceurs + + + DO inu=1, 2 + DO mrfspecies=1,9 + DO k=1, KLEV + DO i=1, KLON + IF (mrfspecies .EQ. 2) THEN ! = total aerosol AER + tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASSO4M,inu)+tau_ae(i,k,id_CSSO4M,inu)+ & + tau_ae(i,k,id_ASBCM,inu)+tau_ae(i,k,id_AIBCM,inu)+ & + tau_ae(i,k,id_ASPOMM,inu)+tau_ae(i,k,id_AIPOMM,inu)+ & + tau_ae(i,k,id_ASSSM,inu)+tau_ae(i,k,id_CSSSM,inu)+tau_ae(i,k,id_SSSSM,inu)+ & + tau_ae(i,k,id_CIDUSTM,inu) + tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20) + + piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSO4M,inu)*piz_ae(i,k,id_ASSO4M,inu)+ & + tau_ae(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu)+ & + tau_ae(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu)+ & + tau_ae(i,k,id_AIBCM,inu)*piz_ae(i,k,id_AIBCM,inu)+ & + tau_ae(i,k,id_ASPOMM,inu)*piz_ae(i,k,id_ASPOMM,inu)+ & + tau_ae(i,k,id_AIPOMM,inu)*piz_ae(i,k,id_AIPOMM,inu)+ & + tau_ae(i,k,id_ASSSM,inu)*piz_ae(i,k,id_ASSSM,inu)+ & + tau_ae(i,k,id_CSSSM,inu)*piz_ae(i,k,id_CSSSM,inu)+ & + tau_ae(i,k,id_SSSSM,inu)*piz_ae(i,k,id_SSSSM,inu)+ & + tau_ae(i,k,id_CIDUSTM,inu)*piz_ae(i,k,id_CIDUSTM,inu))/tau_allaer(i,k,mrfspecies,inu) + piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20) + + cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSO4M,inu)*piz_ae(i,k,id_ASSO4M,inu)*cg_ae(i,k,id_ASSO4M,inu)+ & + tau_ae(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu)*cg_ae(i,k,id_CSSO4M,inu)+ & + tau_ae(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu)*cg_ae(i,k,id_ASBCM,inu)+ & + tau_ae(i,k,id_AIBCM,inu)*piz_ae(i,k,id_AIBCM,inu)*cg_ae(i,k,id_AIBCM,inu)+ & + tau_ae(i,k,id_ASPOMM,inu)*piz_ae(i,k,id_ASPOMM,inu)*cg_ae(i,k,id_ASPOMM,inu)+ & + tau_ae(i,k,id_AIPOMM,inu)*piz_ae(i,k,id_AIPOMM,inu)*cg_ae(i,k,id_AIPOMM,inu)+ & + tau_ae(i,k,id_ASSSM,inu)*piz_ae(i,k,id_ASSSM,inu)*cg_ae(i,k,id_ASSSM,inu)+ & + tau_ae(i,k,id_CSSSM,inu)*piz_ae(i,k,id_CSSSM,inu)*cg_ae(i,k,id_CSSSM,inu)+ & + tau_ae(i,k,id_SSSSM,inu)*piz_ae(i,k,id_SSSSM,inu)*cg_ae(i,k,id_SSSSM,inu)+ & + tau_ae(i,k,id_CIDUSTM,inu)*piz_ae(i,k,id_CIDUSTM,inu)*cg_ae(i,k,id_CIDUSTM,inu))/ & + (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu)) + + + ELSEIF (mrfspecies .EQ. 3) THEN ! = natural aerosol NAT + tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASSO4M,inu)*fractnat_allaer(i,id_ASSO4M)+ & + tau_ae(i,k,id_CSSO4M,inu)*fractnat_allaer(i,id_CSSO4M)+ & + tau_ae(i,k,id_ASBCM,inu)*fractnat_allaer(i,id_ASBCM)+ & + tau_ae(i,k,id_AIBCM,inu)*fractnat_allaer(i,id_AIBCM)+ & + tau_ae(i,k,id_ASPOMM,inu)*fractnat_allaer(i,id_ASPOMM)+ & + tau_ae(i,k,id_AIPOMM,inu)*fractnat_allaer(i,id_AIPOMM)+ & + tau_ae(i,k,id_ASSSM,inu)*fractnat_allaer(i,id_ASSSM)+ & + tau_ae(i,k,id_CSSSM,inu)*fractnat_allaer(i,id_CSSSM)+ & + tau_ae(i,k,id_SSSSM,inu)*fractnat_allaer(i,id_SSSSM)+ & + tau_ae(i,k,id_CIDUSTM,inu)*fractnat_allaer(i,id_CIDUSTM) + tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20) + + piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSO4M,inu)*fractnat_allaer(i,id_ASSO4M)*piz_ae(i,k,id_ASSO4M,inu)+ & + tau_ae(i,k,id_CSSO4M,inu)*fractnat_allaer(i,id_CSSO4M)*piz_ae(i,k,id_CSSO4M,inu)+ & + tau_ae(i,k,id_ASBCM,inu)*fractnat_allaer(i,id_ASBCM)*piz_ae(i,k,id_ASBCM,inu)+ & + tau_ae(i,k,id_AIBCM,inu)*fractnat_allaer(i,id_AIBCM)*piz_ae(i,k,id_AIBCM,inu)+ & + tau_ae(i,k,id_ASPOMM,inu)*fractnat_allaer(i,id_ASPOMM)*piz_ae(i,k,id_ASPOMM,inu)+ & + tau_ae(i,k,id_AIPOMM,inu)*fractnat_allaer(i,id_AIPOMM)*piz_ae(i,k,id_AIPOMM,inu)+ & + tau_ae(i,k,id_ASSSM,inu)*fractnat_allaer(i,id_ASSSM)*piz_ae(i,k,id_ASSSM,inu)+ & + tau_ae(i,k,id_CSSSM,inu)*fractnat_allaer(i,id_CSSSM)*piz_ae(i,k,id_CSSSM,inu)+ & + tau_ae(i,k,id_SSSSM,inu)*fractnat_allaer(i,id_SSSSM)*piz_ae(i,k,id_SSSSM,inu)+ & + tau_ae(i,k,id_CIDUSTM,inu)*fractnat_allaer(i,id_CIDUSTM)*piz_ae(i,k,id_CIDUSTM,inu)) & + /tau_allaer(i,k,mrfspecies,inu) + piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20) + + cg_allaer(i,k,mrfspecies,inu)=(& + tau_ae(i,k,id_ASSO4M,inu)*fractnat_allaer(i,id_ASSO4M)*piz_ae(i,k,id_ASSO4M,inu)*cg_ae(i,k,id_ASSO4M,inu)+ & + tau_ae(i,k,id_CSSO4M,inu)*fractnat_allaer(i,id_CSSO4M)*piz_ae(i,k,id_CSSO4M,inu)*cg_ae(i,k,id_CSSO4M,inu)+ & + tau_ae(i,k,id_ASBCM,inu)*fractnat_allaer(i,id_ASBCM)*piz_ae(i,k,id_ASBCM,inu)*cg_ae(i,k,id_ASBCM,inu)+ & + tau_ae(i,k,id_AIBCM,inu)*fractnat_allaer(i,id_AIBCM)*piz_ae(i,k,id_AIBCM,inu)*cg_ae(i,k,id_AIBCM,inu)+ & + tau_ae(i,k,id_ASPOMM,inu)*fractnat_allaer(i,id_ASPOMM)*piz_ae(i,k,id_ASPOMM,inu)*cg_ae(i,k,id_ASPOMM,inu)+ & + tau_ae(i,k,id_AIPOMM,inu)*fractnat_allaer(i,id_AIPOMM)*piz_ae(i,k,id_AIPOMM,inu)*cg_ae(i,k,id_AIPOMM,inu)+ & + tau_ae(i,k,id_ASSSM,inu)*fractnat_allaer(i,id_ASSSM)*piz_ae(i,k,id_ASSSM,inu)*cg_ae(i,k,id_ASSSM,inu)+ & + tau_ae(i,k,id_CSSSM,inu)*fractnat_allaer(i,id_CSSSM)*piz_ae(i,k,id_CSSSM,inu)*cg_ae(i,k,id_CSSSM,inu)+ & + tau_ae(i,k,id_SSSSM,inu)*fractnat_allaer(i,id_SSSSM)*piz_ae(i,k,id_SSSSM,inu)*cg_ae(i,k,id_SSSSM,inu)+ & + tau_ae(i,k,id_CIDUSTM,inu)*fractnat_allaer(i,id_CIDUSTM)*piz_ae(i,k,id_CIDUSTM,inu)*& + cg_ae(i,k,id_CIDUSTM,inu))/ & + (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu)) + + + ELSEIF (mrfspecies .EQ. 4) THEN ! = BC + tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASBCM,inu)+tau_ae(i,k,id_AIBCM,inu) + tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20) + piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu) & + +tau_ae(i,k,id_AIBCM,inu)*piz_ae(i,k,id_AIBCM,inu))/ & + tau_allaer(i,k,mrfspecies,inu) + piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20) + cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASBCM,inu)*piz_ae(i,k,id_ASBCM,inu) *cg_ae(i,k,id_ASBCM,inu)& + +tau_ae(i,k,id_AIBCM,inu)*piz_ae(i,k,id_AIBCM,inu)*cg_ae(i,k,id_AIBCM,inu))/ & + (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu)) + ELSEIF (mrfspecies .EQ. 5) THEN ! = SO4 + tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASSO4M,inu)+tau_ae(i,k,id_CSSO4M,inu) + tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20) + piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu) & + +tau_ae(i,k,id_ASSO4M,inu)*piz_ae(i,k,id_ASSO4M,inu))/ & + tau_allaer(i,k,mrfspecies,inu) + piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20) + cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_CSSO4M,inu)*piz_ae(i,k,id_CSSO4M,inu) *cg_ae(i,k,id_CSSO4M,inu)& + +tau_ae(i,k,id_ASSO4M,inu)*piz_ae(i,k,id_ASSO4M,inu)*cg_ae(i,k,id_ASSO4M,inu))/ & + (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu)) + + ELSEIF (mrfspecies .EQ. 6) THEN ! = POM + tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASPOMM,inu)+tau_ae(i,k,id_AIPOMM,inu) + tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20) + piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASPOMM,inu)*piz_ae(i,k,id_ASPOMM,inu) & + +tau_ae(i,k,id_AIPOMM,inu)*piz_ae(i,k,id_AIPOMM,inu))/ & + tau_allaer(i,k,mrfspecies,inu) + piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20) + cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASPOMM,inu)*piz_ae(i,k,id_ASPOMM,inu) *cg_ae(i,k,id_ASPOMM,inu)& + +tau_ae(i,k,id_AIPOMM,inu)*piz_ae(i,k,id_AIPOMM,inu)*cg_ae(i,k,id_AIPOMM,inu))/ & + (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu)) + ELSEIF (mrfspecies .EQ. 7) THEN ! = DUST + tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_CIDUSTM,inu) + tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20) + piz_allaer(i,k,mrfspecies,inu)=piz_ae(i,k,id_CIDUSTM,inu) + cg_allaer(i,k,mrfspecies,inu)=cg_ae(i,k,id_CIDUSTM,inu) + + ELSEIF (mrfspecies .EQ. 8) THEN ! = SS + tau_allaer(i,k,mrfspecies,inu)=tau_ae(i,k,id_ASSSM,inu)+tau_ae(i,k,id_CSSSM,inu)+tau_ae(i,k,id_SSSSM,inu) + tau_allaer(i,k,mrfspecies,inu)=MAX(tau_allaer(i,k,mrfspecies,inu),1e-20) + piz_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSSM,inu)*piz_ae(i,k,id_ASSSM,inu) & + +tau_ae(i,k,id_CSSSM,inu)*piz_ae(i,k,id_CSSSM,inu) & + +tau_ae(i,k,id_SSSSM,inu)*piz_ae(i,k,id_SSSSM,inu))/ & + tau_allaer(i,k,mrfspecies,inu) + piz_allaer(i,k,mrfspecies,inu)=MAX(piz_allaer(i,k,mrfspecies,inu),1e-20) + cg_allaer(i,k,mrfspecies,inu)=(tau_ae(i,k,id_ASSSM,inu)*piz_ae(i,k,id_ASSSM,inu) *cg_ae(i,k,id_ASSSM,inu)& + +tau_ae(i,k,id_CSSSM,inu)*piz_ae(i,k,id_CSSSM,inu)*cg_ae(i,k,id_CSSSM,inu) & + +tau_ae(i,k,id_SSSSM,inu)*piz_ae(i,k,id_SSSSM,inu)*cg_ae(i,k,id_SSSSM,inu))/ & + (tau_allaer(i,k,mrfspecies,inu)*piz_allaer(i,k,mrfspecies,inu)) + + ELSEIF (mrfspecies .EQ. 9) THEN ! = NO3 + tau_allaer(i,k,mrfspecies,inu)=0. ! preliminary + piz_allaer(i,k,mrfspecies,inu)=0. + cg_allaer(i,k,mrfspecies,inu)=0. + ENDIF + ENDDO + ENDDO + ENDDO + ENDDO + + DEALLOCATE(aerosol_name) + +END SUBROUTINE AEROPT_2BANDS Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_5wv.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_5wv.F90 (revision 1150) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/aeropt_5wv.F90 (revision 1150) @@ -0,0 +1,515 @@ + + +SUBROUTINE AEROPT_5WV(& + pdel, m_allaer, delt, & + RHcl, ai, flag_aerosol, & + pplay, t_seri) + + USE DIMPHY + + ! + ! Yves Balkanski le 12 avril 2006 + ! Celine Deandreis + ! Anne Cozic Avril 2009 + ! a partir d'une sous-routine de Johannes Quaas pour les sulfates + ! + ! + ! Refractive indices for seasalt come from Shettle and Fenn (1979) + ! + ! Refractive indices from water come from Hale and Querry (1973) + ! + ! Refractive indices from Ammonium Sulfate Toon and Pollack (1976) + ! + ! Refractive indices for Dust, internal mixture of minerals coated with 1.5% hematite + ! by Volume (Balkanski et al., 2006) + ! + ! Refractive indices for POM: Kinne (pers. Communication + ! + ! Refractive index for BC from Shettle and Fenn (1979) + ! + ! Shettle, E. P., & Fenn, R. W. (1979), Models for the aerosols of the lower atmosphere and + ! the effects of humidity variations on their optical properties, U.S. Air Force Geophysics + ! Laboratory Rept. AFGL-TR-79-0214, Hanscomb Air Force Base, MA. + ! + ! Hale, G. M. and M. R. Querry, Optical constants of water in the 200-nm to 200-m + ! wavelength region, Appl. Opt., 12, 555-563, 1973. + ! + ! Toon, O. B. and J. B. Pollack, The optical constants of several atmospheric aerosol species: + ! Ammonium sulfate, aluminum oxide, and sodium chloride, J. Geohys. Res., 81, 5733-5748, + ! 1976. + ! + ! Balkanski, Y., M. Schulz, T. Claquin And O. Boucher, Reevaluation of mineral aerosol + ! radiative forcings suggests a better agreement with satellite and AERONET data, Atmospheric + ! Chemistry and Physics Discussions., 6, pp 8383-8419, 2006. + ! + IMPLICIT NONE + INCLUDE "YOMCST.h" + ! + ! Input arguments: + ! + REAL, DIMENSION(klon,klev), INTENT(in) :: pdel + REAL, INTENT(in) :: delt + REAL, DIMENSION(klon,klev,8), INTENT(in) :: m_allaer + REAL, DIMENSION(klon,klev), INTENT(in) :: RHcl ! humidite relative ciel clair + INTEGER,INTENT(in) :: flag_aerosol + REAL, DIMENSION(klon,klev), INTENT(in) :: pplay + REAL, DIMENSION(klon,klev), INTENT(in) :: t_seri + + ! + ! Output arguments: + ! + REAL, DIMENSION(klon), INTENT(out) :: ai ! POLDER aerosol index + + ! + ! Local + ! + INTEGER, PARAMETER :: las = 5 + LOGICAL :: soluble + + INTEGER :: i, k, m + INTEGER :: spsol, spinsol, la + INTEGER :: RH_num + INTEGER, PARAMETER :: la443 = 1 + INTEGER, PARAMETER :: la550 = 2 + INTEGER, PARAMETER :: la670 = 3 + INTEGER, PARAMETER :: la765 = 4 + INTEGER, PARAMETER :: la865 = 5 + INTEGER, PARAMETER :: nbre_RH=12 + INTEGER, PARAMETER :: nbsol_compaer=5 ! 1- Seasalt AS: 2- Sesalt CS; 3- BC soluble; 4- POM soluble; 5- SO4. + INTEGER, PARAMETER :: nbinsol_compaer=3 ! 1- Dust; 2- BC insoluble; 3- POM insoluble + REAL :: zrho + REAL :: RH_tab(nbre_RH) + REAL :: DELTA, rh + REAL :: tau_ae5wv_int(KLON,KLEV,las) ! Intermediate computation of epaisseur optique aerosol + REAL :: piz_ae5wv_int(KLON,KLEV,las) ! Intermediate single scattering albedo aerosol + REAL :: cg_ae5wv_int(KLON,KLEV,las) ! Intermediate asymmetry parameter aerosol + REAL, PARAMETER :: RH_MAX=95. + DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./ + REAL :: taue670(KLON) ! epaisseur optique aerosol absorption 550 nm + REAL :: taue865(KLON) ! epaisseur optique aerosol extinction 865 nm + REAL :: fac + REAL :: zdp1(klon,klev) + REAL, PARAMETER :: gravit = 9.80616 ! m2/s + INTEGER, ALLOCATABLE, DIMENSION(:) :: aerosol_name + INTEGER, PARAMETER :: id_SSSSM = 1 + INTEGER, PARAMETER :: id_CSSSM = 2 + INTEGER, PARAMETER :: id_ASSSM = 3 + INTEGER, PARAMETER :: id_ASBCM = 4 + INTEGER, PARAMETER :: id_ASPOMM = 5 + INTEGER, PARAMETER :: id_ASSO4M = 6 + INTEGER, PARAMETER :: id_CSSO4M = 7 + INTEGER, PARAMETER :: id_CIDUSTM = 8 + INTEGER, PARAMETER :: id_AIBCM = 9 + INTEGER, PARAMETER :: id_AIPOMM = 10 + INTEGER :: nb_aer + + REAL :: tau3d(KLON,KLEV), piz3d(KLON,KLEV), cg3d(KLON,KLEV) + REAL :: abs3d(KLON,KLEV) ! epaisseur optique d'absorption + + + REAL :: alpha_aers_5wv(nbre_RH,las,nbsol_compaer) ! ext. coeff. Soluble comp. units *** m2/g + ! 1- Seasalt AS: 2- Sesalt CS; 3- BC; 4- POM; 5- SO4. + REAL :: alpha_aeri_5wv(las,nbinsol_compaer) ! ext. coeff. Insoluble comp. 1- Dust: 2- BC; 3- POM + REAL :: cg_aers_5wv(nbre_RH,las,nbsol_compaer) ! Asym. param. soluble comp. + ! 1- Seasalt AS: 2- Sesalt CS; 3- BC; 4- POM; 5- SO4. + REAL :: cg_aeri_5wv(las,nbinsol_compaer) ! Asym. param. insoluble comp. 1- Dust: 2- BC; 3- POM + REAL :: piz_aers_5wv(nbre_RH,las,nbsol_compaer) + ! 1- Seasalt AS: 2- Sesalt CS; 3- BC; 4- POM; 5- SO4. + REAL :: piz_aeri_5wv(las,nbinsol_compaer) ! Insoluble comp. 1- Dust: 2- BC; 3- POM + + REAL, ALLOCATABLE :: TAUSUM(:,:,:) ! Aerosol optical thickness per species + REAL, ALLOCATABLE :: TAU(:,:,:,:) ! Aerosol optical thickness per species + REAL, DIMENSION(klon,klev,8) :: mass_temp + + ! + ! Proprietes optiques + ! + REAL :: radry = 287.054 ! dry air mass constant + + ! + ! + ! + ! From here on we look at the optical parameters at 5 wavelengths: + ! 443nm, 550, 670, 765 and 865 nm + ! le 12 AVRIL 2006 + ! + DATA alpha_aers_5wv/ & + ! seasalt soluble CS + 0.50,0.90,1.05,1.21,1.40,2.41, & + 2.66,3.11,3.88,4.52,5.69,8.84, & + 0.51,0.92,1.07,1.23,1.42,2.45, & + 2.70,3.16,3.94,4.58,5.76,8.94, & + 0.52,0.93,1.08,1.24,1.43,2.47, & + 2.73,3.20,3.99,4.64,5.84,9.04, & + 0.52,0.93,1.09,1.25,1.44,2.50, & + 2.76,3.23,4.03,4.68,5.89,9.14, & + 0.52,0.94,1.09,1.26,1.45,2.51, & + 2.78,3.25,4.06,4.72,5.94,9.22, & + ! seasalt soluble AS + 4.28, 7.17, 8.44, 9.85,11.60,22.44, & + 25.34,30.54,39.38,46.52,59.33,91.77, & + 3.40, 5.67, 6.69, 7.85, 9.32,19.03, & + 21.78,26.88,35.87,43.40,57.33,93.43, & + 2.48, 4.22, 5.02, 5.94, 7.11,15.29, & + 17.70,22.31,30.73,38.06,52.15,90.59, & + 1.90, 3.29, 3.94, 4.69, 5.65, 12.58, & + 14.68,18.77,26.41,33.25,46.77,85.50, & + 1.47, 2.59, 3.12, 3.74, 4.54, 10.42, & + 12.24,15.82,22.66,28.91,41.54,79.33, & + ! bc soluble + 7.930,7.930,7.930,7.930,7.930,7.930, & + 7.930,7.930,10.893,12.618,14.550,16.613, & + 7.658,7.658,7.658,7.658,7.658,7.658, & + 7.658,7.658,10.351,11.879,13.642,15.510, & + 7.195,7.195,7.195,7.195,7.195,7.195, & + 7.195,7.195,9.551,10.847,12.381,13.994, & + 6.736,6.736,6.736,6.736,6.736,6.736, & + 6.736,6.736,8.818,9.938,11.283,12.687, & + 6.277,6.277,6.277,6.277,6.277,6.277, & + 6.277,6.277,8.123,9.094,10.275,11.501, & + ! pom soluble + 6.676,6.676,6.676,6.676,6.710,6.934, & + 7.141,7.569,8.034,8.529,9.456,10.511, & + 5.109,5.109,5.109,5.109,5.189,5.535, & + 5.960,6.852,8.008,9.712,12.897,19.676, & + 3.718,3.718,3.718,3.718,3.779,4.042, & + 4.364,5.052,5.956,7.314,9.896,15.688, & + 2.849,2.849,2.849,2.849,2.897,3.107, & + 3.365,3.916,4.649,5.760,7.900,12.863, & + 2.229,2.229,2.229,2.229,2.268,2.437, & + 2.645,3.095,3.692,4.608,6.391,10.633, & + ! Sulfate + 5.751,6.215,6.690,7.024,7.599,8.195, & + 9.156,10.355,12.660,14.823,18.908,24.508, & + 4.320,4.675,5.052,5.375,5.787,6.274, & + 7.066,8.083,10.088,12.003,15.697,21.133, & + 3.079,3.351,3.639,3.886,4.205,4.584, & + 5.206,6.019,7.648,9.234,12.391,17.220, & + 2.336,2.552,2.781,2.979,3.236,3.540, & + 4.046,4.711,6.056,7.388,10.093,14.313, & + 1.777,1.949,2.134,2.292,2.503,2.751, & + 3.166,3.712,4.828,5.949,8.264,11.922/ + + DATA alpha_aeri_5wv/ & + ! dust insoluble + 0.759, 0.770, 0.775, 0.775, 0.772, & + !!jb bc insoluble + 11.536,10.033, 8.422, 7.234, 6.270, & + ! pom insoluble + 5.042, 3.101, 1.890, 1.294, 0.934/ + + DATA cg_aers_5wv/ & + ! seasalt soluble (CS) + 0.730,0.753,0.760,0.766,0.772,0.793, & + 0.797,0.802,0.809,0.813,0.820,0.830, & + 0.719,0.744,0.751,0.757,0.764,0.786, & + 0.791,0.796,0.803,0.808,0.815,0.826, & + 0.721,0.744,0.750,0.756,0.762,0.784, & + 0.787,0.793,0.800,0.804,0.811,0.822, & + 0.717,0.741,0.747,0.753,0.759,0.780, & + 0.784,0.789,0.795,0.800,0.806,0.817, & + 0.715,0.739,0.745,0.751,0.757,0.777, & + 0.781,0.786,0.793,0.797,0.803,0.814, & + ! seasalt soluble (AS) + 0.698,0.722,0.729,0.736,0.743,0.765, & + 0.768,0.773,0.777,0.779,0.781,0.779, & + 0.682,0.710,0.719,0.727,0.735,0.764, & + 0.769,0.776,0.783,0.787,0.791,0.792, & + 0.658,0.691,0.701,0.710,0.720,0.756, & + 0.763,0.771,0.782,0.788,0.795,0.801, & + 0.632,0.668,0.679,0.690,0.701,0.743, & + 0.750,0.762,0.775,0.783,0.792,0.804, & + 0.605,0.644,0.656,0.669,0.681,0.729, & + 0.737,0.750,0.765,0.775,0.787,0.803, & + ! bc soluble + .651, .651, .651, .651, .651, .651, & + .651, .651, .738, .764, .785, .800, & + .597, .597, .597, .597, .597, .597, & + .597, .597, .695, .725, .751, .770, & + .543, .543, .543, .543, .543, .543, & + .543, .543, .650, .684, .714, .736, & + .504, .504, .504, .504, .504, .504, & + .504, .504, .614, .651, .683, .708, & + .469, .469, .469, .469, .469, .469, & + .469, .469, .582, .620, .655, .681, & + ! pom soluble + .679, .679, .679, .679, .683, .691, & + .703, .720, .736, .751, .766, .784, & + .656, .656, .656, .656, .659, .669, & + .681, .699, .717, .735, .750, .779, & + .623, .623, .623, .623, .627, .637, & + .649, .668, .688, .709, .734, .762, & + .592, .592, .592, .592, .595, .605, & + .618, .639, .660, .682, .711, .743, & + .561, .561, .561, .561, .565, .575, & + .588, .609, .632, .656, .688, .724, & + ! sulfate + .671, .684, .697, .704, .714, .723, & + .734, .746, .762, .771, .781, .789, & + .653, .666, .678, .687, .697, .707, & + .719, .732, .751, .762, .775, .789, & + .622, .635, .648, .657, .667, .678, & + .691, .705, .728, .741, .758, .777, & + .591, .604, .617, .627, .638, .650, & + .664, .679, .704, .719, .739, .761, & + .560, .574, .587, .597, .609, .621, & + .637, .653, .680, .697, .719, .745/ + ! + + DATA cg_aeri_5wv/& + ! dust insoluble + 0.714, 0.697, 0.688, 0.683, 0.679, & + ! bc insoluble + 0.511, 0.445, 0.384, 0.342, 0.307, & + !c pom insoluble + 0.596, 0.536, 0.466, 0.409, 0.359/ + ! + DATA piz_aers_5wv/& + ! seasalt soluble (CS) + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + ! seasalt soluble (AS) + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + ! bc soluble + .445, .445, .445, .445, .445, .445, & + .445, .445, .470, .487, .508, .531, & + .442, .442, .442, .442, .442, .442, & + .442, .442, .462, .481, .506, .533, & + .427, .427, .427, .427, .427, .427, & + .427, .427, .449, .470, .497, .526, & + .413, .413, .413, .413, .413, .413, & + .413, .413, .437, .458, .486, .516, & + .399, .399, .399, .399, .399, .399, & + .399, .399, .423, .445, .473, .506, & + ! pom soluble + .975, .975, .975, .975, .975, .977, & + .979, .982, .984, .987, .990, .994, & + .972, .972, .972, .972, .973, .974, & + .977, .980, .983, .986, .989, .993, & + .963, .963, .963, .963, .964, .966, & + .969, .974, .977, .982, .986, .991, & + .955, .955, .955, .955, .955, .958, & + .962, .967, .972, .977, .983, .989, & + .944, .944, .944, .944, .944, .948, & + .952, .959, .962, .972, .979, .987, & + ! sulfate + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000, & + 1.000,1.000,1.000,1.000,1.000,1.000/ + ! + DATA piz_aeri_5wv/& + ! dust insoluble + 0.944, 0.970, 0.977, 0.982, 0.987, & + ! bc insoluble + 0.415, 0.387, 0.355, 0.328, 0.301, & + ! pom insoluble + 0.972, 0.963, 0.943, 0.923, 0.897/ + + + + DO k=1, klev + DO i=1, klon + IF (t_seri(i,k).EQ.0) stop 'stop aeropt_5wv T ' + IF (pplay(i,k).EQ.0) stop 'stop aeropt_5wv p ' + + zrho=pplay(i,k)/t_seri(i,k)/RD ! kg/m3 + mass_temp(i,k,:) = m_allaer(i,k,:) / zrho / 1.e+9 + + ENDDO + ENDDO + + + + IF (flag_aerosol .EQ. 1) THEN + nb_aer = 1 + ALLOCATE (aerosol_name(nb_aer)) + aerosol_name(1) = id_ASSO4M + ELSEIF (flag_aerosol .EQ. 2) THEN + nb_aer = 2 + ALLOCATE (aerosol_name(nb_aer)) + aerosol_name(1) = id_ASBCM + aerosol_name(2) = id_AIBCM + ELSEIF (flag_aerosol .EQ. 3) THEN + nb_aer = 2 + ALLOCATE (aerosol_name(nb_aer)) + aerosol_name(1) = id_ASPOMM + aerosol_name(2) = id_AIPOMM + ELSEIF (flag_aerosol .EQ. 4) THEN + nb_aer = 5 + ALLOCATE (aerosol_name(nb_aer)) + aerosol_name(1) = id_ASSO4M + aerosol_name(2) = id_ASBCM + aerosol_name(3) = id_AIBCM + aerosol_name(4) = id_ASPOMM + aerosol_name(5) = id_AIPOMM + ELSEIF (flag_aerosol .EQ. 5) THEN + nb_aer = 4 + ALLOCATE (aerosol_name(nb_aer)) + aerosol_name(1) = id_ASBCM + aerosol_name(2) = id_AIBCM + aerosol_name(3) = id_ASPOMM + aerosol_name(4) = id_AIPOMM + ELSEIF (flag_aerosol .EQ. 6) THEN + nb_aer = 3 + ALLOCATE (aerosol_name(nb_aer)) + aerosol_name(1) = id_ASSO4M + aerosol_name(2) = id_ASPOMM + aerosol_name(3) = id_AIPOMM + ENDIF + + ALLOCATE (tausum(klon,las,nb_aer)) + ALLOCATE (tau(klon,klev,las,nb_aer)) + + + + + ! + ! loop over modes, use of precalculated nmd and corresponding sigma + ! loop over wavelengths + ! for each mass species in mode + ! interpolate from Sext to retrieve Sext_at_gridpoint_per_species + ! compute optical_thickness_at_gridpoint_per_species + + ai(:)=0. + + tau_ae5wv_int(:,:,:)=0. + piz_ae5wv_int(:,:,:)=0. + cg_ae5wv_int(:,:,:)=0. + tausum(:,:,:)=0. + + tau(:,:,:,:)=0. + ! + ! Calculations that need to be done since we are not in the subroutines INCA + ! + ! air mass auxiliary variable --> zdp1 [kg/(m^2 *s)] + zdp1=pdel/(gravit*delt) + + DO m=1,nb_aer ! tau is only computed for each mass + + fac=1.0 + IF (aerosol_name(m).EQ.id_SSSSM) THEN ! for now + soluble=.TRUE. + spsol=1 + ELSEIF (aerosol_name(m).EQ.id_CSSSM) THEN + soluble=.TRUE. + spsol=1 + ELSEIF (aerosol_name(m).EQ.id_ASSSM) THEN + soluble=.TRUE. + spsol=2 + ELSEIF (aerosol_name(m).EQ.id_ASBCM) THEN + soluble=.TRUE. + spsol=3 + ELSEIF (aerosol_name(m).EQ.id_ASPOMM) THEN + soluble=.TRUE. + spsol=4 + ELSEIF ((aerosol_name(m).EQ.id_ASSO4M) .OR. (aerosol_name(m).EQ.id_CSSO4M)) THEN + soluble=.TRUE. + spsol=5 + fac=1.375 ! (NH4)2-SO4/SO4 132/96 mass conversion factor for OD + ELSEIF (aerosol_name(m).EQ.id_CIDUSTM) THEN + soluble=.FALSE. + spinsol=1 + ELSEIF (aerosol_name(m).EQ.id_AIBCM) THEN + soluble=.FALSE. + spinsol=2 + ELSEIF (aerosol_name(m).EQ.id_AIPOMM) THEN + soluble=.FALSE. + spinsol=3 + ELSE + CYCLE + ENDIF + + DO la=1,las + tau3d(:,:)=0. + piz3d(:,:)=0. + cg3d(:,:)=0. + abs3d(:,:)=0. + + DO k=1, KLEV + DO i=1, KLON + + rh=MIN(RHcl(i,k)*100.,RH_MAX) + RH_num = INT( rh/10. + 1.) + + IF (rh.GT.85.) RH_num=10 + IF (rh.GT.90.) RH_num=11 + DELTA=(rh-RH_tab(RH_num))/(RH_tab(RH_num+1)-RH_tab(RH_num)) + + IF (soluble) THEN + tau_ae5wv_int(i,k,la) = & + alpha_aers_5wv(RH_num,la,spsol)+DELTA* & + (alpha_aers_5wv(RH_num+1,la,spsol) - & + alpha_aers_5wv(RH_num,la,spsol)) + + piz_ae5wv_int(i,k,la) = & + piz_aers_5wv(RH_num,la,spsol)+DELTA* & + (piz_aers_5wv(RH_num+1,la,spsol) - & + piz_aers_5wv(RH_num,la,spsol)) + + cg_ae5wv_int(i,k,la) = & + cg_aers_5wv(RH_num,la,spsol)+DELTA* & + (cg_aers_5wv(RH_num+1,la,spsol) - & + cg_aers_5wv(RH_num,la,spsol)) + + tau3d(i,k) = & + mass_temp(i,k,spsol)*1000.*zdp1(i,k)*tau_ae5wv_int(i,k,la)*delt*fac + + ELSE + tau_ae5wv_int(i,k,la) = alpha_aeri_5wv(la,spinsol) + piz_ae5wv_int(i,k,la) = piz_aeri_5wv(la,spinsol) + cg_ae5wv_int(i,k,la) = cg_aeri_5wv(la,spinsol) + + tau3d(i,k) = & + mass_temp(i,k,5+spinsol)*1000.*zdp1(i,k)*tau_ae5wv_int(i,k,la)*delt*fac + ENDIF + + + ENDDO ! Boucle sur les points géographiques (grille horizontale) + ENDDO ! Boucle sur les niveaux verticaux + + tau(:,:,la,m)=tau3d(:,:) + + DO k=1, KLEV + DO i=1,KLON + tausum(i,la,m)=tausum(i,la,m)+tau3d(i,k) + ENDDO + ENDDO + + ENDDO ! boucle sur les longueurs d'onde + ENDDO ! Boucle sur les masses de traceurs + + + taue670(:) = SUM(tausum(:,la670,:),dim=2) + taue865(:) = SUM(tausum(:,la865,:),dim=2) + + DO i=1, klon + ai(i)=-LOG(MAX(taue670(i),0.0001)/ & + MAX(taue865(i),0.0001))/LOG(670./865.) + ENDDO + + +END SUBROUTINE AEROPT_5WV Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/aerosol_optic.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/aerosol_optic.F90 (revision 1150) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/aerosol_optic.F90 (revision 1150) @@ -0,0 +1,145 @@ +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/clesphys.h =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/clesphys.h (revision 1149) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/clesphys.h (revision 1150) @@ -40,5 +40,5 @@ INTEGER lev_histhf, lev_histday, lev_histmth CHARACTER*4 type_run -! aer_type: pour utiliser un fichier constant dans readsulfate +! aer_type: pour utiliser un fichier constant dans readaerosol CHARACTER*8 :: aer_type LOGICAL ok_isccp, ok_regdyn Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/conf_phys.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/conf_phys.F90 (revision 1149) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/conf_phys.F90 (revision 1150) @@ -12,5 +12,6 @@ & iflag_cldcon, & & iflag_ratqs,ratqsbas,ratqshaut, & - & ok_ade, ok_aie, aerosol_couple, & + & ok_ade, ok_aie, aerosol_couple, & + & flag_aerosol, new_aod, & & bl95_b0, bl95_b1,& & iflag_thermals,nsplit_thermals,tau_thermals, & @@ -60,4 +61,6 @@ logical :: ok_LES LOGICAL :: ok_ade, ok_aie, aerosol_couple + INTEGER :: flag_aerosol + LOGICAL :: new_aod REAL :: bl95_b0, bl95_b1 real :: fact_cldcon, facttemps,ratqsbas,ratqshaut @@ -71,4 +74,6 @@ logical,SAVE :: ok_LES_omp LOGICAL,SAVE :: ok_ade_omp, ok_aie_omp, aerosol_couple_omp + INTEGER, SAVE :: flag_aerosol_omp + LOGICAL, SAVE :: new_aod_omp REAL,SAVE :: bl95_b0_omp, bl95_b1_omp REAL,SAVE :: freq_ISCCP_omp, ecrit_ISCCP_omp @@ -239,9 +244,34 @@ CALL getin('aerosol_couple',aerosol_couple_omp) +! +!Config Key = flag_aerosol +!Config Desc = which aerosol is use for coupled model +!Config Def = 1 +!Config Help = Used in physiq.F +! +! - flag_aerosol=1 => so4 only (defaut) +! - flag_aerosol=2 => bc only +! - flag_aerosol=3 => pom only +! - flag_aerosol=4 => all aerosol +! - flag_aerosol=5 => bcpom +! - flag_aerosol=6 => pomsulf + + flag_aerosol_omp = 1 + CALL getin('flag_aerosol',flag_aerosol_omp) + +! Temporary variable for testing purpose!! +!Config Key = new_aod +!Config Desc = which calcul of aeropt +!Config Def = false +!Config Help = Used in physiq.F +! + new_aod_omp = .true. + CALL getin('new_aod',new_aod_omp) + ! !Config Key = aer_type !Config Desc = Use a constant field for the aerosols !Config Def = scenario -!Config Help = Used in readsulfate.F +!Config Help = Used in readaerosol.F90 ! aer_type_omp = 'scenario' @@ -918,5 +948,5 @@ !Config Help = ! - iflag_coupl = 0 + iflag_coupl_omp = 0 call getin('iflag_coupl',iflag_coupl_omp) @@ -927,5 +957,5 @@ !Config Help = ! - iflag_clos = 1 + iflag_clos_omp = 1 call getin('iflag_clos',iflag_clos_omp) ! @@ -935,5 +965,5 @@ !Config Help = ! - iflag_cvl_sigd = 0 + iflag_cvl_sigd_omp = 0 call getin('iflag_cvl_sigd',iflag_cvl_sigd_omp) @@ -943,5 +973,5 @@ !Config Help = ! - iflag_wake = 0 + iflag_wake_omp = 0 call getin('iflag_wake',iflag_wake_omp) @@ -1266,4 +1296,6 @@ ok_aie = ok_aie_omp aerosol_couple = aerosol_couple_omp + flag_aerosol=flag_aerosol_omp + new_aod=new_aod_omp aer_type = aer_type_omp bl95_b0 = bl95_b0_omp @@ -1326,4 +1358,12 @@ WRITE(numout,*)' ERROR version_ocean=',version_ocean,' not valid with slab ocean' CALL abort_gcm('conf_phys','version_ocean not valid',1) + END IF + +! Test sur new_aod. Ce flag permet de retrouver les resultats de l'AR4 +! il n'est utilisable que lors du couplage avec le SO4 seul + IF (ok_ade .OR. ok_aie) THEN + IF ( .NOT. new_aod .AND. flag_aerosol .NE. 1) THEN + CALL abort_gcm('conf_phys','new_aod=.FALSE. not compatible avec flag_aerosol=1',1) + END IF END IF @@ -1396,4 +1436,6 @@ write(numout,*)' ok_aie = ',ok_aie write(numout,*)' aerosol_couple = ', aerosol_couple + write(numout,*)' flag_aerosol = ', flag_aerosol + write(numout,*)' new_aod = ', new_aod write(numout,*)' aer_type = ',aer_type write(numout,*)' bl95_b0 = ',bl95_b0 @@ -1409,5 +1451,5 @@ write(numout,*)' type_run = ',type_run write(numout,*)' ok_isccp = ',ok_isccp - WRITE(numout,*)' solarlong0 = ', solarlong0 + write(numout,*)' solarlong0 = ', solarlong0 write(numout,*)' qsol0 = ', qsol0 write(numout,*)' inertie_sol = ', inertie_sol Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_output_write.h =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_output_write.h (revision 1149) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_output_write.h (revision 1150) @@ -832,19 +832,23 @@ IF (ok_ade) THEN - IF (o_topswad%flag(iff)<=lev_files(iff)) THEN - CALL histwrite_phy(nid_files(iff),o_topswad%name,itau_w,topswad) - ENDIF - IF (o_solswad%flag(iff)<=lev_files(iff)) THEN - CALL histwrite_phy(nid_files(iff),o_solswad%name,itau_w,solswad) - ENDIF + IF (o_topswad%flag(iff)<=lev_files(iff)) THEN + CALL histwrite_phy(nid_files(iff),o_topswad%name,itau_w, + $ topswad_aero) + ENDIF + IF (o_solswad%flag(iff)<=lev_files(iff)) THEN + CALL histwrite_phy(nid_files(iff),o_solswad%name,itau_w, + $ solswad_aero) + ENDIF ENDIF IF (ok_aie) THEN - IF (o_topswai%flag(iff)<=lev_files(iff)) THEN - CALL histwrite_phy(nid_files(iff),o_topswai%name,itau_w,topswai) - ENDIF - IF (o_solswai%flag(iff)<=lev_files(iff)) THEN - CALL histwrite_phy(nid_files(iff),o_solswai%name,itau_w,solswai) - ENDIF + IF (o_topswai%flag(iff)<=lev_files(iff)) THEN + CALL histwrite_phy(nid_files(iff),o_topswai%name,itau_w, + $ topswai_aero) + ENDIF + IF (o_solswai%flag(iff)<=lev_files(iff)) THEN + CALL histwrite_phy(nid_files(iff),o_solswai%name,itau_w, + $ solswai_aero) + ENDIF ENDIF Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_state_var_mod.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_state_var_mod.F90 (revision 1149) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/phys_state_var_mod.F90 (revision 1150) @@ -255,8 +255,4 @@ !$OMP THREADPRIVATE(snow_con) ! -! sulfate_pi : SO4 aerosol concentration [ug/m3] (pre-industrial value) - - REAL,SAVE,ALLOCATABLE :: sulfate_pi(:, :) -!$OMP THREADPRIVATE(sulfate_pi) REAL,SAVE,ALLOCATABLE :: rlonPOS(:) !$OMP THREADPRIVATE(rlonPOS) @@ -269,5 +265,5 @@ ! ok_aie=T -> ! ok_ade=T -AIE=topswai-topswad -! ok_ade=F -AIE=topswai-topsw +! ok_ade=F -AIE=topswai-topsw ! !topswad, solswad : Aerosol direct effect @@ -277,29 +273,5 @@ REAL,SAVE,ALLOCATABLE :: topswai(:), solswai(:) !$OMP THREADPRIVATE(topswai,solswai) - REAL,SAVE,ALLOCATABLE :: tau_ae(:,:,:), piz_ae(:,:,:) -!$OMP THREADPRIVATE(tau_ae,piz_ae) - REAL,SAVE,ALLOCATABLE :: cg_ae(:,:,:) -!$OMP THREADPRIVATE(cg_ae) - -! Les variables suivants uniquement pour un configuration avec INCA -! topswad_inca, solswad_inca : Aerosol direct effect - REAL,SAVE,ALLOCATABLE :: topswad_inca(:), solswad_inca(:) -!$OMP THREADPRIVATE(topswad_inca,solswad_inca) -! topswad0_inca, solswad0_inca : Aerosol direct effect - REAL,SAVE,ALLOCATABLE :: topswad0_inca(:), solswad0_inca(:) -!$OMP THREADPRIVATE(topswad0_inca,solswad0_inca) -! topswai_inca, solswai_inca : Aerosol indirect effect - REAL,SAVE,ALLOCATABLE :: topswai_inca(:), solswai_inca(:) -!$OMP THREADPRIVATE(topswai_inca,solswai_inca) - REAL,SAVE,ALLOCATABLE :: topsw_inca(:,:), solsw_inca(:,:) -!$OMP THREADPRIVATE(topsw_inca,solsw_inca) - REAL,SAVE,ALLOCATABLE :: topsw0_inca(:,:), solsw0_inca(:,:) -!$OMP THREADPRIVATE(topsw0_inca,solsw0_inca) - REAL,SAVE,ALLOCATABLE :: tau_inca(:,:,:,:) -!$OMP THREADPRIVATE(tau_inca) - REAL,SAVE,ALLOCATABLE :: piz_inca(:,:,:,:) -!$OMP THREADPRIVATE(piz_inca) - REAL,SAVE,ALLOCATABLE :: cg_inca(:,:,:,:) -!$OMP THREADPRIVATE(cg_inca) + REAL,SAVE,ALLOCATABLE :: ccm(:,:,:) !$OMP THREADPRIVATE(ccm) @@ -402,6 +374,4 @@ ALLOCATE(ibas_con(klon), itop_con(klon)) ALLOCATE(rain_con(klon), snow_con(klon)) -! - ALLOCATE(sulfate_pi(klon, klev)) ALLOCATE(rlonPOS(klon)) ALLOCATE(newsst(klon)) @@ -409,19 +379,5 @@ ALLOCATE(topswad(klon), solswad(klon)) ALLOCATE(topswai(klon), solswai(klon)) - ALLOCATE(tau_ae(klon,klev,2), piz_ae(klon,klev,2)) - ALLOCATE(cg_ae(klon,klev,2)) - - IF (config_inca /= 'none') THEN - ALLOCATE(topswad_inca(klon), solswad_inca(klon)) - ALLOCATE(topswad0_inca(klon), solswad0_inca(klon)) - ALLOCATE(topswai_inca(klon), solswai_inca(klon)) - ALLOCATE(topsw_inca(klon,9), solsw_inca(klon,9)) - ALLOCATE(topsw0_inca(klon,9), solsw0_inca(klon,9)) - END IF - ! Following 4 variables are needed only by INCA but must be - ! allocated as they exist in the phytrac argument list - ALLOCATE(tau_inca(klon,klev,9,2)) - ALLOCATE(piz_inca(klon,klev,9,2)) - ALLOCATE(cg_inca(klon,klev,9,2)) + ALLOCATE(ccm(klon,klev,2)) @@ -505,6 +461,4 @@ deallocate(ibas_con, itop_con) deallocate(rain_con, snow_con) -! - deallocate(sulfate_pi) deallocate(rlonPOS) deallocate(newsst) @@ -513,17 +467,4 @@ deallocate(topswai, solswai) - deallocate(tau_ae, piz_ae) - deallocate(cg_ae) - - IF (config_inca /= 'none') THEN - deallocate(topswad_inca, solswad_inca) - deallocate(topswad0_inca, solswad0_inca) - deallocate(topswai_inca, solswai_inca) - deallocate(topsw_inca, solsw_inca) - deallocate(topsw0_inca, solsw0_inca) - END IF - deallocate(tau_inca) - deallocate(piz_inca) - deallocate(cg_inca) deallocate(ccm) Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F (revision 1149) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/physiq.F (revision 1150) @@ -50,5 +50,5 @@ c c nlon----input-I-nombre de points horizontaux -c nlev----input-I-nombre de couches verticales +c nlev----input-I-nombre de couches verticales, doit etre egale a klev c debut---input-L-variable logique indiquant le premier passage c lafin---input-L-variable logique indiquant le dernier passage @@ -736,5 +736,4 @@ EXTERNAL phyetat0 ! lire l'etat initial de la physique EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique - EXTERNAL radlwsw ! rayonnements solaire et infrarouge EXTERNAL suphel ! initialiser certaines constantes EXTERNAL transp ! transport total de l'eau et de l'energie @@ -1056,7 +1055,4 @@ CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max CHARACTER*40 tinst, tave, typeval -cjq Aerosol effects (Johannes Quaas, 27/11/2003) - REAL sulfate(klon, klev) ! SO4 aerosol concentration [ug/m3] - REAL cldtaupi(klon,klev) ! Cloud optical thickness for pre-industrial (pi) aerosols @@ -1067,9 +1063,16 @@ ! Aerosol optical properties - - ! Aerosol optical properties by INCA model - CHARACTER*4 :: rfname(9) - REAL aerindex(klon) ! POLDER aerosol index - + 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 + + ! Parameters LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not @@ -1080,5 +1083,9 @@ ! false : lecture des aerosol dans un fichier c$OMP THREADPRIVATE(aerosol_couple) - + INTEGER, SAVE :: flag_aerosol +c$OMP THREADPRIVATE(flag_aerosol) + LOGICAL, SAVE :: new_aod +c$OMP THREADPRIVATE(new_aod) + c c Declaration des constantes et des fonctions thermodynamiques @@ -1106,15 +1113,15 @@ write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!' write(lunout,*) - s 'nlon,nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys' + s 'nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys' write(lunout,*) - s nlon,nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys + s nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys write(lunout,*) 'papers, play, phi, u, v, t, omega' - do k=1,nlev + do k=1,klev write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), s u(igout,k),v(igout,k),t(igout,k),omega(igout,k) enddo write(lunout,*) 'ovap (g/kg), oliq (g/kg)' - do k=1,nlev + do k=1,klev write(lunout,*) qx(igout,k,1)*1000,qx(igout,k,2)*1000. enddo @@ -1190,7 +1197,4 @@ u10m(:,:)=0. v10m(:,:)=0. - piz_ae(:,:,:)=0. - tau_ae(:,:,:)=0. - cg_ae(:,:,:)=0. rain_con(:)=0. snow_con(:)=0. @@ -1205,20 +1209,6 @@ wmax_th(:)=0. tau_overturning_th(:)=0. - IF (config_inca /= 'none') THEN - tau_inca(:,:,:,:) = 0. - piz_inca(:,:,:,:) = 0. - cg_inca(:,:,:,:) = 0. - ccm(:,:,:) = 0. - topswai_inca(:) = 0. - topswad_inca(:) = 0. - topswad0_inca(:) = 0. - topsw_inca(:,:) = 0. - topsw0_inca(:,:) = 0. - solswai_inca(:) = 0. - solswad_inca(:) = 0. - solswad0_inca(:) = 0. - solsw_inca(:,:) = 0. - solsw0_inca(:,:) = 0. - END IF + + IF (config_inca /= 'none') ccm(:,:,:) = 0. rnebcon0(:,:) = 0.0 @@ -1239,4 +1229,5 @@ . iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut, . ok_ade, ok_aie, aerosol_couple, + . flag_aerosol, new_aod, . bl95_b0, bl95_b1, . iflag_thermals,nsplit_thermals,tau_thermals, @@ -2658,17 +2649,14 @@ cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr) IF (ok_ade.OR.ok_aie) THEN - IF ( .NOT. aerosol_couple ) THEN - ! Get sulfate aerosol distribution - CALL readsulfate(rjourvrai, debut, sulfate) - CALL readsulfate_preind(rjourvrai, debut, sulfate_pi) - - ! Calculate aerosol optical properties (Olivier Boucher) - CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl, - . tau_ae, piz_ae, cg_ae, aerindex) - ENDIF + IF (.NOT. aerosol_couple) + & CALL aerosol_optic( + & debut, new_aod, flag_aerosol, rjourvrai, pdtphys, + & pplay, paprs, t_seri, rhcl, + & maerosol, maerosol_pi, + & tau_aero, piz_aero, cg_aero ) ELSE - tau_ae(:,:,:)=0.0 - piz_ae(:,:,:)=0.0 - cg_ae(:,:,:)=0.0 + tau_aero(:,:,:,:) = 0. + piz_aero(:,:,:,:) = 0. + cg_aero(:,:,:,:) = 0. ENDIF @@ -2847,7 +2835,7 @@ IF (aerosol_couple) THEN - sulfate(:,:) = ccm(:,:,1) - sulfate_pi(:,:) = ccm(:,:,2) - ENDIF + maerosol(:,:) = ccm(:,:,1) + maerosol_pi(:,:) = ccm(:,:,2) + END IF if (ok_newmicro) then @@ -2857,5 +2845,5 @@ . flwp, fiwp, flwc, fiwc, e ok_aie, - e sulfate, sulfate_pi, + e maerosol, maerosol_pi, e bl95_b0, bl95_b1, s cldtaupi, re, fl) @@ -2865,5 +2853,5 @@ . cldh, cldl, cldm, cldt, cldq, e ok_aie, - e sulfate, sulfate_pi, + e maerosol, maerosol_pi, e bl95_b0, bl95_b1, s cldtaupi, re, fl) @@ -2895,41 +2883,48 @@ IF (aerosol_couple) THEN #ifdef INCA - CALL radlwsw_inca - e (kdlon,kflev,dist, rmu0, fract, solaire, - e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, - e wo, - e cldfra, cldemi, cldtau, - s heat,heat0,cool,cool0,radsol,albpla, - s topsw,toplw,solsw,sollw, - s sollwdown, - s topsw0,toplw0,solsw0,sollw0, - s lwdn0, lwdn, lwup0, lwup, - s swdn0, swdn, swup0, swup, - e ok_ade, ok_aie, - e tau_inca, piz_inca, cg_inca, - s topswad_inca, solswad_inca, - s topswad0_inca, solswad0_inca, - s topsw_inca, topsw0_inca, - s solsw_inca, solsw0_inca, - e cldtaupi, - s topswai_inca, solswai_inca) + CALL radlwsw_inca + e (kdlon,kflev,dist, rmu0, fract, solaire, + e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, + e wo, + e cldfra, cldemi, cldtau, + s heat,heat0,cool,cool0,radsol,albpla, + s topsw,toplw,solsw,sollw, + s sollwdown, + s topsw0,toplw0,solsw0,sollw0, + s lwdn0, lwdn, lwup0, lwup, + s swdn0, swdn, swup0, swup, + e ok_ade, ok_aie, + e tau_aero, piz_aero, cg_aero, + s topswad_aero, solswad_aero, + s topswad0_aero, solswad0_aero, + s topsw_aero, topsw0_aero, + s solsw_aero, solsw0_aero, + e cldtaupi, + s topswai_aero, solswai_aero) + #endif ELSE - CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS) - e (dist, rmu0, fract, - e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, - e wo, - e cldfra, cldemi, cldtau, - s heat,heat0,cool,cool0,radsol,albpla, - s topsw,toplw,solsw,sollw, - s sollwdown, - s topsw0,toplw0,solsw0,sollw0, - s lwdn0, lwdn, lwup0, lwup, - s swdn0, swdn, swup0, swup, - e ok_ade, ok_aie, ! new for aerosol radiative effects - e tau_ae, piz_ae, cg_ae, ! ="= - s topswad, solswad, ! ="= - e cldtaupi, ! ="= - s topswai, solswai,zqsat,flwc,fiwc) ! ="= + + CALL radlwsw_aero + e (dist, rmu0, fract, solaire, + e paprs, pplay,zxtsol,albsol1, albsol2, + e t_seri,q_seri,wo, + e cldfra, cldemi, cldtau, + e ok_ade, ok_aie, + e tau_aero, piz_aero, cg_aero, + e cldtaupi,new_aod, + s heat,heat0,cool,cool0,radsol,albpla, + s topsw,toplw,solsw,sollw, + s sollwdown, + s topsw0,toplw0,solsw0,sollw0, + s lwdn0, lwdn, lwup0, lwup, + s swdn0, swdn, swup0, swup, + s topswad_aero, solswad_aero, + s topswai_aero, solswai_aero, + o topswad0_aero, solswad0_aero, + o topsw_aero, topsw0_aero, + o solsw_aero, solsw0_aero) + + ENDIF itaprad = 0 @@ -3159,5 +3154,5 @@ I lafin, I nlon, - I nlev, + I klev, I dtime, I u, @@ -3207,7 +3202,7 @@ I aerosol_couple, I flxmass_w, - I tau_inca, - I piz_inca, - I cg_inca, + I tau_aero, + I piz_aero, + I cg_aero, I ccm, I rfname, @@ -3218,5 +3213,5 @@ print*,'Attention on met a 0 les thermiques pour phystoke' call phystokenc ( - I nlon,nlev,pdtphys,rlon,rlat, + I nlon,klev,pdtphys,rlon,rlat, I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, I fm_therm,entr_therm, @@ -3415,11 +3410,11 @@ write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!' write(lunout,*) - s 'nlon,nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos' + s 'nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos' write(lunout,*) - s nlon,nlev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys, + s nlon,klev,nqtot,debut,lafin,rjourvrai,gmtime,pdtphys, s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), s pctsrf(igout,is_sic) write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva' - do k=1,nlev + do k=1,klev write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), @@ -3427,10 +3422,10 @@ enddo write(lunout,*) 'cool,heat' - do k=1,nlev + do k=1,klev write(lunout,*) cool(igout,k),heat(igout,k) enddo write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec' - do k=1,nlev + do k=1,klev write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k) @@ -3439,5 +3434,5 @@ write(lunout,*) 'd_ps ',d_ps(igout) write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 ' - do k=1,nlev + do k=1,klev write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), s d_qx(igout,k,1),d_qx(igout,k,2) Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F (revision 1149) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/phytrac.F (revision 1150) @@ -57,7 +57,7 @@ I aerosol_couple, I flxmass_w, - I tau_inca, - I piz_inca, - I cg_inca, + I tau_aero, + I piz_aero, + I cg_aero, I ccm, I rfname, @@ -138,7 +138,7 @@ CHARACTER(len=8) :: solsym(nbtr) integer la - REAL :: tau_inca(klon,klev,9,2) - REAL :: piz_inca(klon,klev,9,2) - REAL :: cg_inca(klon,klev,9,2) + 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) @@ -458,7 +458,7 @@ $ t_seri, ! for chimiaq $ rh, - $ tau_inca, - $ piz_inca, - $ cg_inca, + $ tau_aero, + $ piz_aero, + $ cg_aero, $ rfname, $ ccm, @@ -497,7 +497,4 @@ $ sh, !sh $ rh, !rh - $ .false., !wrhstts - $ hbuf, !hbuf - $ obuf, !obuf $ iip1, !nx $ jjp1, !ny Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/radlwsw.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/radlwsw.F (revision 1149) +++ (revision ) @@ -1,401 +1,0 @@ -! -! $Header$ -! - SUBROUTINE radlwsw(dist, rmu0, fract, - . paprs, pplay,tsol,alb1, alb2, t,q,wo, - . cldfra, cldemi, cldtaupd, - . heat,heat0,cool,cool0,radsol,albpla, - . topsw,toplw,solsw,sollw, - . sollwdown, - . topsw0,toplw0,solsw0,sollw0, - . lwdn0, lwdn, lwup0, lwup, - . swdn0, swdn, swup0, swup, - . ok_ade, ok_aie, - . tau_ae, piz_ae, cg_ae, - . topswad, solswad, - . cldtaupi, topswai, solswai,qsat,flwc,fiwc) -c - USE dimphy - IMPLICIT none -c====================================================================== -c Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719 -c Objet: interface entre le modele et les rayonnements -c Arguments: -c dist-----input-R- distance astronomique terre-soleil -c rmu0-----input-R- cosinus de l'angle zenithal -c fract----input-R- duree d'ensoleillement normalisee -c co2_ppm--input-R- concentration du gaz carbonique (en ppm) -c solaire--input-R- constante solaire (W/m**2) -c paprs----input-R- pression a inter-couche (Pa) -c pplay----input-R- pression au milieu de couche (Pa) -c tsol-----input-R- temperature du sol (en K) -c alb1-----input-R- albedo du sol(entre 0 et 1) dans l'interval visible -c alb2-----input-R- albedo du sol(entre 0 et 1) dans l'interval proche infra-rouge -c t--------input-R- temperature (K) -c q--------input-R- vapeur d'eau (en kg/kg) -c wo-------input-R- contenu en ozone (en kg/kg) correction MPL 100505 -c cldfra---input-R- fraction nuageuse (entre 0 et 1) -c cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value) -c cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1) -c ok_ade---input-L- apply the Aerosol Direct Effect or not? -c ok_aie---input-L- apply the Aerosol Indirect Effect or not? -c tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F) -c cldtaupi-input-R- epaisseur optique des nuages dans le visible -c calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller -c droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd -c it is needed for the diagnostics of the aerosol indirect radiative forcing -c -c heat-----output-R- echauffement atmospherique (visible) (K/jour) -c cool-----output-R- refroidissement dans l'IR (K/jour) -c radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas) -c albpla---output-R- albedo planetaire (entre 0 et 1) -c topsw----output-R- flux solaire net au sommet de l'atm. -c toplw----output-R- ray. IR montant au sommet de l'atmosphere -c solsw----output-R- flux solaire net a la surface -c sollw----output-R- ray. IR montant a la surface -c solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir) -c topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir) -c solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind) -c topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind) -c -c ATTENTION: swai and swad have to be interpreted in the following manner: -c --------- -c ok_ade=F & ok_aie=F -both are zero -c ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad -c indirect is zero -c ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai -c direct is zero -c ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai -c aerosol direct forcing is F_{AD} = topswai-topswad -c - -c====================================================================== -cym#include "dimensions.h" -cym#include "dimphy.h" -cym#include "raddim.h" -#include "YOETHF.h" -c - real rmu0(klon), fract(klon), dist -cIM real co2_ppm -cIM real solaire -#include "clesphys.h" -c - real paprs(klon,klev+1), pplay(klon,klev) - real alb1(klon), alb2(klon), tsol(klon) - real t(klon,klev), q(klon,klev), wo(klon,klev) - real cldfra(klon,klev), cldemi(klon,klev), cldtaupd(klon,klev) - real heat(klon,klev), cool(klon,klev) - real heat0(klon,klev), cool0(klon,klev) - real radsol(klon), topsw(klon), toplw(klon) - real solsw(klon), sollw(klon), albpla(klon) - real topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon) - real sollwdown(klon) -cIM output 3D - REAL*8 ZFSUP(KDLON,KFLEV+1) - REAL*8 ZFSDN(KDLON,KFLEV+1) - REAL*8 ZFSUP0(KDLON,KFLEV+1) - REAL*8 ZFSDN0(KDLON,KFLEV+1) -c - REAL*8 ZFLUP(KDLON,KFLEV+1) - REAL*8 ZFLDN(KDLON,KFLEV+1) - REAL*8 ZFLUP0(KDLON,KFLEV+1) - REAL*8 ZFLDN0(KDLON,KFLEV+1) -c - REAL*8 zx_alpha1, zx_alpha2 -c -#include "YOMCST.h" -c - INTEGER k, kk, i, j, iof, nb_gr - EXTERNAL LW_LMDAR4,SW_LMDAR4 -c -cIM ctes ds clesphys.h REAL*8 RCO2, RCH4, RN2O, RCFC11, RCFC12 - REAL*8 PSCT -c - REAL*8 PALBD(kdlon,2), PALBP(kdlon,2) - REAL*8 PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon) - REAL*8 PPSOL(kdlon), PDP(kdlon,klev) - REAL*8 PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1) - REAL*8 PTAVE(kdlon,kflev) - REAL*8 PWV(kdlon,kflev), PQS(kdlon,kflev), POZON(kdlon,kflev) - REAL*8 PAER(kdlon,kflev,5) - REAL*8 PCLDLD(kdlon,kflev) - REAL*8 PCLDLU(kdlon,kflev) - REAL*8 PCLDSW(kdlon,kflev) - REAL*8 PTAU(kdlon,2,kflev) - REAL*8 POMEGA(kdlon,2,kflev) - REAL*8 PCG(kdlon,2,kflev) -c - REAL*8 zfract(kdlon), zrmu0(kdlon), zdist -c - REAL*8 zheat(kdlon,kflev), zcool(kdlon,kflev) - REAL*8 zheat0(kdlon,kflev), zcool0(kdlon,kflev) - REAL*8 ztopsw(kdlon), ztoplw(kdlon) - REAL*8 zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon) -cIM - REAL*8 zsollwdown(kdlon) -c - REAL*8 ztopsw0(kdlon), ztoplw0(kdlon) - REAL*8 zsolsw0(kdlon), zsollw0(kdlon) - REAL*8 zznormcp -cIM output 3D : SWup, SWdn, LWup, LWdn - REAL swdn(klon,kflev+1),swdn0(klon,kflev+1) - REAL swup(klon,kflev+1),swup0(klon,kflev+1) - REAL lwdn(klon,kflev+1),lwdn0(klon,kflev+1) - REAL lwup(klon,kflev+1),lwup0(klon,kflev+1) - REAL qsat(klon,klev),flwc(klon,klev),fiwc(klon,klev) -c-OB -cjq the following quantities are needed for the aerosol radiative forcings - - real topswad(klon), solswad(klon) ! output: aerosol direct forcing at TOA and surface - real topswai(klon), solswai(klon) ! output: aerosol indirect forcing atTOA and surface - real tau_ae(klon,klev,2), piz_ae(klon,klev,2), cg_ae(klon,klev,2) ! aerosol optical properties (see aeropt.F) - real cldtaupi(klon,klev) ! cloud optical thickness for pre-industrial aerosol concentrations - ! (i.e., with a smaller droplet concentrationand thus larger droplet radii) - logical ok_ade, ok_aie ! switches whether to use aerosol direct (indirect) effects or not - real*8 tauae(kdlon,kflev,2) ! aer opt properties - real*8 pizae(kdlon,kflev,2) - real*8 cgae(kdlon,kflev,2) - REAL*8 PTAUA(kdlon,2,kflev) ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use - REAL*8 POMEGAA(kdlon,2,kflev) ! dito for single scatt albedo - REAL*8 ztopswad(kdlon), zsolswad(kdlon) ! Aerosol direct forcing at TOAand surface - REAL*8 ztopswai(kdlon), zsolswai(kdlon) ! dito, indirect -cjq-end -!rv - tauae(:,:,:)=0. - pizae(:,:,:)=0. - cgae(:,:,:)=0. -!rv - -c -c------------------------------------------- - nb_gr = klon / kdlon - IF (nb_gr*kdlon .NE. klon) THEN - PRINT*, "kdlon mauvais:", klon, kdlon, nb_gr - CALL abort - ENDIF - IF (kflev .NE. klev) THEN - PRINT*, "kflev differe de klev, kflev, klev" - CALL abort - ENDIF -c------------------------------------------- - DO k = 1, klev - DO i = 1, klon - heat(i,k)=0. - cool(i,k)=0. - heat0(i,k)=0. - cool0(i,k)=0. - ENDDO - ENDDO -c - zdist = dist -c -cIM anciennes valeurs -c RCO2 = co2_ppm * 1.0e-06 * 44.011/28.97 -c -cIM : on met RCO2, RCH4, RN2O, RCFC11 et RCFC12 dans clesphys.h /lecture ds conf_phys.F90 -c RCH4 = 1.65E-06* 16.043/28.97 -c RN2O = 306.E-09* 44.013/28.97 -c RCFC11 = 280.E-12* 137.3686/28.97 -c RCFC12 = 484.E-12* 120.9140/28.97 -cIM anciennes valeurs -c RCH4 = 1.72E-06* 16.043/28.97 -c RN2O = 310.E-09* 44.013/28.97 -c -c PRINT*,'IMradlwsw : solaire, co2= ', solaire, co2_ppm - PSCT = solaire/zdist/zdist -c - DO 99999 j = 1, nb_gr - iof = kdlon*(j-1) -c - DO i = 1, kdlon - zfract(i) = fract(iof+i) - zrmu0(i) = rmu0(iof+i) - PALBD(i,1) = alb1(iof+i) -! PALBD(i,2) = alb1(iof+i) - PALBD(i,2) = alb2(iof+i) - PALBP(i,1) = alb1(iof+i) -! PALBP(i,2) = alb1(iof+i) - PALBP(i,2) = alb2(iof+i) -cIM cf. JLD pour etre en accord avec ORCHIDEE il faut mettre PEMIS(i) = 0.96 - PEMIS(i) = 1.0 - PVIEW(i) = 1.66 - PPSOL(i) = paprs(iof+i,1) - zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2)) - . / (pplay(iof+i,1)-pplay(iof+i,2)) - zx_alpha2 = 1.0 - zx_alpha1 - PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2 - PTL(i,klev+1) = t(iof+i,klev) - PDT0(i) = tsol(iof+i) - PTL(i,1) - ENDDO - DO k = 2, kflev - DO i = 1, kdlon - PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5 - ENDDO - ENDDO - DO k = 1, kflev - DO i = 1, kdlon - PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1) - PTAVE(i,k) = t(iof+i,k) - PWV(i,k) = MAX (q(iof+i,k), 1.0e-12) - PQS(i,k) = PWV(i,k) -c wo: cm.atm (epaisseur en cm dans la situation standard) -c POZON: kg/kg - POZON(i,k) = MAX(wo(iof+i,k),1.0e-12)*RG/46.6968 - . /(paprs(iof+i,k)-paprs(iof+i,k+1)) - . *(paprs(iof+i,1)/101325.0) - PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k) - PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k) - PCLDSW(i,k) = cldfra(iof+i,k) - PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable - PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines - POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k)) - POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k)) - PCG(i,1,k) = 0.865 - PCG(i,2,k) = 0.910 -c-OB -cjq Introduced for aerosol indirect forcings. -cjq The following values use the cloud optical thickness calculated from -cjq present-day aerosol concentrations whereas the quantities without the -cjq "A" at the end are for pre-industial (natural-only) aerosol concentrations -cjq - PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable - PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines - POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k)) - POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k)) -cjq-end - ENDDO - ENDDO -c - DO k = 1, kflev+1 - DO i = 1, kdlon - PPMB(i,k) = paprs(iof+i,k)/100.0 - ENDDO - ENDDO -c - DO kk = 1, 5 - DO k = 1, kflev - DO i = 1, kdlon - PAER(i,k,kk) = 1.0E-15 - ENDDO - ENDDO - ENDDO -c-OB - DO k = 1, kflev - DO i = 1, kdlon - tauae(i,k,1)=tau_ae(iof+i,k,1) - pizae(i,k,1)=piz_ae(iof+i,k,1) - cgae(i,k,1) =cg_ae(iof+i,k,1) - tauae(i,k,2)=tau_ae(iof+i,k,2) - pizae(i,k,2)=piz_ae(iof+i,k,2) - cgae(i,k,2) =cg_ae(iof+i,k,2) - ENDDO - ENDDO -c -c===== si iflag_rrtm=0 ================================================ -cIM ctes ds clesphys.h CALL LW(RCO2,RCH4,RN2O,RCFC11,RCFC12, -cIM ctes ds clesphys.h CALL SW(PSCT, RCO2, zrmu0, zfract, -c - if (iflag_rrtm.eq.0) then - CALL LW_LMDAR4( - . PPMB, PDP, - . PPSOL,PDT0,PEMIS, - . PTL, PTAVE, PWV, POZON, PAER, - . PCLDLD,PCLDLU, - . PVIEW, - . zcool, zcool0, - . ztoplw,zsollw,ztoplw0,zsollw0, - . zsollwdown, - . ZFLUP, ZFLDN, ZFLUP0,ZFLDN0) - CALL SW_LMDAR4(PSCT, zrmu0, zfract, - S PPMB, PDP, - S PPSOL, PALBD, PALBP, - S PTAVE, PWV, PQS, POZON, PAER, - S PCLDSW, PTAU, POMEGA, PCG, - S zheat, zheat0, - S zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0, - S ZFSUP,ZFSDN,ZFSUP0,ZFSDN0, - S tauae, pizae, cgae, ! aerosol optical properties - s PTAUA, POMEGAA, - s ztopswad,zsolswad,ztopswai,zsolswai, ! diagnosed aerosol forcing - J ok_ade, ok_aie) ! apply aerosol effects or not? - else -c===== si iflag_rrtm=1, on passe dans SW via RECMWFL =============== - PRINT*, "Cette option ne fonctionne pas encore !!!" - CALL abort - endif ! if(iflag_rrtm=0) - -c====================================================================== - DO i = 1, kdlon - radsol(iof+i) = zsolsw(i) + zsollw(i) - topsw(iof+i) = ztopsw(i) - toplw(iof+i) = ztoplw(i) - solsw(iof+i) = zsolsw(i) - sollw(iof+i) = zsollw(i) - sollwdown(iof+i) = zsollwdown(i) -cIM - DO k = 1, kflev+1 - lwdn0 ( iof+i,k) = ZFLDN0 ( i,k) - lwdn ( iof+i,k) = ZFLDN ( i,k) - lwup0 ( iof+i,k) = ZFLUP0 ( i,k) - lwup ( iof+i,k) = ZFLUP ( i,k) - ENDDO -c - topsw0(iof+i) = ztopsw0(i) - toplw0(iof+i) = ztoplw0(i) - solsw0(iof+i) = zsolsw0(i) - sollw0(iof+i) = zsollw0(i) - albpla(iof+i) = zalbpla(i) -cIM - DO k = 1, kflev+1 - swdn0 ( iof+i,k) = ZFSDN0 ( i,k) - swdn ( iof+i,k) = ZFSDN ( i,k) - swup0 ( iof+i,k) = ZFSUP0 ( i,k) - swup ( iof+i,k) = ZFSUP ( i,k) - ENDDO !k=1, kflev+1 - ENDDO -cjq-transform the aerosol forcings, if they have -cjq to be calculated - IF (ok_ade) THEN - DO i = 1, kdlon - topswad(iof+i) = ztopswad(i) - solswad(iof+i) = zsolswad(i) - ENDDO - ELSE - DO i = 1, kdlon - topswad(iof+i) = 0.0 - solswad(iof+i) = 0.0 - ENDDO - ENDIF - IF (ok_aie) THEN - DO i = 1, kdlon - topswai(iof+i) = ztopswai(i) - solswai(iof+i) = zsolswai(i) - ENDDO - ELSE - DO i = 1, kdlon - topswai(iof+i) = 0.0 - solswai(iof+i) = 0.0 - ENDDO - ENDIF -cjq-end - DO k = 1, kflev -c DO i = 1, kdlon -c heat(iof+i,k) = zheat(i,k) -c cool(iof+i,k) = zcool(i,k) -c heat0(iof+i,k) = zheat0(i,k) -c cool0(iof+i,k) = zcool0(i,k) -c ENDDO - DO i = 1, kdlon -C scale factor to take into account the difference between -C dry air and watter vapour scpecific heat capacity - zznormcp=1.0+RVTMP2*PWV(i,k) - heat(iof+i,k) = zheat(i,k)/zznormcp - cool(iof+i,k) = zcool(i,k)/zznormcp - heat0(iof+i,k) = zheat0(i,k)/zznormcp - cool0(iof+i,k) = zcool0(i,k)/zznormcp - ENDDO - ENDDO -c -99999 CONTINUE - RETURN - END Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/radlwsw_aero.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/radlwsw_aero.F90 (revision 1150) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/radlwsw_aero.F90 (revision 1150) @@ -0,0 +1,413 @@ +SUBROUTINE radlwsw_aero( & + dist, rmu0, fract, solaire, & + paprs, pplay,tsol,albedo, alblw, & + t,q,wo,& + cldfra, cldemi, cldtaupd,& + ok_ade, ok_aie,& + tau_aero, piz_aero, cg_aero,& + cldtaupi, new_aod, & + heat,heat0,cool,cool0,radsol,albpla,& + topsw,toplw,solsw,sollw,& + sollwdown,& + topsw0,toplw0,solsw0,sollw0,& + lwdn0, lwdn, lwup0, lwup,& + swdn0, swdn, swup0, swup,& + topswad_aero, solswad_aero,& + topswai_aero, solswai_aero, & + topswad0_aero, solswad0_aero,& + topsw_aero, topsw0_aero,& + solsw_aero, solsw0_aero) + + + + USE DIMPHY + + IMPLICIT NONE + !====================================================================== + ! Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719 + ! Objet: interface entre le modele et les rayonnements + ! Arguments: + ! dist-----input-R- distance astronomique terre-soleil + ! rmu0-----input-R- cosinus de l'angle zenithal + ! fract----input-R- duree d'ensoleillement normalisee + ! co2_ppm--input-R- concentration du gaz carbonique (en ppm) + ! solaire--input-R- constante solaire (W/m**2) + ! paprs----input-R- pression a inter-couche (Pa) + ! pplay----input-R- pression au milieu de couche (Pa) + ! tsol-----input-R- temperature du sol (en K) + ! albedo---input-R- albedo du sol (entre 0 et 1) + ! t--------input-R- temperature (K) + ! q--------input-R- vapeur d'eau (en kg/kg) + ! wo-------input-R- contenu en ozone (en kg/kg) correction MPL 100505 + ! cldfra---input-R- fraction nuageuse (entre 0 et 1) + ! cldtaupd---input-R- epaisseur optique des nuages dans le visible (present-day value) + ! cldemi---input-R- emissivite des nuages dans l'IR (entre 0 et 1) + ! ok_ade---input-L- apply the Aerosol Direct Effect or not? + ! ok_aie---input-L- apply the Aerosol Indirect Effect or not? + ! tau_ae, piz_ae, cg_ae-input-R- aerosol optical properties (calculated in aeropt.F) + ! cldtaupi-input-R- epaisseur optique des nuages dans le visible + ! calculated for pre-industrial (pi) aerosol concentrations, i.e. with smaller + ! droplet concentration, thus larger droplets, thus generally cdltaupi cldtaupd + ! it is needed for the diagnostics of the aerosol indirect radiative forcing + ! + ! heat-----output-R- echauffement atmospherique (visible) (K/jour) + ! cool-----output-R- refroidissement dans l'IR (K/jour) + ! radsol---output-R- bilan radiatif net au sol (W/m**2) (+ vers le bas) + ! albpla---output-R- albedo planetaire (entre 0 et 1) + ! topsw----output-R- flux solaire net au sommet de l'atm. + ! toplw----output-R- ray. IR montant au sommet de l'atmosphere + ! solsw----output-R- flux solaire net a la surface + ! sollw----output-R- ray. IR montant a la surface + ! solswad---output-R- ray. solaire net absorbe a la surface (aerosol dir) + ! topswad---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol dir) + ! solswai---output-R- ray. solaire net absorbe a la surface (aerosol ind) + ! topswai---output-R- ray. solaire absorbe au sommet de l'atm. (aerosol ind) + ! + ! ATTENTION: swai and swad have to be interpreted in the following manner: + ! --------- + ! ok_ade=F & ok_aie=F -both are zero + ! ok_ade=T & ok_aie=F -aerosol direct forcing is F_{AD} = topsw-topswad + ! indirect is zero + ! ok_ade=F & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai + ! direct is zero + ! ok_ade=T & ok_aie=T -aerosol indirect forcing is F_{AI} = topsw-topswai + ! aerosol direct forcing is F_{AD} = topswai-topswad + ! + + !====================================================================== + + ! ==================================================================== + ! Adapte au modele de chimie INCA par Celine Deandreis & Anne Cozic -- 2009 + ! 1 = ZERO + ! 2 = AER total + ! 3 = NAT + ! 4 = BC + ! 5 = SO4 + ! 6 = POM + ! 7 = DUST + ! 8 = SS + ! 9 = NO3 + ! + ! ==================================================================== + include "YOETHF.h" + include "YOMCST.h" + +! Input arguments + REAL, INTENT(in) :: solaire + REAL, INTENT(in) :: dist + REAL, INTENT(in) :: rmu0(KLON), fract(KLON) + REAL, INTENT(in) :: paprs(KLON,KLEV+1), pplay(KLON,KLEV) + REAL, INTENT(in) :: albedo(KLON), alblw(KLON), tsol(KLON) + REAL, INTENT(in) :: t(KLON,KLEV), q(KLON,KLEV), wo(KLON,KLEV) + LOGICAL, INTENT(in) :: ok_ade, ok_aie ! switches whether to use aerosol direct (indirect) effects or not + REAL, INTENT(in) :: cldfra(KLON,KLEV), cldemi(KLON,KLEV), cldtaupd(KLON,KLEV) + REAL, INTENT(in) :: tau_aero(KLON,KLEV,9,2) ! aerosol optical properties (see aeropt.F) + REAL, INTENT(in) :: piz_aero(KLON,KLEV,9,2) ! aerosol optical properties (see aeropt.F) + REAL, INTENT(in) :: cg_aero(KLON,KLEV,9,2) ! aerosol optical properties (see aeropt.F) + REAL, INTENT(in) :: cldtaupi(KLON,KLEV) ! cloud optical thickness for pre-industrial aerosol concentrations + LOGICAL, INTENT(in) :: new_aod ! flag pour retrouver les resultats exacts de l'AR4 dans le cas ou l'on ne travaille qu'avec les sulfates + +! Output arguments + REAL, INTENT(out) :: heat(KLON,KLEV), cool(KLON,KLEV) + REAL, INTENT(out) :: heat0(KLON,KLEV), cool0(KLON,KLEV) + REAL, INTENT(out) :: radsol(KLON), topsw(KLON), toplw(KLON) + REAL, INTENT(out) :: solsw(KLON), sollw(KLON), albpla(KLON) + REAL, INTENT(out) :: topsw0(KLON), toplw0(KLON), solsw0(KLON), sollw0(KLON) + REAL, INTENT(out) :: sollwdown(KLON) + REAL, INTENT(out) :: swdn(KLON,kflev+1),swdn0(KLON,kflev+1) + REAL, INTENT(out) :: swup(KLON,kflev+1),swup0(KLON,kflev+1) + REAL, INTENT(out) :: lwdn(KLON,kflev+1),lwdn0(KLON,kflev+1) + REAL, INTENT(out) :: lwup(KLON,kflev+1),lwup0(KLON,kflev+1) + REAL, INTENT(out) :: topswad_aero(KLON), solswad_aero(KLON) ! output: aerosol direct forcing at TOA and surface + REAL, INTENT(out) :: topswai_aero(KLON), solswai_aero(KLON) ! output: aerosol indirect forcing atTOA and surface + REAL, DIMENSION(klon), INTENT(out) :: topswad0_aero + REAL, DIMENSION(klon), INTENT(out) :: solswad0_aero + REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw_aero + REAL, DIMENSION(kdlon,9), INTENT(out) :: topsw0_aero + REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw_aero + REAL, DIMENSION(kdlon,9), INTENT(out) :: solsw0_aero + +! Local variables + REAL*8 ZFSUP(KDLON,KFLEV+1) + REAL*8 ZFSDN(KDLON,KFLEV+1) + REAL*8 ZFSUP0(KDLON,KFLEV+1) + REAL*8 ZFSDN0(KDLON,KFLEV+1) + REAL*8 ZFLUP(KDLON,KFLEV+1) + REAL*8 ZFLDN(KDLON,KFLEV+1) + REAL*8 ZFLUP0(KDLON,KFLEV+1) + REAL*8 ZFLDN0(KDLON,KFLEV+1) + REAL*8 zx_alpha1, zx_alpha2 + INTEGER k, kk, i, j, iof, nb_gr + REAL*8 PSCT + REAL*8 PALBD(kdlon,2), PALBP(kdlon,2) + REAL*8 PEMIS(kdlon), PDT0(kdlon), PVIEW(kdlon) + REAL*8 PPSOL(kdlon), PDP(kdlon,KLEV) + REAL*8 PTL(kdlon,kflev+1), PPMB(kdlon,kflev+1) + REAL*8 PTAVE(kdlon,kflev) + REAL*8 PWV(kdlon,kflev), PQS(kdlon,kflev), POZON(kdlon,kflev) + REAL*8 PAER(kdlon,kflev,5) + REAL*8 PCLDLD(kdlon,kflev) + REAL*8 PCLDLU(kdlon,kflev) + REAL*8 PCLDSW(kdlon,kflev) + REAL*8 PTAU(kdlon,2,kflev) + REAL*8 POMEGA(kdlon,2,kflev) + REAL*8 PCG(kdlon,2,kflev) + REAL*8 zfract(kdlon), zrmu0(kdlon), zdist + REAL*8 zheat(kdlon,kflev), zcool(kdlon,kflev) + REAL*8 zheat0(kdlon,kflev), zcool0(kdlon,kflev) + REAL*8 ztopsw(kdlon), ztoplw(kdlon) + REAL*8 zsolsw(kdlon), zsollw(kdlon), zalbpla(kdlon) + REAL*8 zsollwdown(kdlon) + REAL*8 ztopsw0(kdlon), ztoplw0(kdlon) + REAL*8 zsolsw0(kdlon), zsollw0(kdlon) + REAL*8 zznormcp + REAL*8 tauaero(kdlon,kflev,9,2) ! aer opt properties + REAL*8 pizaero(kdlon,kflev,9,2) + REAL*8 cgaero(kdlon,kflev,9,2) + REAL*8 PTAUA(kdlon,2,kflev) ! present-day value of cloud opt thickness (PTAU is pre-industrial value), local use + REAL*8 POMEGAA(kdlon,2,kflev) ! dito for single scatt albedo + REAL*8 ztopswadaero(kdlon), zsolswadaero(kdlon) ! Aerosol direct forcing at TOAand surface + REAL*8 ztopswad0aero(kdlon), zsolswad0aero(kdlon) ! Aerosol direct forcing at TOAand surface + REAL*8 ztopswaiaero(kdlon), zsolswaiaero(kdlon) ! dito, indirect + REAL*8 ztopsw_aero(kdlon,9), ztopsw0_aero(kdlon,9) + REAL*8 zsolsw_aero(kdlon,9), zsolsw0_aero(kdlon,9) + + ! initialisation + tauaero(:,:,:,:)=0. + pizaero(:,:,:,:)=0. + cgaero(:,:,:,:)=0. + + ! + !------------------------------------------- + nb_gr = KLON / kdlon + IF (nb_gr*kdlon .NE. KLON) THEN + PRINT*, "kdlon mauvais:", KLON, kdlon, nb_gr + CALL abort + ENDIF + IF (kflev .NE. KLEV) THEN + PRINT*, "kflev differe de KLEV, kflev, KLEV" + CALL abort + ENDIF + !------------------------------------------- + DO k = 1, KLEV + DO i = 1, KLON + heat(i,k)=0. + cool(i,k)=0. + heat0(i,k)=0. + cool0(i,k)=0. + ENDDO + ENDDO + ! + zdist = dist + ! + PSCT = solaire/zdist/zdist + DO j = 1, nb_gr + iof = kdlon*(j-1) + DO i = 1, kdlon + zfract(i) = fract(iof+i) + zrmu0(i) = rmu0(iof+i) + PALBD(i,1) = albedo(iof+i) + PALBD(i,2) = alblw(iof+i) + PALBP(i,1) = albedo(iof+i) + PALBP(i,2) = alblw(iof+i) + PEMIS(i) = 1.0 + PVIEW(i) = 1.66 + PPSOL(i) = paprs(iof+i,1) + zx_alpha1 = (paprs(iof+i,1)-pplay(iof+i,2))/(pplay(iof+i,1)-pplay(iof+i,2)) + zx_alpha2 = 1.0 - zx_alpha1 + PTL(i,1) = t(iof+i,1) * zx_alpha1 + t(iof+i,2) * zx_alpha2 + PTL(i,KLEV+1) = t(iof+i,KLEV) + PDT0(i) = tsol(iof+i) - PTL(i,1) + ENDDO + DO k = 2, kflev + DO i = 1, kdlon + PTL(i,k) = (t(iof+i,k)+t(iof+i,k-1))*0.5 + ENDDO + ENDDO + DO k = 1, kflev + DO i = 1, kdlon + PDP(i,k) = paprs(iof+i,k)-paprs(iof+i,k+1) + PTAVE(i,k) = t(iof+i,k) + PWV(i,k) = MAX (q(iof+i,k), 1.0e-12) + PQS(i,k) = PWV(i,k) + ! wo: cm.atm (epaisseur en cm dans la situation standard) + ! POZON: kg/kg + POZON(i,k) = MAX(wo(iof+i,k),1.0e-12)*RG/46.6968 & + /(paprs(iof+i,k)-paprs(iof+i,k+1))& + *(paprs(iof+i,1)/101325.0) + PCLDLD(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k) + PCLDLU(i,k) = cldfra(iof+i,k)*cldemi(iof+i,k) + PCLDSW(i,k) = cldfra(iof+i,k) + PTAU(i,1,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! 1e-12 serait instable + PTAU(i,2,k) = MAX(cldtaupi(iof+i,k), 1.0e-05)! pour 32-bit machines + POMEGA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAU(i,1,k)) + POMEGA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAU(i,2,k)) + PCG(i,1,k) = 0.865 + PCG(i,2,k) = 0.910 + !- + ! Introduced for aerosol indirect forcings. + ! The following values use the cloud optical thickness calculated from + ! present-day aerosol concentrations whereas the quantities without the + ! "A" at the end are for pre-industial (natural-only) aerosol concentrations + ! + PTAUA(i,1,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! 1e-12 serait instable + PTAUA(i,2,k) = MAX(cldtaupd(iof+i,k), 1.0e-05)! pour 32-bit machines + POMEGAA(i,1,k) = 0.9999 - 5.0e-04 * EXP(-0.5 * PTAUA(i,1,k)) + POMEGAA(i,2,k) = 0.9988 - 2.5e-03 * EXP(-0.05 * PTAUA(i,2,k)) + ENDDO + ENDDO + ! + DO k = 1, kflev+1 + DO i = 1, kdlon + PPMB(i,k) = paprs(iof+i,k)/100.0 + ENDDO + ENDDO + ! + DO kk = 1, 5 + DO k = 1, kflev + DO i = 1, kdlon + PAER(i,k,kk) = 1.0E-15 + ENDDO + ENDDO + ENDDO + DO k = 1, kflev + DO i = 1, kdlon + tauaero(i,k,:,1)=tau_aero(iof+i,k,:,1) + pizaero(i,k,:,1)=piz_aero(iof+i,k,:,1) + cgaero(i,k,:,1) =cg_aero(iof+i,k,:,1) + tauaero(i,k,:,2)=tau_aero(iof+i,k,:,2) + pizaero(i,k,:,2)=piz_aero(iof+i,k,:,2) + cgaero(i,k,:,2) =cg_aero(iof+i,k,:,2) + ENDDO + ENDDO + ! + !====================================================================== + CALL LW_LMDAR4(& + PPMB, PDP,& + PPSOL,PDT0,PEMIS,& + PTL, PTAVE, PWV, POZON, PAER,& + PCLDLD,PCLDLU,& + PVIEW,& + zcool, zcool0,& + ztoplw,zsollw,ztoplw0,zsollw0,& + zsollwdown,& + ZFLUP, ZFLDN, ZFLUP0,ZFLDN0) + + + IF (.NOT. new_aod) THEN + ! use old version + CALL SW_LMDAR4(PSCT, zrmu0, zfract,& + PPMB, PDP, & + PPSOL, PALBD, PALBP,& + PTAVE, PWV, PQS, POZON, PAER,& + PCLDSW, PTAU, POMEGA, PCG,& + zheat, zheat0,& + zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,& + ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,& + tau_aero(:,:,5,:), piz_aero(:,:,5,:), cg_aero(:,:,5,:),& + PTAUA, POMEGAA,& + ztopswadaero,zsolswadaero,& + ztopswaiaero,zsolswaiaero,& + ok_ade, ok_aie) + ELSE + + CALL SW_AERO(PSCT, zrmu0, zfract,& + PPMB, PDP,& + PPSOL, PALBD, PALBP,& + PTAVE, PWV, PQS, POZON, PAER,& + PCLDSW, PTAU, POMEGA, PCG,& + zheat, zheat0,& + zalbpla,ztopsw,zsolsw,ztopsw0,zsolsw0,& + ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,& + tauaero, pizaero, cgaero, & + PTAUA, POMEGAA,& + ztopswadaero,zsolswadaero,& + ztopswad0aero,zsolswad0aero,& + ztopswaiaero,zsolswaiaero, & + ztopsw_aero,ztopsw0_aero,& + zsolsw_aero,zsolsw0_aero,& + ok_ade, ok_aie) + + ENDIF + !====================================================================== + DO i = 1, kdlon + radsol(iof+i) = zsolsw(i) + zsollw(i) + topsw(iof+i) = ztopsw(i) + toplw(iof+i) = ztoplw(i) + solsw(iof+i) = zsolsw(i) + sollw(iof+i) = zsollw(i) + sollwdown(iof+i) = zsollwdown(i) + DO k = 1, kflev+1 + lwdn0 ( iof+i,k) = ZFLDN0 ( i,k) + lwdn ( iof+i,k) = ZFLDN ( i,k) + lwup0 ( iof+i,k) = ZFLUP0 ( i,k) + lwup ( iof+i,k) = ZFLUP ( i,k) + ENDDO + topsw0(iof+i) = ztopsw0(i) + toplw0(iof+i) = ztoplw0(i) + solsw0(iof+i) = zsolsw0(i) + sollw0(iof+i) = zsollw0(i) + albpla(iof+i) = zalbpla(i) + + DO k = 1, kflev+1 + swdn0 ( iof+i,k) = ZFSDN0 ( i,k) + swdn ( iof+i,k) = ZFSDN ( i,k) + swup0 ( iof+i,k) = ZFSUP0 ( i,k) + swup ( iof+i,k) = ZFSUP ( i,k) + ENDDO + ENDDO + !-transform the aerosol forcings, if they have + ! to be calculated + IF (ok_ade) THEN + DO i = 1, kdlon + topswad_aero(iof+i) = ztopswadaero(i) + topswad0_aero(iof+i) = ztopswad0aero(i) + solswad_aero(iof+i) = zsolswadaero(i) + solswad0_aero(iof+i) = zsolswad0aero(i) + topsw_aero(iof+i,:) = ztopsw_aero(iof+i,:) + topsw0_aero(iof+i,:) = ztopsw0_aero(iof+i,:) + solsw_aero(iof+i,:) = zsolsw_aero(iof+i,:) + solsw0_aero(iof+i,:) = zsolsw0_aero(iof+i,:) + + ENDDO + ELSE + DO i = 1, kdlon + topswad_aero(iof+i) = 0.0 + solswad_aero(iof+i) = 0.0 + topswad0_aero(iof+i) = 0.0 + solswad0_aero(iof+i) = 0.0 + topsw_aero(iof+i,:) = 0. + topsw0_aero(iof+i,:) =0. + solsw_aero(iof+i,:) = 0. + solsw0_aero(iof+i,:) = 0. + ENDDO + ENDIF + IF (ok_aie) THEN + DO i = 1, kdlon + topswai_aero(iof+i) = ztopswaiaero(i) + solswai_aero(iof+i) = zsolswaiaero(i) + ENDDO + ELSE + DO i = 1, kdlon + topswai_aero(iof+i) = 0.0 + solswai_aero(iof+i) = 0.0 + ENDDO + ENDIF + DO k = 1, kflev + DO i = 1, kdlon + ! scale factor to take into account the difference between + ! dry air and watter vapour scpecifi! heat capacity + zznormcp=1.0+RVTMP2*PWV(i,k) + heat(iof+i,k) = zheat(i,k)/zznormcp + cool(iof+i,k) = zcool(i,k)/zznormcp + heat0(iof+i,k) = zheat0(i,k)/zznormcp + cool0(iof+i,k) = zcool0(i,k)/zznormcp + ENDDO + ENDDO + ! + ENDDO + + +ENDSUBROUTINE radlwsw_aero + + Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90 (revision 1150) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90 (revision 1150) @@ -0,0 +1,701 @@ +! +! $Header$ +! +SUBROUTINE readaerosol (id_aero, r_day, first, 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. +! +! +! 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 in data: + ! 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)) + + ENDIF + ! 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 ! 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 + + + + + + + + + + + + + + + + Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/readsulfate.F =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/readsulfate.F (revision 1149) +++ (revision ) @@ -1,660 +1,0 @@ -! -! $Header$ -! - SUBROUTINE readsulfate (r_day, first, sulfate_p) - USE dimphy, ONLY : klev - USE mod_grid_phy_lmdz, klon=>klon_glo - USE mod_phys_lmdz_para - IMPLICIT none - -c Content: -c -------- -c This routine reads in monthly mean values of sulfate aerosols and -c returns a linearly interpolated dayly-mean field. -c -c -c Author: -c ------- -c Johannes Quaas (quaas@lmd.jussieu.fr) -c 26/04/01 -c -c Modifications: -c -------------- -c 21/06/01: Make integrations of more than one year possible ;-) -c ATTENTION!! runs are supposed to start with Jan, 1. 1930 -c (rday=1) -c -c 27/06/01: Correction: The model always has 360 days per year! -c 27/06/01: SO4 concentration rather than mixing ratio -c 27/06/01: 10yr-mean-values to interpolate -c 20/08/01: Correct the error through integer-values in interpolations -c 21/08/01: Introduce flag to read in just one decade -c -c Include-files: -c -------------- -#include "YOMCST.h" -#include "chem.h" -#include "dimensions.h" -#include "temps.h" -#include "clesphys.h" -#include "iniprint.h" -c -c Input: -c ------ - REAL r_day ! Day of integration - LOGICAL first ! First timestep - ! (and therefore initialization necessary) -c -c Output: -c ------- - REAL sulfate_p(klon_omp,klev) - REAL sulfate (klon, klev) ! Mass of sulfate (monthly mean data, - ! from file) [ug SO4/m3] -c -c Local Variables: -c ---------------- - INTEGER i, ig, k, it - INTEGER j, iday, ny, iyr, iyr1, iyr2 - parameter (ny=jjm+1) - -CJLD INTEGER idec1, idec2 ! The two decadal data read ini - CHARACTER*4 cyear - - INTEGER im, day1, day2, im2 - REAL so4_1(iim, jjm+1, klev, 12) - REAL so4_2(iim, jjm+1, klev, 12) ! The sulfate distributions - - REAL, allocatable,save :: so4(:, :, :) ! SO4 in right dimension - REAL, allocatable,save :: so4_out(:, :) -c$OMP THREADPRIVATE(so4,so4_out) - - LOGICAL lnewday - LOGICAL lonlyone - PARAMETER (lonlyone=.FALSE.) - logical,save :: first2=.true. -c$OMP THREADPRIVATE(first2) - -c$OMP MASTER - if (first2) then - - allocate( so4(klon, klev, 12) ) - allocate( so4_out(klon, klev)) - first2=.false. - - endif - - if (is_mpi_root) then - - IF (aer_type /= 'actuel ' .AND. aer_type /= 'preind ' .AND. & - & aer_type /= 'scenario') THEN - WRITE(lunout,*)' *** Warning ***' - WRITE(lunout,*)'Option aer_type pour les aerosols = ', & - & aer_type - WRITE(lunout,*)'Cas non prevu, force a preind' - aer_type = 'preind ' - ENDIF - - 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 - so4(i,k,it)=0. - ENDDO - ENDDO - ENDDO - - - - IF (aer_type == 'actuel ') then - cyear='1980' - CALL getso4fromfile(cyear, so4_1) - ELSE IF (aer_type == 'preind ') THEN - cyear='.nat' - CALL getso4fromfile(cyear, so4_1) - ELSE - IF (iyr .lt. 1850) THEN - cyear='.nat' - WRITE(*,*) 'getso4 iyr=', iyr,' ',cyear - CALL getso4fromfile(cyear, so4_1) - ELSE IF (iyr .ge. 2100) THEN - cyear='2100' - WRITE(*,*) 'getso4 iyr=', iyr,' ',cyear - CALL getso4fromfile(cyear, so4_1) - ELSE - - ! Read in data: - ! 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 - ENDIF - WRITE(*,*) 'getso4 iyr=', iyr,' ',cyear - CALL getso4fromfile(cyear, so4_1) - - - ! If to read two decades: - IF (.NOT.lonlyone) THEN - - ! b) from the next following one - WRITE(cyear,'(I4)') iyr2 - WRITE(*,*) 'getso4 iyr=', iyr,' ',cyear - CALL getso4fromfile(cyear, so4_2) - - - ! Interpolate linarily to the actual year: - DO it=1,12 - DO k=1,klev - DO j=1,jjm - DO i=1,iim - so4_1(i,j,k,it)=so4_1(i,j,k,it) - . - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1) - . * (so4_1(i,j,k,it) - so4_2(i,j,k,it)) - ENDDO - ENDDO - ENDDO - ENDDO - - - ENDIF !lonlyone - 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 - so4(1,k,it)=so4(1,k,it)+so4_1(i,jjm+1,klev+1-k,it) - ! South pole - so4(klon,k,it)=so4(klon,k,it)+so4_1(i,1,klev+1-k,it) - ENDDO - so4(1,k,it)=so4(1,k,it)/FLOAT(iim) - so4(klon,k,it)=so4(klon,k,it)/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) write (*,*) 'shit' - so4(ig,k,it) = so4_1(i,jjm+1+1-j,klev+1-k,it) - ENDDO - ENDDO - IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)' - 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 - sulfate(i,k) = so4(i,k,im2) - . - FLOAT(iday-day2)/FLOAT(day1-day2) - . * (so4(i,k,im2) - so4(i,k,im)) - IF (sulfate(i,k).LT.0.) THEN - IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2 - IF (so4(i,k,im2) - so4(i,k,im).LT.0.) - . write(*,*) 'so4(i,k,im2) - so4(i,k,im)', - . so4(i,k,im2) - so4(i,k,im) - IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2 - stop 'sulfate' - 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 - sulfate(i,k) = so4(i,k,im2) - . - FLOAT(iday-day2)/FLOAT(day1-day2) - . * (so4(i,k,im2) - so4(i,k,im)) - IF (sulfate(i,k).LT.0.) THEN - IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2 - IF (so4(i,k,im2) - so4(i,k,im).LT.0.) - . write(*,*) 'so4(i,k,im2) - so4(i,k,im)', - . so4(i,k,im2) - so4(i,k,im) - IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2 - stop 'sulfate' - endif - ENDDO - ENDDO - ENDIF - - -CJLD ! The sulfate concentration [molec cm-3] is read in. -CJLD ! Convert it into mass [ug SO4/m3] -CJLD ! masse_so4 in [g/mol], n_avogadro in [molec/mol] - ! The sulfate mass [ug SO4/m3] is read in. - DO k=1,klev - DO i=1,klon -CJLD sulfate(i,k) = sulfate(i,k)*masse_so4 -CJLD . /n_avogadro*1.e+12 - so4_out(i,k) = sulfate(i,k) - IF (so4_out(i,k).LT.0) - . stop 'WAS SOLL DER SCHEISS ? ' - ENDDO - ENDDO - ELSE ! if no new day, use old data: - DO k=1,klev - DO i=1,klon - sulfate(i,k) = so4_out(i,k) - IF (so4_out(i,k).LT.0) - . stop 'WAS SOLL DER SCHEISS ? ' - ENDDO - ENDDO - - - ENDIF ! Did I have to do anything (was it a new day?) - - endif ! phy_rank==0 - -c$OMP END MASTER - call Scatter(sulfate,sulfate_p) - - RETURN - END - - - - - -c----------------------------------------------------------------------------- -c Read in /calculate pre-industrial values of sulfate -c----------------------------------------------------------------------------- - - SUBROUTINE readsulfate_preind (r_day, first, pi_sulfate_p) - USE dimphy, ONLY : klev - USE mod_grid_phy_lmdz, klon=>klon_glo - USE mod_phys_lmdz_para - IMPLICIT none - -c Content: -c -------- -c This routine reads in monthly mean values of sulfate aerosols and -c returns a linearly interpolated dayly-mean field. -c -c It does so for the preindustriel values of the sulfate, to a large part -c analogous to the routine readsulfate above. -c -c Only Pb: Variables must be saved and don t have to be overwritten! -c -c Author: -c ------- -c Johannes Quaas (quaas@lmd.jussieu.fr) -c 26/06/01 -c -c Modifications: -c -------------- -c see above -c -c Include-files: -c -------------- -#include "YOMCST.h" -#include "chem.h" -#include "dimensions.h" -#include "temps.h" -c -c Input: -c ------ - REAL r_day ! Day of integration - LOGICAL first ! First timestep - ! (and therefore initialization necessary) -c -c Output: -c ------- - REAL pi_sulfate_p (klon_omp, klev) - - REAL pi_sulfate (klon, klev) ! Number conc. sulfate (monthly mean data, - ! from fil -c -c Local Variables: -c ---------------- - INTEGER i, ig, k, it - INTEGER j, iday, ny, iyr - parameter (ny=jjm+1) - - INTEGER im, day1, day2, im2 - REAL pi_so4_1(iim, jjm+1, klev, 12) - - REAL, allocatable,save :: pi_so4(:, :, :) ! SO4 in right dimension - REAL, allocatable,save :: pi_so4_out(:, :) -c$OMP THREADPRIVATE(pi_so4,pi_so4_out) - - CHARACTER*4 cyear - LOGICAL lnewday - logical,save :: first2=.true. -c$OMP THREADPRIVATE(first2) - -c$OMP MASTER - if (first2) then - - allocate( pi_so4(klon, klev, 12) ) - allocate( pi_so4_out(klon, klev)) - 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 getso4fromfile(cyear,pi_so4_1) - - ! 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_so4(i,k,it)=0. - ENDDO - ENDDO - ENDDO - - write (*,*) '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_so4(1,k,it)=pi_so4(1,k,it)+pi_so4_1(i,jjm+1,klev+1-k,it) - ! South pole - pi_so4(klon,k,it)=pi_so4(klon,k,it)+pi_so4_1(i,1,klev+1-k,it) - ENDDO - pi_so4(1,k,it)=pi_so4(1,k,it)/FLOAT(iim) - pi_so4(klon,k,it)=pi_so4(klon,k,it)/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) write (*,*) 'shit' - pi_so4(ig,k,it) = pi_so4_1(i,jjm+1+1-j,klev+1-k,it) - ENDDO - ENDDO - IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)' - 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_sulfate(i,k) = pi_so4(i,k,im2) - . - FLOAT(iday-day2)/FLOAT(day1-day2) - . * (pi_so4(i,k,im2) - pi_so4(i,k,im)) - IF (pi_sulfate(i,k).LT.0.) THEN - IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2 - IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) - . write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', - . pi_so4(i,k,im2) - pi_so4(i,k,im) - IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2 - stop 'pi_sulfate' - 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_sulfate(i,k) = pi_so4(i,k,im2) - . - FLOAT(iday-day2)/FLOAT(day1-day2) - . * (pi_so4(i,k,im2) - pi_so4(i,k,im)) - IF (pi_sulfate(i,k).LT.0.) THEN - IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2 - IF (pi_so4(i,k,im2) - pi_so4(i,k,im).LT.0.) - . write(*,*) 'pi_so4(i,k,im2) - pi_so4(i,k,im)', - . pi_so4(i,k,im2) - pi_so4(i,k,im) - IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2 - stop 'pi_sulfate' - endif - ENDDO - ENDDO - ENDIF - - -CJLD ! The sulfate concentration [molec cm-3] is read in. -CJLD ! Convert it into mass [ug SO4/m3] -CJLD ! masse_so4 in [g/mol], n_avogadro in [molec/mol] - DO k=1,klev - DO i=1,klon -CJLD pi_sulfate(i,k) = pi_sulfate(i,k)*masse_so4 -CJLD . /n_avogadro*1.e+12 - pi_so4_out(i,k) = pi_sulfate(i,k) - ENDDO - ENDDO - - ELSE ! If no new day, use old data: - DO k=1,klev - DO i=1,klon - pi_sulfate(i,k) = pi_so4_out(i,k) - ENDDO - ENDDO - - - ENDIF ! Was this the beginning of a new day? - - endif ! is_mpi_root==0 - -c$OMP END MASTER - call Scatter(pi_sulfate,pi_sulfate_p) - - RETURN - END - - - - - - - - - - -c----------------------------------------------------------------------------- -c Routine for reading SO4 data from files -c----------------------------------------------------------------------------- - - - SUBROUTINE getso4fromfile (cyr, so4) - USE dimphy -#include "netcdf.inc" -#include "dimensions.h" - CHARACTER*15 fname - CHARACTER*4 cyr - - CHARACTER*6 cvar - INTEGER START(3), COUNT(3) - INTEGER STATUS, NCID, VARID - INTEGER imth, i, j, k, ny - PARAMETER (ny=jjm+1) - - - REAL so4mth(iim, ny, klev) - REAL so4(iim, ny, klev, 12) - - - fname = 'so4.run'//cyr//'.cdf' - - write (*,*) 'reading ', fname - STATUS = NF_OPEN (fname, NF_NOWRITE, NCID) - IF (STATUS .NE. NF_NOERR) write (*,*) 'err in open ',status - - DO imth=1, 12 - IF (imth.eq.1) THEN - cvar='SO4JAN' - ELSEIF (imth.eq.2) THEN - cvar='SO4FEB' - ELSEIF (imth.eq.3) THEN - cvar='SO4MAR' - ELSEIF (imth.eq.4) THEN - cvar='SO4APR' - ELSEIF (imth.eq.5) THEN - cvar='SO4MAY' - ELSEIF (imth.eq.6) THEN - cvar='SO4JUN' - ELSEIF (imth.eq.7) THEN - cvar='SO4JUL' - ELSEIF (imth.eq.8) THEN - cvar='SO4AUG' - ELSEIF (imth.eq.9) THEN - cvar='SO4SEP' - ELSEIF (imth.eq.10) THEN - cvar='SO4OCT' - ELSEIF (imth.eq.11) THEN - cvar='SO4NOV' - ELSEIF (imth.eq.12) THEN - cvar='SO4DEC' - ENDIF - start(1)=1 - start(2)=1 - start(3)=1 - count(1)=iim - count(2)=ny - count(3)=klev -c write(*,*) 'here i am' - STATUS = NF_INQ_VARID (NCID, cvar, VARID) - write (*,*) ncid,imth,cvar, varid - - IF (STATUS .NE. NF_NOERR) write (*,*) 'err in read ',status - -#ifdef NC_DOUBLE - status = NF_GET_VARA_DOUBLE(NCID, VARID, START, COUNT, so4mth) -#else - status = NF_GET_VARA_REAL(NCID, VARID, START, COUNT, so4mth) -#endif - IF (STATUS .NE. NF_NOERR) write (*,*) 'err in read data',status - - DO k=1,klev - DO j=1,jjm+1 - DO i=1,iim - IF (so4mth(i,j,k).LT.0.) then - write(*,*) 'this is shit' - write(*,*) 'so4(',i,j,k,') =',so4mth(i,j,k) - endif - so4(i,j,k,imth)=so4mth(i,j,k) - ENDDO - ENDDO - ENDDO - ENDDO - - STATUS = NF_CLOSE(NCID) - IF (STATUS .NE. NF_NOERR) write (*,*) 'err in closing file',status - - - END ! subroutine getso4fromfile - - - - - - - - - - - - - - - - Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/sw_aero.F90 =================================================================== --- LMDZ4/branches/LMDZ4-dev/libf/phylmd/sw_aero.F90 (revision 1150) +++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/sw_aero.F90 (revision 1150) @@ -0,0 +1,750 @@ +SUBROUTINE SW_AERO(PSCT, PRMU0, PFRAC, & + PPMB, PDP, & + PPSOL, PALBD, PALBP,& + PTAVE, PWV, PQS, POZON, PAER,& + PCLDSW, PTAU, POMEGA, PCG,& + PHEAT, PHEAT0,& + PALBPLA,PTOPSW,PSOLSW,PTOPSW0,PSOLSW0,& + ZFSUP,ZFSDN,ZFSUP0,ZFSDN0,& + tauaero, pizaero, cgaero,& + PTAUA, POMEGAA,& + PTOPSWADAERO,PSOLSWADAERO,& + PTOPSWAD0AERO,PSOLSWAD0AERO,& + PTOPSWAIAERO,PSOLSWAIAERO,& + PTOPSWAERO,PTOPSW0AERO,& + PSOLSWAERO,PSOLSW0AERO,& + ok_ade, ok_aie ) + + USE dimphy + IMPLICIT NONE + +#include "YOMCST.h" + ! + ! ------------------------------------------------------------------ + ! + ! PURPOSE. + ! -------- + ! + ! THIS ROUTINE COMPUTES THE SHORTWAVE RADIATION FLUXES IN TWO + ! SPECTRAL INTERVALS FOLLOWING FOUQUART AND BONNEL (1980). + ! + ! METHOD. + ! ------- + ! + ! 1. COMPUTES ABSORBER AMOUNTS (SWU) + ! 2. COMPUTES FLUXES IN 1ST SPECTRAL INTERVAL (SW1S) + ! 3. COMPUTES FLUXES IN 2ND SPECTRAL INTERVAL (SW2S) + ! + ! REFERENCE. + ! ---------- + ! + ! SEE RADIATION'S PART OF THE ECMWF RESEARCH DEPARTMENT + ! DOCUMENTATION, AND FOUQUART AND BONNEL (1980) + ! + ! AUTHOR. + ! ------- + ! JEAN-JACQUES MORCRETTE *ECMWF* + ! + ! MODIFICATIONS. + ! -------------- + ! ORIGINAL : 89-07-14 + ! 95-01-01 J.-J. MORCRETTE Direct/Diffuse Albedo + ! 03-11-27 J. QUAAS Introduce aerosol forcings (based on BOUCHER) + ! 09-04 A. COZIC - C.DEANDREIS Indroduce NAT/BC/POM/DUST/SS aerosol forcing + ! ------------------------------------------------------------------ + ! + !* ARGUMENTS: + ! + REAL*8 PSCT ! constante solaire (valeur conseillee: 1370) + + REAL*8 PPSOL(KDLON) ! SURFACE PRESSURE (PA) + REAL*8 PDP(KDLON,KFLEV) ! LAYER THICKNESS (PA) + REAL*8 PPMB(KDLON,KFLEV+1) ! HALF-LEVEL PRESSURE (MB) + + REAL*8 PRMU0(KDLON) ! COSINE OF ZENITHAL ANGLE + REAL*8 PFRAC(KDLON) ! fraction de la journee + + REAL*8 PTAVE(KDLON,KFLEV) ! LAYER TEMPERATURE (K) + REAL*8 PWV(KDLON,KFLEV) ! SPECIFI! HUMIDITY (KG/KG) + REAL*8 PQS(KDLON,KFLEV) ! SATURATED WATER VAPOUR (KG/KG) + REAL*8 POZON(KDLON,KFLEV) ! OZONE CONCENTRATION (KG/KG) + REAL*8 PAER(KDLON,KFLEV,5) ! AEROSOLS' OPTICAL THICKNESS + + REAL*8 PALBD(KDLON,2) ! albedo du sol (lumiere diffuse) + REAL*8 PALBP(KDLON,2) ! albedo du sol (lumiere parallele) + + REAL*8 PCLDSW(KDLON,KFLEV) ! CLOUD FRACTION + REAL*8 PTAU(KDLON,2,KFLEV) ! CLOUD OPTICAL THICKNESS + REAL*8 PCG(KDLON,2,KFLEV) ! ASYMETRY FACTOR + REAL*8 POMEGA(KDLON,2,KFLEV) ! SINGLE SCATTERING ALBEDO + + REAL*8 PHEAT(KDLON,KFLEV) ! SHORTWAVE HEATING (K/DAY) + REAL*8 PHEAT0(KDLON,KFLEV)! SHORTWAVE HEATING (K/DAY) clear-sky + REAL*8 PALBPLA(KDLON) ! PLANETARY ALBEDO + REAL*8 PTOPSW(KDLON) ! SHORTWAVE FLUX AT T.O.A. + REAL*8 PSOLSW(KDLON) ! SHORTWAVE FLUX AT SURFACE + REAL*8 PTOPSW0(KDLON) ! SHORTWAVE FLUX AT T.O.A. (CLEAR-SKY) + REAL*8 PSOLSW0(KDLON) ! SHORTWAVE FLUX AT SURFACE (CLEAR-SKY) + ! + !* LOCAL VARIABLES: + ! + REAL*8 ZOZ(KDLON,KFLEV) + REAL*8 ZAKI(KDLON,2) + REAL*8 ZCLD(KDLON,KFLEV) + REAL*8 ZCLEAR(KDLON) + REAL*8 ZDSIG(KDLON,KFLEV) + REAL*8 ZFACT(KDLON) + REAL*8 ZFD(KDLON,KFLEV+1) + REAL*8 ZFDOWN(KDLON,KFLEV+1) + REAL*8 ZFU(KDLON,KFLEV+1) + REAL*8 ZFUP(KDLON,KFLEV+1) + REAL*8 ZRMU(KDLON) + REAL*8 ZSEC(KDLON) + REAL*8 ZUD(KDLON,5,KFLEV+1) + REAL*8 ZCLDSW0(KDLON,KFLEV) + + REAL*8 ZFSUP(KDLON,KFLEV+1) + REAL*8 ZFSDN(KDLON,KFLEV+1) + REAL*8 ZFSUP0(KDLON,KFLEV+1) + REAL*8 ZFSDN0(KDLON,KFLEV+1) + + INTEGER inu, jl, jk, i, k, kpl1 + + INTEGER swpas ! Every swpas steps, sw is calculated + PARAMETER(swpas=1) + + INTEGER itapsw + LOGICAL appel1er + DATA itapsw /0/ + DATA appel1er /.TRUE./ + SAVE itapsw,appel1er + !$OMP THREADPRIVATE(appel1er) + !$OMP THREADPRIVATE(itapsw) + !jq-Introduced for aerosol forcings + REAL*8 flag_aer + LOGICAL ok_ade, ok_aie ! use aerosol forcings or not? + REAL*8 tauaero(kdlon,kflev,9,2) ! aerosol optical properties + REAL*8 pizaero(kdlon,kflev,9,2) ! (see aeropt.F) + REAL*8 cgaero(kdlon,kflev,9,2) ! -"- + REAL*8 PTAUA(KDLON,2,KFLEV) ! CLOUD OPTICAL THICKNESS (pre-industrial value) + REAL*8 POMEGAA(KDLON,2,KFLEV) ! SINGLE SCATTERING ALBEDO + REAL*8 PTOPSWADAERO(KDLON) ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR) + REAL*8 PSOLSWADAERO(KDLON) ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR) + REAL*8 PTOPSWAD0AERO(KDLON) ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL DIR) + REAL*8 PSOLSWAD0AERO(KDLON) ! SHORTWAVE FLUX AT SURFACE(+AEROSOL DIR) + REAL*8 PTOPSWAIAERO(KDLON) ! SHORTWAVE FLUX AT T.O.A.(+AEROSOL IND) + REAL*8 PSOLSWAIAERO(KDLON) ! SHORTWAVE FLUX AT SURFACE(+AEROSOL IND) + REAL*8 PTOPSWAERO(KDLON,9) + REAL*8 PTOPSW0AERO(KDLON,9) + REAL*8 PSOLSWAERO(KDLON,9) + REAL*8 PSOLSW0AERO(KDLON,9) + + !jq - Fluxes including aerosol effects + REAL*8,ALLOCATABLE,SAVE :: ZFSUPAD_AERO(:,:) + !$OMP THREADPRIVATE(ZFSUPAD_AERO) + REAL*8,ALLOCATABLE,SAVE :: ZFSDNAD_AERO(:,:) + !$OMP THREADPRIVATE(ZFSDNAD_AERO) + !jq - Fluxes including aerosol effects + REAL*8,ALLOCATABLE,SAVE :: ZFSUPAD0_AERO(:,:) + !$OMP THREADPRIVATE(ZFSUPAD0_AERO) + REAL*8,ALLOCATABLE,SAVE :: ZFSDNAD0_AERO(:,:) + !$OMP THREADPRIVATE(ZFSDNAD0_AERO) + REAL*8,ALLOCATABLE,SAVE :: ZFSUPAI_AERO(:,:) + !$OMP THREADPRIVATE(ZFSUPAI_AERO) + REAL*8,ALLOCATABLE,SAVE :: ZFSDNAI_AERO(:,:) + !$OMP THREADPRIVATE(ZFSDNAI_AERO) + REAL*8,ALLOCATABLE,SAVE :: ZFSUP_AERO(:,:,:) + !$OMP THREADPRIVATE(ZFSUP_AERO) + REAL*8,ALLOCATABLE,SAVE :: ZFSDN_AERO(:,:,:) + !$OMP THREADPRIVATE(ZFSDN_AERO) + REAL*8,ALLOCATABLE,SAVE :: ZFSUP0_AERO(:,:,:) + !$OMP THREADPRIVATE(ZFSUP0_AERO) + REAL*8,ALLOCATABLE,SAVE :: ZFSDN0_AERO(:,:,:) + !$OMP THREADPRIVATE(ZFSDN0_AERO) + + LOGICAL initialized + !rv + SAVE flag_aer + !$OMP THREADPRIVATE(flag_aer) + DATA initialized/.FALSE./ + SAVE initialized + !$OMP THREADPRIVATE(initialized) + + + IF(.NOT.initialized) THEN + flag_aer=0. + initialized=.TRUE. + ALLOCATE(ZFSUPAD_AERO(KDLON,KFLEV+1)) + ALLOCATE(ZFSDNAD_AERO(KDLON,KFLEV+1)) + ALLOCATE(ZFSUPAD0_AERO(KDLON,KFLEV+1)) + ALLOCATE(ZFSDNAD0_AERO(KDLON,KFLEV+1)) + ALLOCATE(ZFSUPAI_AERO(KDLON,KFLEV+1)) + ALLOCATE(ZFSDNAI_AERO(KDLON,KFLEV+1)) + ALLOCATE(ZFSUP_AERO (KDLON,KFLEV+1,9)) + ALLOCATE(ZFSDN_AERO (KDLON,KFLEV+1,9)) + ALLOCATE(ZFSUP0_AERO(KDLON,KFLEV+1,9)) + ALLOCATE(ZFSDN0_AERO(KDLON,KFLEV+1,9)) + ZFSUPAD_AERO(:,:)=0. + ZFSDNAD_AERO(:,:)=0. + ZFSUPAD0_AERO(:,:)=0. + ZFSDNAD0_AERO(:,:)=0. + ZFSUPAI_AERO(:,:)=0. + ZFSDNAI_AERO(:,:)=0. + ZFSUP_AERO (:,:,:)=0. + ZFSDN_AERO (:,:,:)=0. + ZFSUP0_AERO(:,:,:)=0. + ZFSDN0_AERO(:,:,:)=0. + ENDIF + !rv + + + IF (appel1er) THEN + PRINT*, 'SW calling frequency : ', swpas + PRINT*, " In general, it should be 1" + appel1er = .FALSE. + ENDIF + ! ------------------------------------------------------------------ + IF (MOD(itapsw,swpas).EQ.0) THEN + + DO JK = 1 , KFLEV + DO JL = 1, KDLON + ZCLDSW0(JL,JK) = 0.0 + ZOZ(JL,JK) = POZON(JL,JK)*46.6968/RG & + *PDP(JL,JK)*(101325.0/PPSOL(JL)) + ENDDO + ENDDO + + + ! clear-sky: + flag_aer=0.0 + CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU,PAER, flag_aer, & + tauaero(:,:,1,:), pizaero(:,:,1,:), cgaero(:,:,1,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, ZCLDSW0,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer, & + tauaero(:,:,1,:), pizaero(:,:,1,:), cgaero(:,:,1,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, ZCLDSW0,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUP0(JL,JK) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN0(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ZFSUP0_AERO(JL,JK,1) = ZFSUP0(JL,JK) + ZFSDN0_AERO(JL,JK,1) = ZFSDN0(JL,JK) + ENDDO + ENDDO + + + ! cloudy-sky: + flag_aer=0.0 + CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU, PAER, flag_aer, & + tauaero(:,:,1,:), pizaero(:,:,1,:), cgaero(:,:,1,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer, & + tauaero(:,:,1,:), pizaero(:,:,1,:), cgaero(:,:,1,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUP(JL,JK) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ZFSUP_AERO(JL,JK,1) = ZFSUP(JL,JK) + ZFSDN_AERO(JL,JK,1) = ZFSDN(JL,JK) + ENDDO + ENDDO + + ZFSUP0_AERO(:,:,2) = ZFSUP0_AERO(:,:,1) + ZFSDN0_AERO(:,:,2) = ZFSDN0_AERO(:,:,1) + ZFSUP_AERO(:,:,2) = ZFSUP_AERO(:,:,1) + ZFSDN_AERO(:,:,2) = ZFSDN_AERO(:,:,1) + + + IF (ok_ade) THEN + + ! clear sky (Anne Cozic 03/07/2007) + ! CAS AER (2) + flag_aer=1.0 + CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUPAD0_AERO(JL,JK) = ZFSUP0(JL,JK) + ZFSDNAD0_AERO(JL,JK) = ZFSDN0(JL,JK) + ZFSUP0(JL,JK) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN0(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ZFSUP0_AERO(JL,JK,2) = ZFSUP0(JL,JK) + ZFSDN0_AERO(JL,JK,2) = ZFSDN0(JL,JK) + ENDDO + ENDDO + + ! cloudy-sky + aerosol dir OB + ! ACo AER + flag_aer=1.0 + CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUPAD_AERO(JL,JK) = ZFSUP(JL,JK) + ZFSDNAD_AERO(JL,JK) = ZFSDN(JL,JK) + ZFSUP(JL,JK) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ZFSUP_AERO(JL,JK,2) = ZFSUP(JL,JK) + ZFSDN_AERO(JL,JK,2) = ZFSDN(JL,JK) + ENDDO + ENDDO + + !CAS NAT + ! clear sky + flag_aer=1.0 + CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUP0_AERO(JL,JK,3) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN0_AERO(JL,JK,3) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ENDDO + ENDDO + + ! cloudy-sky + ! ACo NAT + flag_aer=1.0 + CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,3,:), pizaero(:,:,3,:), cgaero(:,:,3,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUP_AERO(JL,JK,3) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN_AERO(JL,JK,3) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ENDDO + ENDDO + + ! clear sky (Yves 01/12/2007) + ! CAS BC (4) + flag_aer=1.0 + CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,4,:), pizaero(:,:,4,:), cgaero(:,:,4,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,4,:), pizaero(:,:,4,:), cgaero(:,:,4,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUP0_AERO(JL,JK,4) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN0_AERO(JL,JK,4) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ENDDO + ENDDO + + ! cloudy-sky + aerosol dir OB + ! CAS BC (4) + flag_aer=1.0 + CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,4,:), pizaero(:,:,4,:), cgaero(:,:,4,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,4,:), pizaero(:,:,4,:), cgaero(:,:,4,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUP_AERO(JL,JK,4) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN_AERO(JL,JK,4) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ENDDO + ENDDO + + ! clear sky (Yves 13/12/2007) + ! CAS SO4 (5) + flag_aer=1.0 + CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUP0_AERO(JL,JK,5) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN0_AERO(JL,JK,5) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ENDDO + ENDDO + + ! cloudy-sky + aerosol dir OB + ! CAS SO4 (5) + flag_aer=1.0 + CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,5,:), pizaero(:,:,5,:), cgaero(:,:,5,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUP_AERO(JL,JK,5) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN_AERO(JL,JK,5) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ENDDO + ENDDO + + ! clear sky (Yves 13/12/2007) + ! CAS POM (6) + flag_aer=1.0 + CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,6,:), pizaero(:,:,6,:), cgaero(:,:,6,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,6,:), pizaero(:,:,6,:), cgaero(:,:,6,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUP0_AERO(JL,JK,6) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN0_AERO(JL,JK,6) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ENDDO + ENDDO + + ! cloudy-sky + aerosol dir OB + ! CAS POM (6) + flag_aer=1.0 + CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,6,:), pizaero(:,:,6,:), cgaero(:,:,6,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,6,:), pizaero(:,:,6,:), cgaero(:,:,6,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUP_AERO(JL,JK,6) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN_AERO(JL,JK,6) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ENDDO + ENDDO + + ! clear sky (Yves 13/12/2007) + ! CAS DUST (7) + flag_aer=1.0 + CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,7,:), pizaero(:,:,7,:), cgaero(:,:,7,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,7,:), pizaero(:,:,7,:), cgaero(:,:,7,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUP0_AERO(JL,JK,7) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN0_AERO(JL,JK,7) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ENDDO + ENDDO + + ! cloudy-sky + aerosol dir OB + ! CAS DUST (7) + flag_aer=1.0 + CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,7,:), pizaero(:,:,7,:), cgaero(:,:,7,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,7,:), pizaero(:,:,7,:), cgaero(:,:,7,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUP_AERO(JL,JK,7) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN_AERO(JL,JK,7) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ENDDO + ENDDO + + ! clear sky (Yves 13/12/2007) + ! CAS Seasalt SS (8) + flag_aer=1.0 + CALL SWU_LMDAR4(PSCT,ZCLDSW0,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,8,:), pizaero(:,:,8,:), cgaero(:,:,8,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,8,:), pizaero(:,:,8,:), cgaero(:,:,8,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUP0_AERO(JL,JK,8) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN0_AERO(JL,JK,8) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ENDDO + ENDDO + + ! cloudy-sky + aerosol dir OB + ! CAS Seasalt SS (8) + flag_aer=1.0 + CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,8,:), pizaero(:,:,8,:), cgaero(:,:,8,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,8,:), pizaero(:,:,8,:), cgaero(:,:,8,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGA, ZOZ, ZRMU, ZSEC, PTAU, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUP_AERO(JL,JK,8) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN_AERO(JL,JK,8) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ENDDO + ENDDO + + ENDIF ! ok_ade + + + IF (ok_aie) THEN + !jq cloudy-sky + aerosol direct + aerosol indirect + flag_aer=1.0 + CALL SWU_LMDAR4(PSCT,PCLDSW,PPMB,PPSOL,& + PRMU0,PFRAC,PTAVE,PWV,& + ZAKI,ZCLD,ZCLEAR,ZDSIG,ZFACT,ZRMU,ZSEC,ZUD) + INU = 1 + CALL SW1S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),& + PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGAA, ZOZ, ZRMU, ZSEC, PTAUA, ZUD,& + ZFD, ZFU) + INU = 2 + CALL SW2S_LMDAR4(INU, PAER, flag_aer,& + tauaero(:,:,2,:), pizaero(:,:,2,:), cgaero(:,:,2,:),& + ZAKI, PALBD, PALBP, PCG, ZCLD, ZCLEAR, PCLDSW,& + ZDSIG, POMEGAA, ZOZ, ZRMU, ZSEC, PTAUA, ZUD,& + PWV, PQS,& + ZFDOWN, ZFUP) + DO JK = 1 , KFLEV+1 + DO JL = 1, KDLON + ZFSUPAI_AERO(JL,JK) = ZFSUP(JL,JK) + ZFSDNAI_AERO(JL,JK) = ZFSDN(JL,JK) + ZFSUP(JL,JK) = (ZFUP(JL,JK) + ZFU(JL,JK)) * ZFACT(JL) + ZFSDN(JL,JK) = (ZFDOWN(JL,JK) + ZFD(JL,JK)) * ZFACT(JL) + ZFSUP_AERO(JL,JK,2) = ZFSUP(JL,JK) + ZFSDN_AERO(JL,JK,2) = ZFSDN(JL,JK) + ENDDO + ENDDO + ENDIF ! ok_aie + + itapsw = 0 + ENDIF + itapsw = itapsw + 1 + + DO k = 1, KFLEV + kpl1 = k+1 + DO i = 1, KDLON + + PHEAT(i,k) = -(ZFSUP_AERO(i,kpl1,2)-ZFSUP_AERO(i,k,2)) & + -(ZFSDN_AERO(i,k,2)-ZFSDN_AERO(i,kpl1,2)) + PHEAT(i,k) = PHEAT(i,k) * RDAY*RG/RCPD / PDP(i,k) + PHEAT0(i,k) = -(ZFSUP0_AERO(i,kpl1,2)-ZFSUP0_AERO(i,k,2)) & + -(ZFSDN0_AERO(i,k,2)-ZFSDN0_AERO(i,kpl1,2)) + PHEAT0(i,k) = PHEAT0(i,k) * RDAY*RG/RCPD / PDP(i,k) + + ENDDO + ENDDO + DO i = 1, KDLON + PALBPLA(i) = ZFSUP(i,KFLEV+1)/(ZFSDN(i,KFLEV+1)+1.0e-20) + ! clear sky + PSOLSW0(i) = ZFSDN0(i,1) - ZFSUP0(i,1) + PTOPSW0(i) = ZFSDN0(i,KFLEV+1) - ZFSUP0(i,KFLEV+1) + + PSOLSW(i) = ZFSDN(i,1) - ZFSUP(i,1) + PTOPSW(i) = ZFSDN(i,KFLEV+1) - ZFSUP(i,KFLEV+1) + + PSOLSW0AERO(i,:) = ZFSDN0_AERO(i,1,:) - ZFSUP0_AERO(i,1,:) + PTOPSW0AERO(i,:) = & + ZFSDN0_AERO(i,KFLEV+1,:) - ZFSUP0_AERO(i,KFLEV+1,:) + + PSOLSWAERO(i,:) = ZFSDN_AERO(i,1,:) - ZFSUP_AERO(i,1,:) + PTOPSWAERO(i,:) = & + ZFSDN_AERO(i,KFLEV+1,:) - ZFSUP_AERO(i,KFLEV+1,:) + + PSOLSWADAERO(i) = ZFSDNAD_AERO(i,1) - ZFSUPAD_AERO(i,1) + PTOPSWADAERO(i) = & + ZFSDNAD_AERO(i,KFLEV+1) - ZFSUPAD_AERO(i,KFLEV+1) + + PSOLSWAD0AERO(i) = ZFSDNAD0_AERO(i,1) - ZFSUPAD0_AERO(i,1) + PTOPSWAD0AERO(i) = & + ZFSDNAD0_AERO(i,KFLEV+1) - ZFSUPAD0_AERO(i,KFLEV+1) + + PSOLSWAIAERO(i) = ZFSDNAI_AERO(i,1) - ZFSUPAI_AERO(i,1) + PTOPSWAIAERO(i) = & + ZFSDNAI_AERO(i,KFLEV+1) - ZFSUPAI_AERO(i,KFLEV+1) + + ENDDO +END SUBROUTINE SW_AERO +