Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/aeropt.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/aeropt.F	(revision 520)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/aeropt.F	(revision 520)
@@ -0,0 +1,130 @@
+      SUBROUTINE aeropt(pplay, paprs, t_seri, msulfate, RHcl,
+     .            tau_ae, piz_ae, cg_ae, ai        )
+c
+      IMPLICIT none
+c
+c
+c     
+#include "dimensions.h"
+#include "dimphy.h"
+#include "YOMCST.h"
+c
+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 
+c
+c Local
+c
+      INTEGER i, k, inu
+      INTEGER RH_num, nbre_RH
+      PARAMETER (nbre_RH=12)
+      REAL RH_tab(nbre_RH)
+      REAL RH_MAX, DELTA, rh 
+      PARAMETER (RH_MAX=95.)
+      DATA RH_tab/0.,10.,20.,30.,40.,50.,60.,70.,80.,85.,90.,95./
+      REAL zrho, zdz
+      REAL taue670(klon)       ! epaisseur optique aerosol absorption 550 nm
+      REAL taue865(klon)       ! epaisseur optique aerosol extinction 865 nm
+      REAL alpha_aer_sulfate(nbre_RH,5)   !--unit m2/g SO4
+      REAL alphasulfate      
+c
+c Proprietes optiques
+c
+      REAL alpha_aer(nbre_RH,2)   !--unit m2/g SO4
+      REAL cg_aer(nbre_RH,2)
+      DATA alpha_aer/.500130E+01,  .500130E+01,  .500130E+01,  
+     .               .500130E+01,  .500130E+01,  .616710E+01,  
+     .               .826850E+01,  .107687E+02,  .136976E+02,  
+     .               .162972E+02,  .211690E+02,  .354833E+02,  
+     .               .139460E+01,  .139460E+01,  .139460E+01,  
+     .               .139460E+01,  .139460E+01,  .173910E+01,  
+     .               .244380E+01,  .332320E+01,  .440120E+01,  
+     .               .539570E+01,  .734580E+01,  .136038E+02 / 
+      DATA cg_aer/.619800E+00,  .619800E+00,  .619800E+00,  
+     .            .619800E+00,  .619800E+00,  .662700E+00,  
+     .            .682100E+00,  .698500E+00,  .712500E+00,  
+     .            .721800E+00,  .734600E+00,  .755800E+00,  
+     .            .545600E+00,  .545600E+00,  .545600E+00,  
+     .            .545600E+00,  .545600E+00,  .583700E+00,  
+     .            .607100E+00,  .627700E+00,  .645800E+00,  
+     .            .658400E+00,  .676500E+00,  .708500E+00 / 
+      DATA alpha_aer_sulfate/
+     . 4.910,4.910,4.910,4.910,6.547,7.373,
+     . 8.373,9.788,12.167,14.256,17.924,28.433,
+     . 1.453,1.453,1.453,1.453,2.003,2.321,
+     . 2.711,3.282,4.287,5.210,6.914,12.305,
+     . 4.308,4.308,4.308,4.308,5.753,6.521,
+     . 7.449,8.772,11.014,12.999,16.518,26.772,
+     . 3.265,3.265,3.265,3.265,4.388,5.016,
+     . 5.775,6.868,8.745,10.429,13.457,22.538,
+     . 2.116,2.116,2.116,2.116,2.882,3.330,
+     . 3.876,4.670,6.059,7.327,9.650,16.883/
+c
+      DO i=1, klon
+         taue670(i)=0.0
+         taue865(i)=0.0
+      ENDDO
+c      
+      DO k=1, klev
+      DO i=1, klon
+         if (t_seri(i,k).eq.0) write (*,*) 'aeropt T ',i,k,t_seri(i,k)
+         if (pplay(i,k).eq.0) write (*,*) 'aeropt p ',i,k,pplay(i,k)         
+        zrho=pplay(i,k)/t_seri(i,k)/RD                  ! kg/m3
+        zdz=(paprs(i,k)-paprs(i,k+1))/zrho/RG           ! m
+        rh=MIN(RHcl(i,k)*100.,RH_MAX)
+        RH_num = INT( rh/10. + 1.)
+        IF (rh.LT.0.) STOP 'aeropt: RH < 0 not possible'
+        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))
+c                
+        inu=1
+        tau_ae(i,k,inu)=alpha_aer(RH_num,inu) +
+     .             DELTA*(alpha_aer(RH_num+1,inu)-alpha_aer(RH_num,inu))
+        tau_ae(i,k,inu)=tau_ae(i,k,inu)*msulfate(i,k)*zdz*1.e-6
+        piz_ae(i,k,inu)=1.0
+        cg_ae(i,k,inu)=cg_aer(RH_num,inu) +
+     .                 DELTA*(cg_aer(RH_num+1,inu)-cg_aer(RH_num,inu))
+c
+        inu=2
+        tau_ae(i,k,inu)=alpha_aer(RH_num,inu) +
+     .             DELTA*(alpha_aer(RH_num+1,inu)-alpha_aer(RH_num,inu))
+        tau_ae(i,k,inu)=tau_ae(i,k,inu)*msulfate(i,k)*zdz*1.e-6
+        piz_ae(i,k,inu)=1.0
+        cg_ae(i,k,inu)=cg_aer(RH_num,inu) +
+     .                 DELTA*(cg_aer(RH_num+1,inu)-cg_aer(RH_num,inu))
+cjq
+cjq for aerosol index
+c
+        alphasulfate=alpha_aer_sulfate(RH_num,4) +
+     .       DELTA*(alpha_aer_sulfate(RH_num+1,4)-
+     .       alpha_aer_sulfate(RH_num,4)) !--m2/g 
+c     
+        taue670(i)=taue670(i)+
+     .       alphasulfate*msulfate(i,k)*zdz*1.e-6
+c
+        alphasulfate=alpha_aer_sulfate(RH_num,5) +
+     .       DELTA*(alpha_aer_sulfate(RH_num+1,5)-
+     .       alpha_aer_sulfate(RH_num,5)) !--m2/g 
+c
+        taue865(i)=taue865(i)+
+     .         alphasulfate*msulfate(i,k)*zdz*1.e-6
+        
+      ENDDO
+      ENDDO
+c      
+      DO i=1, klon
+        ai(i)=(-log(MAX(taue670(i),0.0001)/
+     .                MAX(taue865(i),0.0001))/log(670./865.)) * 
+     .        taue865(i)
+      ENDDO     
+c
+      RETURN
+      END
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/chem.h
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/chem.h	(revision 520)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/chem.h	(revision 520)
@@ -0,0 +1,11 @@
+      INTEGER idms, iso2, iso4, ih2s, idmso, imsa, ih2o2
+      PARAMETER (idms=1, iso2=2, iso4=3)
+      PARAMETER (ih2s=4, idmso=5, imsa=6, ih2o2=7)
+
+      REAL n_avogadro, masse_s, masse_so4, rho_water, rho_ice
+      PARAMETER (n_avogadro=6.02E23)                  !--molec mol-1
+      PARAMETER (masse_s=32.0)                        !--g mol-1
+      PARAMETER (masse_so4=96.0)                      !--g mol-1
+      PARAMETER (rho_water=1000.0)                    !--kg m-3
+      PARAMETER (rho_ice=500.0)                       !--kg m-3
+
Index: LMDZ.3.3/branches/rel-LF/libf/phylmd/readsulfate.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/phylmd/readsulfate.F	(revision 520)
+++ LMDZ.3.3/branches/rel-LF/libf/phylmd/readsulfate.F	(revision 520)
@@ -0,0 +1,576 @@
+      SUBROUTINE readsulfate (r_day, first, sulfate)
+      
+      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 "dimphy.h"      
+#include "temps.h"      
+c 
+c Input:
+c ------
+      REAL*8  r_day                   ! Day of integration
+      LOGICAL first                 ! First timestep 
+                                    ! (and therefore initialization necessary)
+c      
+c Output:      
+c -------     
+      REAL*8  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
+      parameter (ny=jjm+1)
+      
+      INTEGER ismaller
+      INTEGER idec1, idec2 ! The two decadal data read ini
+      CHARACTER*4 cyear
+      
+      INTEGER im, day1, day2, im2
+      REAL*8 so4_1(iim, jjm+1, klev, 12)
+      REAL*8 so4_2(iim, jjm+1, klev, 12)   ! The sulfate distributions
+      
+      REAL*8 so4(klon, klev, 12)  ! SO4 in right dimension
+      SAVE so4
+      REAL*8 so4_out(klon, klev)
+      SAVE so4_out
+      
+      LOGICAL lnewday 
+      LOGICAL lonlyone
+      PARAMETER (lonlyone=.FALSE.)
+
+      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
+
+      ! Read in data:
+      ! a) from actual 10-yr-period
+
+      idec1 = (iyr-1900)/10
+      IF (idec1.LT.10) THEN
+         cyear='19'//char(idec1+48)//'0'
+      ELSE         
+         cyear='20'//char(idec1-10+48)//'0'
+      ENDIF
+      CALL getso4fromfile(cyear, so4_1)
+
+      
+      ! If to read two decades:
+      IF (.NOT.lonlyone) THEN
+      idec2=idec1+1
+         
+      ! b) from the next following one
+      IF (idec2.LT.10) THEN
+         cyear='19'//char(idec2+48)//'0'
+      ELSE
+         cyear='20'//char(idec2-10+48)//'0'
+      ENDIF
+      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-1900-10*idec1)/10.
+     .                 * (so4_1(i,j,k,it) - so4_2(i,j,k,it))
+               ENDDO
+            ENDDO
+         ENDDO
+      ENDDO                           
+      
+      ENDIF !lonlyone
+      
+      ! 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-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
+
+      
+      ! The sulfate concentration [molec cm-3] is read in. 
+      ! Convert it into mass [ug SO4/m3]
+      ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
+      DO k=1,klev
+         DO i=1,klon
+            sulfate(i,k) = sulfate(i,k)*masse_so4
+     .           /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?)
+      
+      RETURN
+      END
+
+      
+      
+      
+      
+c-----------------------------------------------------------------------------
+c Read in /calculate pre-industrial values of sulfate      
+c-----------------------------------------------------------------------------
+      
+      SUBROUTINE readsulfate_preind (r_day, first, pi_sulfate)
+      
+      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 "dimphy.h"      
+#include "temps.h"      
+c 
+c Input:
+c ------
+      REAL*8  r_day                   ! Day of integration
+      LOGICAL first                 ! First timestep 
+                                    ! (and therefore initialization necessary)
+c      
+c Output:      
+c -------     
+      REAL*8  pi_sulfate (klon, klev)  ! Number conc. sulfate (monthly mean data, 
+                                  !  from file)
+c      
+c Local Variables:
+c ----------------      
+      INTEGER i, ig, k, it
+      INTEGER j, iday, ny, iyr
+      parameter (ny=jjm+1)
+      
+      INTEGER im, day1, day2, im2, ismaller
+      REAL*8 pi_so4_1(iim, jjm+1, klev, 12)
+      
+      REAL*8 pi_so4(klon, klev, 12)  ! SO4 in right dimension
+      SAVE pi_so4
+      REAL*8 pi_so4_out(klon, klev)
+      SAVE pi_so4_out
+      
+      CHARACTER*4 cyear
+      LOGICAL lnewday
+
+      
+
+      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-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
+
+      
+      ! The sulfate concentration [molec cm-3] is read in. 
+      ! Convert it into mass [ug SO4/m3]
+      ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
+      DO k=1,klev
+         DO i=1,klon
+            pi_sulfate(i,k) = pi_sulfate(i,k)*masse_so4
+     .           /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?
+      RETURN
+      END
+
+      
+      
+      
+      
+      
+      
+      
+      
+      
+c-----------------------------------------------------------------------------
+c Routine for reading SO4 data from files
+c-----------------------------------------------------------------------------
+            
+
+      SUBROUTINE getso4fromfile (cyr, so4)
+#include "netcdf.inc"
+#include "dimensions.h"      
+#include "dimphy.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*8 so4mth(iim, ny, klev)
+c      REAL*8 so4mth(klev, ny, iim)
+      REAL*8 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
+c         STATUS = NF_INQ_VARID (NCID, VARMONTHS(i), VARID(i))
+         IF (STATUS .NE. NF_NOERR) write (*,*) 'err in read ',status      
+         STATUS = NF_GET_VARA_DOUBLE
+     .    (NCID, VARID, START,COUNT, so4mth)
+         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)
+c                  so4(i,j,k,imth)=so4mth(k,j,i)
+               ENDDO
+            ENDDO
+         ENDDO
+      ENDDO
+      
+      STATUS = NF_CLOSE(NCID)
+      END ! subroutine getso4fromfile
+      
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
