/[lmdze]/trunk/Sources/phylmd/readsulfate.f
ViewVC logotype

Diff of /trunk/Sources/phylmd/readsulfate.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/phylmd/readsulfate.f revision 51 by guez, Tue Sep 20 09:14:34 2011 UTC trunk/phylmd/readsulfate.f revision 105 by guez, Thu Sep 4 10:40:24 2014 UTC
# Line 1  Line 1 
1  !  module readsulfate_m
 ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/readsulfate.F,v 1.2 2005/05/19 08:27:15 fairhead Exp $  
 !  
       SUBROUTINE readsulfate(r_day, first, sulfate)  
         
       use dimens_m        
       use dimphy        
       use temps        
       use SUPHEC_M  
       use chem        
       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 --------------      
 c  
 c Input:  
 c ------  
       REAL*8, intent(in):: r_day                   ! Day of integration  
       LOGICAL, intent(in):: 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, iyr1, iyr2  
       parameter (ny=jjm+1)  
         
       INTEGER ismaller  
 CJLD      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  
   
   
       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  
       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)  
   
       ENDIF  
   
       ! 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  
         
       ! 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  
   
         
 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?)  
         
       RETURN  
       END  
   
         
         
         
         
 c-----------------------------------------------------------------------------  
 c Read in /calculate pre-industrial values of sulfate        
 c-----------------------------------------------------------------------------  
         
       SUBROUTINE readsulfate_preind (r_day, first, pi_sulfate)  
         
       use dimens_m        
       use dimphy        
       use temps        
       use SUPHEC_M  
       use chem        
       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 --------------      
 c  
 c Input:  
 c ------  
       REAL*8, intent(in)::  r_day                   ! Day of integration  
       LOGICAL, intent(in):: 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  
   
         
 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?  
       RETURN  
       END  
   
         
         
         
         
         
         
         
         
         
 c-----------------------------------------------------------------------------  
 c Routine for reading SO4 data from files  
 c-----------------------------------------------------------------------------  
               
   
       SUBROUTINE getso4fromfile (cyr, so4)  
       use dimens_m        
       use dimphy  
       include "netcdf.inc"  
       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  
         
   
   
   
   
   
   
   
   
   
2    
3      IMPLICIT none
4    
5    contains
6    
7      SUBROUTINE readsulfate(r_day, first, sulfate)
8    
9        ! From LMDZ4/libf/phylmd/readsulfate.F, version 1.2 2005/05/19
10        ! 08:27:15 fairhead
11    
12        ! This routine reads in monthly mean values of sulfate aerosols and
13        ! returns a linearly interpolated daily-mean field.
14    
15        ! Author: Johannes Quaas (quaas@lmd.jussieu.fr)
16        ! 26/04/01
17    
18        ! Modifications:
19        ! 21/06/01: Make integrations of more than one year possible ;-)
20        ! ATTENTION!! runs are supposed to start with Jan, 1. 1930
21        ! (rday=1)
22    
23        ! 27/06/01: Correction: The model always has 360 days per year!
24        ! 27/06/01: SO4 concentration rather than mixing ratio
25        ! 27/06/01: 10yr-mean-values to interpolate
26        ! 20/08/01: Correct the error through integer-values in interpolations
27        ! 21/08/01: Introduce flag to read in just one decade
28    
29        USE dimens_m, ONLY: iim, jjm
30        USE dimphy, ONLY: klev, klon
31        use getso4fromfile_m, only: getso4fromfile
32        USE temps, ONLY: annee_ref
33    
34        ! Input:
35    
36        real, intent(in):: r_day                   ! Day of integration
37        LOGICAL, intent(in):: first                 ! First timestep
38        ! (and therefore initialization necessary)
39    
40        ! Output:
41    
42        real sulfate (klon, klev)  ! Mass of sulfate (monthly mean data,
43        !  from file) [ug SO4/m3]
44    
45        ! Local Variables:
46    
47        INTEGER i, ig, k, it
48        INTEGER j, iday, ny, iyr, iyr1, iyr2
49        parameter (ny=jjm+1)
50    
51        CHARACTER(len=4) cyear
52    
53        INTEGER im, day1, day2, im2
54        double precision so4_1(iim, jjm+1, klev, 12)
55        double precision so4_2(iim, jjm+1, klev, 12)   ! The sulfate distributions
56    
57        double precision so4(klon, klev, 12)  ! SO4 in right dimension
58        SAVE so4
59        double precision so4_out(klon, klev)
60        SAVE so4_out
61    
62        LOGICAL lnewday
63        LOGICAL lonlyone
64        PARAMETER (lonlyone=.FALSE.)
65    
66        !--------------------------------------------------------------------
67    
68        iday = INT(r_day)
69    
70        ! Get the year of the run
71        iyr  = iday/360
72    
73        ! Get the day of the actual year:
74        iday = iday-iyr*360
75    
76        ! 0.02 is about 0.5/24, namly less than half an hour
77        lnewday = (r_day-FLOAT(iday).LT.0.02)
78    
79        ! All has to be done only, if a new day begins!
80    
81        IF (lnewday.OR.first) THEN
82           im = iday/30 +1 ! the actual month
83           ! annee_ref is the initial year (defined in temps.h)
84           iyr = iyr + annee_ref
85    
86           ! Do I have to read new data? (Is this the first day of a year?)
87           IF (first.OR.iday.EQ.1.) THEN
88              ! Initialize values
89              DO it=1,12
90                 DO k=1,klev
91                    DO i=1,klon
92                       so4(i,k,it)=0.
93                    ENDDO
94                 ENDDO
95              ENDDO
96    
97              IF (iyr .lt. 1850) THEN
98                 cyear='.nat'
99                 WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear
100                 CALL getso4fromfile(cyear, so4_1)
101              ELSE IF (iyr .ge. 2100) THEN
102                 cyear='2100'
103                 WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear
104                 CALL getso4fromfile(cyear, so4_1)
105              ELSE
106    
107                 ! Read in data:
108                 ! a) from actual 10-yr-period
109    
110                 IF (iyr.LT.1900) THEN
111                    iyr1 = 1850
112                    iyr2 = 1900
113                 ELSE IF (iyr.ge.1900.and.iyr.lt.1920) THEN
114                    iyr1 = 1900
115                    iyr2 = 1920
116                 ELSE
117                    iyr1 = INT(iyr/10)*10
118                    iyr2 = INT(1+iyr/10)*10
119                 ENDIF
120                 WRITE(cyear,'(I4)') iyr1
121                 WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear
122                 CALL getso4fromfile(cyear, so4_1)
123    
124                 ! If to read two decades:
125                 IF (.NOT.lonlyone) THEN
126    
127                    ! b) from the next following one
128                    WRITE(cyear,'(I4)') iyr2
129                    WRITE(*,*) 'getso4  iyr=', iyr,'   ',cyear
130                    CALL getso4fromfile(cyear, so4_2)
131    
132                 ENDIF
133    
134                 ! Interpolate linarily to the actual year:
135                 DO it=1,12
136                    DO k=1,klev
137                       DO j=1,jjm
138                          DO i=1,iim
139                             so4_1(i,j,k,it)=so4_1(i,j,k,it) &
140                                  - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1) &
141                                  * (so4_1(i,j,k,it) - so4_2(i,j,k,it))
142                          ENDDO
143                       ENDDO
144                    ENDDO
145                 ENDDO
146    
147              ENDIF !lonlyone
148    
149              ! Transform the horizontal 2D-field into the physics-field
150              ! (Also the levels and the latitudes have to be inversed)
151    
152              DO it=1,12
153                 DO k=1, klev
154                    ! a) at the poles, use the zonal mean:
155                    DO i=1,iim
156                       ! North pole
157                       so4(1,k,it)=so4(1,k,it)+so4_1(i,jjm+1,klev+1-k,it)
158                       ! South pole
159                       so4(klon,k,it)=so4(klon,k,it)+so4_1(i,1,klev+1-k,it)
160                    ENDDO
161                    so4(1,k,it)=so4(1,k,it)/FLOAT(iim)
162                    so4(klon,k,it)=so4(klon,k,it)/FLOAT(iim)
163    
164                    ! b) the values between the poles:
165                    ig=1
166                    DO j=2,jjm
167                       DO i=1,iim
168                          ig=ig+1
169                          if (ig.gt.klon) write (*,*) 'shit'
170                          so4(ig,k,it) = so4_1(i,jjm+1-j,klev+1-k,it)
171                       ENDDO
172                    ENDDO
173                    IF (ig.NE.klon-1) STOP 'Error in readsulfate (var conversion)'
174                 ENDDO ! Loop over k (vertical)
175              ENDDO ! Loop over it (months)
176    
177           ENDIF ! Had to read new data?
178    
179           ! Interpolate to actual day:
180           IF (iday.LT.im*30-15) THEN
181              ! in the first half of the month use month before and actual month
182              im2=im-1
183              day2 = im2*30-15
184              day1 = im2*30+15
185              IF (im2.LE.0) THEN
186                 ! the month is january, thus the month before december
187                 im2=12
188              ENDIF
189              DO k=1,klev
190                 DO i=1,klon
191                    sulfate(i,k) = so4(i,k,im2)   &
192                         - FLOAT(iday-day2)/FLOAT(day1-day2) &
193                         * (so4(i,k,im2) - so4(i,k,im))
194                    IF (sulfate(i,k).LT.0.) THEN
195                       IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
196                       IF (so4(i,k,im2) - so4(i,k,im).LT.0.) &
197                            write(*,*) 'so4(i,k,im2) - so4(i,k,im)', &
198                            so4(i,k,im2) - so4(i,k,im)
199                       IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
200                       stop 'sulfate'
201                    endif
202                 ENDDO
203              ENDDO
204           ELSE
205              ! the second half of the month
206              im2=im+1
207              IF (im2.GT.12) THEN
208                 ! the month is december, the following thus january
209                 im2=1
210              ENDIF
211              day2 = im*30-15
212              day1 = im*30+15
213              DO k=1,klev
214                 DO i=1,klon
215                    sulfate(i,k) = so4(i,k,im2)   &
216                         - FLOAT(iday-day2)/FLOAT(day1-day2) &
217                         * (so4(i,k,im2) - so4(i,k,im))
218                    IF (sulfate(i,k).LT.0.) THEN
219                       IF (iday-day2.LT.0.) write(*,*) 'iday-day2',iday-day2
220                       IF (so4(i,k,im2) - so4(i,k,im).LT.0.) &
221                            write(*,*) 'so4(i,k,im2) - so4(i,k,im)', &
222                            so4(i,k,im2) - so4(i,k,im)
223                       IF (day1-day2.LT.0.) write(*,*) 'day1-day2',day1-day2
224                       stop 'sulfate'
225                    endif
226                 ENDDO
227              ENDDO
228           ENDIF
229    
230           !JLD      ! The sulfate concentration [molec cm-3] is read in.
231           !JLD      ! Convert it into mass [ug SO4/m3]
232           !JLD      ! masse_so4 in [g/mol], n_avogadro in [molec/mol]
233           ! The sulfate mass [ug SO4/m3] is read in.
234           DO k=1,klev
235              DO i=1,klon
236                 !JLD            sulfate(i,k) = sulfate(i,k)*masse_so4
237                 !JLD     .           /n_avogadro*1.e+12
238                 so4_out(i,k) = sulfate(i,k)
239                 IF (so4_out(i,k).LT.0)  &
240                      stop 'WAS SOLL DER SCHEISS ? '
241              ENDDO
242           ENDDO
243        ELSE ! if no new day, use old data:
244           DO k=1,klev
245              DO i=1,klon
246                 sulfate(i,k) = so4_out(i,k)
247                 IF (so4_out(i,k).LT.0)  &
248                      stop 'WAS SOLL DER SCHEISS ? '
249              ENDDO
250           ENDDO
251        ENDIF ! Did I have to do anything (was it a new day?)
252    
253      END SUBROUTINE readsulfate
254    
255    end module readsulfate_m
Legend:
Removed from v.51  
changed lines
  Added in v.105
  ViewVC Help
Powered by ViewVC 1.1.21