source: LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90 @ 1176

Last change on this file since 1176 was 1151, checked in by jghattas, 15 years ago
  • readaerosol : ENDIF mal place. Bug existant depuis la merge avec la branche LMDZ4_V3_patches. Le bug existe sur cette branche.
  • etat0_netcdf : manque des variables dans call conf_phys.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.0 KB
Line 
1!
2! $Header$
3!
4SUBROUTINE readaerosol (id_aero, r_day, first, massvar_p)
5 
6  USE dimphy, ONLY : klev
7  USE mod_grid_phy_lmdz, klon=>klon_glo
8  USE mod_phys_lmdz_para
9 
10  IMPLICIT NONE
11     
12! Content:
13! --------
14! This routine reads in monthly mean values of massvar aerosols and
15! returns a linearly interpolated dayly-mean field.     
16!
17!
18! Author:
19! -------
20! Johannes Quaas (quaas@lmd.jussieu.fr)
21! Celine Deandreis & Anne Cozic
22! 19/02/09
23!
24 
25!     
26! Include-files:
27! --------------     
28  INCLUDE "YOMCST.h"
29  INCLUDE "chem.h"     
30  INCLUDE "dimensions.h"     
31  INCLUDE "temps.h"     
32  INCLUDE "clesphys.h"
33  INCLUDE "iniprint.h"
34!
35! Input:
36! ------
37  INTEGER, INTENT(in) :: id_aero
38  REAL, INTENT(in)    :: r_day        ! Day of integration
39  LOGICAL, INTENT(in) :: first        ! First timestep
40                                      ! (and therefore initialization necessary)
41!     
42! Output:     
43! -------     
44  REAL, INTENT(out) ::   massvar_p(klon_omp,klev) ! Mass of massvar (monthly mean data,from file) [ug AIBCM/m3]
45
46!     
47! Local Variables:
48! ----------------     
49  INTEGER                         ::  i, ig, k, it
50  INTEGER                         :: j, iday, iyr, iyr1, iyr2
51  INTEGER                         :: im, day1, day2, im2
52  INTEGER, PARAMETER              :: ny=jjm+1
53  CHARACTER(len=4)                :: cyear
54  CHARACTER(len=7),DIMENSION(8)   :: name_aero
55  REAL                            :: var_1(iim, jjm+1, klev, 12)
56  REAL, DIMENSION(klon,klev)      ::  massvar
57  REAL, DIMENSION(iim, jjm+1, klev, 12)        :: var_2  ! The massvar distributions
58  REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE  :: var ! VAR in right dimension
59  REAL, ALLOCATABLE, DIMENSION(:,:,:),SAVE     :: var_out
60!$OMP THREADPRIVATE(var,var_out)
61
62  LOGICAL            :: lnewday
63  LOGICAL, PARAMETER :: lonlyone=.FALSE.
64  LOGICAL,SAVE       :: first2=.TRUE.
65!$OMP THREADPRIVATE(first2)
66
67! Variable pour definir a partir d'un numero d'aerosol, son nom
68  name_aero(1) = "SSSSM  "
69  name_aero(2) = "ASSSM  "
70  name_aero(3) = "ASBCM  "
71  name_aero(4) = "ASPOMM "
72  name_aero(5) = "SO4    "
73  name_aero(6) = "CIDUSTM"
74  name_aero(7) = "AIBCM  "
75  name_aero(8) = "AIPOMM "
76
77!$OMP MASTER
78  IF (first2) THEN
79      ALLOCATE( var(klon, klev, 12,8) )
80      ALLOCATE( var_out(klon, klev,8))
81      first2=.FALSE.
82
83      IF (aer_type /= 'actuel  ' .AND. aer_type /= 'preind  ' .AND.   &
84           aer_type /= 'scenario') THEN
85         WRITE(lunout,*)' *** Warning ***'
86         WRITE(lunout,*)'Option aer_type non prevu pour les aerosols = ',       &
87              aer_type
88         CALL abort_gcm('readaerosol','Error : aer_type parameter not accepted',1)
89      ENDIF
90   ENDIF
91
92  IF (is_mpi_root) THEN
93     
94      iday = INT(r_day)
95      ! Get the year of the run
96      iyr  = iday/360
97      ! Get the day of the actual year:
98      iday = iday-iyr*360
99      ! 0.02 is about 0.5/24, namly less than half an hour
100      lnewday = (r_day-FLOAT(iday).LT.0.02)
101     
102      !     ---------------------------------------------
103      !     All has to be done only, if a new day begins!       
104      !     ---------------------------------------------
105         
106      IF (lnewday.OR.first) THEN
107         
108          im = iday/30 +1     ! the actual month
109          ! annee_ref is the initial year (defined in temps.h)
110          iyr = iyr + annee_ref
111         
112          ! Do I have to read new data? (Is this the first day of a year?)
113          IF (first.OR.iday.EQ.1.) THEN
114              ! Initialize values
115              DO it=1,12
116                DO k=1,klev
117                  DO i=1,klon
118                    var(i,k,it,id_aero)=0.
119                  ENDDO
120                ENDDO
121              ENDDO
122             
123
124              IF (aer_type == 'actuel  ') then
125
126                 cyear='1980'
127                 CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero))
128
129              ELSE IF (aer_type == 'preind  ') THEN
130
131                 cyear='.nat'
132                 CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero))
133
134              ELSE ! aer_type=scenario
135
136                 IF (iyr .LT. 1850) THEN
137                    cyear='.nat'
138                    WRITE(lunout,*) 'get_aero  iyr=', iyr,'   ',cyear
139                    CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero))
140
141                 ELSE IF (iyr .GE. 2100) THEN
142                    cyear='2100'
143                    WRITE(lunout,*) 'get_aero  iyr=', iyr,'   ',cyear
144                    CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero))
145
146                 ELSE ! Read data from 2 decades
147                    ! a) from actual 10-yr-period
148                    IF (iyr.LT.1900) THEN
149                       iyr1 = 1850
150                       iyr2 = 1900
151                    ELSE IF (iyr.GE.1900.AND.iyr.LT.1920) THEN
152                       iyr1 = 1900
153                       iyr2 = 1920
154                    ELSE
155                       iyr1 = INT(iyr/10)*10
156                       iyr2 = INT(1+iyr/10)*10
157                    ENDIF
158                   
159                    WRITE(cyear,'(I4)') iyr1
160                    WRITE(lunout,*) 'get_aero  iyr=', iyr,'   ',cyear
161                    CALL get_aero_fromfile(cyear, var_1, name_aero(id_aero))
162                   
163                    ! If to read two decades:
164                    IF (.NOT.lonlyone) THEN
165                       
166                       ! b) from the next following one
167                       WRITE(cyear,'(I4)') iyr2
168                       WRITE(lunout,*) 'get_aero  iyr=', iyr,'   ',cyear
169                       
170                       CALL get_aero_fromfile(cyear, var_2, name_aero(id_aero))
171                       
172                       ! Interpolate linarily to the actual year:
173                       DO it=1,12
174                          DO k=1,klev
175                             DO j=1,jjm
176                                DO i=1,iim
177                                   var_1(i,j,k,it) = &
178                                        var_1(i,j,k,it) - FLOAT(iyr-iyr1)/FLOAT(iyr2-iyr1) * &
179                                        (var_1(i,j,k,it) - var_2(i,j,k,it))
180                                ENDDO
181                             ENDDO
182                          ENDDO
183                       ENDDO
184                   
185                    ENDIF ! lonlyone
186                 ENDIF ! iyr .LT. 1850       
187              ENDIF ! aer_type
188               
189              ! Transform the horizontal 2D-field into the physics-field
190              ! (Also the levels and the latitudes have to be inversed)
191             
192              DO it=1,12
193                DO k=1, klev         
194                  ! a) at the poles, use the zonal mean:
195                  DO i=1,iim
196                    ! North pole
197                    var(1,k,it,id_aero) = &
198                       var(1,k,it,id_aero)+ var_1(i,jjm+1,klev+1-k,it)
199                    ! South pole
200                    var(klon,k,it,id_aero)= &
201                       var(klon,k,it,id_aero)+ var_1(i,1,klev+1-k,it)
202                  ENDDO
203                  var(1,k,it,id_aero)   = var(1,k,it,id_aero)/FLOAT(iim)
204                  var(klon,k,it,id_aero)= var(klon,k,it,id_aero)/FLOAT(iim)
205                   
206                  ! b) the values between the poles:
207                  ig=1
208                  DO j=2,jjm
209                    DO i=1,iim
210                      ig=ig+1
211
212                      IF (ig.GT.klon)  THEN
213                          WRITE(lunout,*) 'Attention dans readaerosol', &
214                             name_aero(id_aero), 'ig > klon', ig, klon
215                      ENDIF
216
217                      var(ig,k,it,id_aero) =  var_1(i,jjm+1+1-j,klev+1-k,it)
218
219                    ENDDO
220                  ENDDO
221                  IF (ig.NE.klon-1) THEN
222                      PRINT *,'Error in readaerosol', name_aero(id_aero), 'conversion'
223                      CALL abort_gcm('readaerosol','Error in readaerosol 1',1)
224                  ENDIF
225                ENDDO         ! Loop over k (vertical)
226              ENDDO           ! Loop over it (months)
227             
228          ENDIF               ! Had to read new data?
229           
230     
231          ! Interpolate to actual day:
232          IF (iday.LT.im*30-15) THEN         
233              ! in the first half of the month use month before and actual month
234              im2=im-1
235              day2 = im2*30-15
236              day1 = im2*30+15
237              IF (im2.LE.0) THEN
238                  ! the month is january, thus the month before december
239                  im2=12
240              ENDIF
241              DO k=1,klev
242                DO i=1,klon
243                  massvar(i,k) = &
244                     var(i,k,im2,id_aero)- FLOAT(iday-day2)/FLOAT(day1-day2) * &
245                     (var(i,k,im2,id_aero) - var(i,k,im,id_aero))
246
247                  IF (massvar(i,k).LT.0.) THEN
248                      IF (iday-day2.LT.0.) WRITE(lunout,*) 'iday-day2',iday-day2
249                      IF (var(i,k,im2,id_aero) - var(i,k,im,id_aero).LT.0.) THEN
250                          PRINT*, name_aero(id_aero),'(i,k,im2)- ', &
251                             name_aero(id_aero),'(i,k,im)',         &
252                             var(i,k,im2,id_aero) - var(i,k,im,id_aero)
253                      ENDIF
254
255                      IF (day1-day2.LT.0.) WRITE(lunout,*) 'day1-day2',day1-day2
256                      WRITE(lunout,*) 'stop',name_aero(id_aero)
257                      CALL abort_gcm('readaerosol','Error in interpolation 2',1)
258                  ENDIF
259                ENDDO
260              ENDDO
261          ELSE
262              ! the second half of the month
263              im2=im+1
264              IF (im2.GT.12) THEN
265                  ! the month is december, the following thus january
266                  im2=1
267              ENDIF
268              day2 = im*30-15
269              day1 = im*30+15
270              DO k=1,klev
271                DO i=1,klon
272                  massvar(i,k) = &
273                     var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
274                     (var(i,k,im2,id_aero) - var(i,k,im,id_aero))
275
276                  IF (massvar(i,k).LT.0.) THEN
277                      IF (iday-day2.LT.0.) &
278                         WRITE(lunout,*) 'iday-day2',iday-day2
279                      IF (var(i,k,im2,id_aero) -  var(i,k,im,id_aero).LT.0.) THEN
280                          WRITE(lunout,*) name_aero(id_aero),'(i,k,im2)-', &
281                             name_aero(id_aero),'(i,k,im)',           &
282                             var(i,k,im2,id_aero) - var(i,k,im,id_aero)
283                      ENDIF
284                      IF (day1-day2.LT.0.) &
285                         WRITE(lunout,*) 'day1-day2',day1-day2
286                      WRITE(lunout,*) 'stop',name_aero(id_aero)
287                      CALL abort_gcm('readaerosol','Error in interpolation 3',1)
288                  ENDIF
289                ENDDO
290              ENDDO
291          ENDIF
292         
293     
294          ! The massvar concentration [molec cm-3] is read in.
295          ! Convert it into mass [ug VAR/m3]
296          ! masse_var in [g/mol], n_avogadro in [molec/mol]
297          ! The sulfate mass [ug VAR/m3] is read in.
298          DO k=1,klev
299            DO i=1,klon
300              var_out(i,k,id_aero) = massvar(i,k)
301              IF (var_out(i,k,id_aero).LT.0) THEN
302                  PRINT*, 'Attention il y a un probleme dans readaerosol', &
303                     name_aero(id_aero), 'la masse est negative'
304                  CALL abort_gcm('readaerosol','Error in readaerosol 4',1)
305              ENDIF
306            ENDDO
307          ENDDO
308      ELSE                    ! if no new day, use old data:
309          DO k=1,klev
310            DO i=1,klon
311              massvar(i,k) = var_out(i,k,id_aero)
312              IF (var_out(i,k,id_aero).LT.0)   THEN
313                  PRINT*, 'Attention il y a un probleme dans readaerosol', &
314                     name_aero(id_aero), 'la masse est negative'
315                  CALL abort_gcm('readaerosol','Error in readaerosol 5',1)
316              ENDIF
317            ENDDO
318          ENDDO
319         
320           
321      ENDIF                   ! Did I have to do anything (was it a new day?)
322     
323  ENDIF                      ! phy_rank==0
324 
325!$OMP END MASTER
326  CALL Scatter(massvar,massvar_p)             
327
328END SUBROUTINE readaerosol
329     
330     
331     
332     
333     
334!-----------------------------------------------------------------------------
335! Read in /calculate pre-industrial values of sulfate     
336!-----------------------------------------------------------------------------
337     
338SUBROUTINE readaerosol_preind (id_aero, r_day, first, pi_massvar_p)
339 
340  USE dimphy, ONLY : klev
341  USE mod_grid_phy_lmdz, klon=>klon_glo
342  USE mod_phys_lmdz_para
343  IMPLICIT NONE
344     
345  ! Content:
346  ! --------
347  ! This routine reads in monthly mean values of massvar aerosols and
348  ! returns a linearly interpolated dayly-mean field.     
349  !
350  ! It does so for the preindustriel values of the massvar, to a large part
351  ! analogous to the routine readmassvar above.     
352  !
353  ! Only Pb: Variables must be saved and don t have to be overwritten!
354  !     
355  ! Author:
356  ! -------
357  ! Celine Deandreis & Anne Cozic (LSCE) 2009
358  ! Johannes Quaas (quaas@lmd.jussieu.fr)
359  ! 26/06/01
360  !
361  ! Modifications:
362  ! --------------
363  ! see above
364
365     
366  ! Include-files:
367  ! --------------     
368
369  INCLUDE "YOMCST.h"
370  INCLUDE "chem.h"     
371  INCLUDE "dimensions.h"     
372  INCLUDE "temps.h"     
373  INCLUDE "iniprint.h"
374 
375  ! Input:
376  ! ------
377  INTEGER, INTENT(in) ::  id_aero
378  REAL, INTENT(in)    ::  r_day          ! Day of integration
379  LOGICAL, INTENT(in) ::  first          ! First timestep (and therefore initialization necessary)
380
381
382  !     
383  ! Output:     
384  ! -------     
385  REAL, DIMENSION(klon_omp,klev), INTENT(out) ::  pi_massvar_p ! Number conc. massvar (monthly mean data, from file)
386 
387
388  !     
389  ! Local Variables:
390  ! ----------------     
391  INTEGER                                      :: i, ig, k, it
392  INTEGER                                      :: j, iday, iyr
393  INTEGER, PARAMETER                           :: ny=jjm+1
394  INTEGER                                      :: im, day1, day2, im2
395
396  REAL, DIMENSION(iim, jjm+1, klev, 12)        :: pi_var_1
397  REAL, DIMENSION(klon,klev)                   :: pi_massvar
398  REAL, ALLOCATABLE, DIMENSION(:,:,:,:), SAVE  :: pi_var ! VAR in right dimension
399  REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE    :: pi_var_out
400!$OMP THREADPRIVATE(pi_var,pi_var_out)           
401
402  CHARACTER(len=4)                :: cyear
403  CHARACTER(len=7),DIMENSION(8)   :: name_aero
404  LOGICAL                         :: lnewday
405  LOGICAL,SAVE                    :: first2=.TRUE.
406!$OMP THREADPRIVATE(first2)
407 
408
409! Variable pour definir a partir d'un numero d'aerosol, son nom
410  name_aero(1) = "SSSSM  "
411  name_aero(2) = "ASSSM  "
412  name_aero(3) = "ASBCM  "
413  name_aero(4) = "ASPOMM "
414  name_aero(5) = "SO4    "
415  name_aero(6) = "CIDUSTM"
416  name_aero(7) = "AIBCM  "
417  name_aero(8) = "AIPOMM "
418
419
420!$OMP MASTER
421  IF (first2) THEN
422      ALLOCATE( pi_var(klon, klev, 12,8) )
423      ALLOCATE( pi_var_out(klon, klev,8))
424      first2=.FALSE.
425  ENDIF
426
427  IF (is_mpi_root) THEN
428      iday = INT(r_day)
429      ! Get the year of the run
430      iyr  = iday/360
431      ! Get the day of the actual year:
432      iday = iday-iyr*360
433      ! 0.02 is about 0.5/24, namly less than half an hour
434      lnewday = (r_day-FLOAT(iday).LT.0.02)
435     
436      !     ---------------------------------------------
437      !     All has to be done only, if a new day begins!       
438      !     ---------------------------------------------
439     
440      IF (lnewday.OR.first) THEN
441         
442          im = iday/30 +1     ! the actual month
443          ! annee_ref is the initial year (defined in temps.h)
444          iyr = iyr + annee_ref     
445
446          IF (first) THEN
447              cyear='.nat'
448              CALL get_aero_fromfile(cyear,pi_var_1, name_aero(id_aero))
449             
450              ! Transform the horizontal 2D-field into the physics-field
451              ! (Also the levels and the latitudes have to be inversed)
452             
453              ! Initialize field
454              DO it=1,12
455                DO k=1,klev
456                  DO i=1,klon
457                    pi_var(i,k,it,id_aero)=0.
458                  ENDDO
459                ENDDO
460              ENDDO
461               
462              WRITE(lunout,*) 'preind: finished reading', FLOAT(iim)
463              DO it=1,12
464                DO k=1, klev         
465                  ! a) at the poles, use the zonal mean:
466                  DO i=1,iim
467                    ! North pole
468                    pi_var(1,k,it,id_aero)    = &
469                       pi_var(1,k,it,id_aero) + pi_var_1(i,jjm+1,klev+1-k,it)
470                    ! South pole
471                    pi_var(klon,k,it,id_aero) = &
472                       pi_var(klon,k,it,id_aero) + pi_var_1(i,1,klev+1-k,it)
473                  ENDDO
474                  pi_var(1,k,it,id_aero)    = pi_var(1,k,it,id_aero)/FLOAT(iim)
475                  pi_var(klon,k,it,id_aero) = pi_var(klon,k,it,id_aero)/FLOAT(iim)
476                 
477                  ! b) the values between the poles:
478                  ig=1
479                  DO j=2,jjm
480                    DO i=1,iim
481                      ig=ig+1
482                      IF (ig.GT.klon)  THEN
483                          WRITE(lunout,*) 'Attention dans readaerosol_preind', &
484                             name_aero(id_aero), 'ig > klon', ig, klon
485                      ENDIF
486                      pi_var(ig,k,it,id_aero) = &
487                         pi_var_1(i,jjm+1+1-j,klev+1-k,it)
488                    ENDDO
489                  ENDDO
490                  IF (ig.NE.klon-1) THEN
491                      WRITE(lunout,*) 'Error in readaerosol_preind (', name_aero(id_aero), ' conversion)'
492                      CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 1',1)
493                  ENDIF
494                ENDDO         ! Loop over k (vertical)
495              ENDDO             ! Loop over it (months)
496               
497          ENDIF                ! Had to read new data?
498           
499           
500                                ! Interpolate to actual day:
501          IF (iday.LT.im*30-15) THEN         
502              ! in the first half of the month use month before and actual month
503              im2=im-1
504              day1 = im2*30+15
505              day2 = im2*30-15
506              IF (im2.LE.0) THEN 
507                  ! the month is january, thus the month before december
508                  im2=12
509              ENDIF
510              DO k=1,klev
511                DO i=1,klon
512                  pi_massvar(i,k) = &
513                     pi_var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
514                     (pi_var(i,k,im2,id_aero) -  pi_var(i,k,im,id_aero))
515
516                  IF (pi_massvar(i,k).LT.0.) THEN
517                      IF (iday-day2.LT.0.)  WRITE(lunout,*) 'iday-day2',iday-day2
518                      IF (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero).LT.0.) THEN
519                          WRITE(lunout,*) 'pi_',name_aero(id_aero),'(i,k,im2) - pi_', &
520                             name_aero(id_aero),'(i,k,im)', &
521                             pi_var(i,k,im2,id_aero) -pi_var(i,k,im,id_aero)
522                      ENDIF
523                      IF (day1-day2.LT.0.)WRITE(lunout,*) 'day1-day2',day1-day2
524                      PRINT *, 'pi_',name_aero(id_aero)
525                      CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 2',1)
526                  ENDIF
527                ENDDO
528              ENDDO
529          ELSE
530
531              ! the second half of the month
532              im2=im+1
533              day1 = im*30+15
534              IF (im2.GT.12) THEN
535                  ! the month is december, the following thus january
536                  im2=1
537              ENDIF
538              day2 = im*30-15
539             
540              DO k=1,klev
541                DO i=1,klon
542                  pi_massvar(i,k) = &
543                     pi_var(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
544                     (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero))
545                  IF (pi_massvar(i,k).LT.0.) THEN
546                      IF (iday-day2.LT.0.)  WRITE(lunout,*) 'iday-day2',iday-day2
547                      IF (pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero).LT.0.) THEN
548                          WRITE(lunout,*) 'pi_', name_aero(id_aero),'(i,k,im2) - pi_',&
549                             name_aero(id_aero), '(i,k,im)',&
550                             pi_var(i,k,im2,id_aero) - pi_var(i,k,im,id_aero)
551                      ENDIF
552                      IF (day1-day2.LT.0.) &
553                         WRITE(lunout,*) 'day1-day2',day1-day2
554                      PRINT *, 'pi_',name_aero(id_aero)
555                      CALL abort_gcm('readaerosol_preind','Error in readaerosol_preind 3',1)
556                  ENDIF
557                ENDDO
558              ENDDO
559          ENDIF
560         
561         
562          ! The massvar concentration [molec cm-3] is read in.
563          ! Convert it into mass [ug VAR/m3]
564          ! masse_var in [g/mol], n_avogadro in [molec/mol]
565          DO k=1,klev
566            DO i=1,klon
567              pi_var_out(i,k,id_aero) = pi_massvar(i,k)
568            ENDDO
569          ENDDO
570         
571      ELSE                   ! If no new day, use old data:
572          DO k=1,klev
573            DO i=1,klon
574              pi_massvar(i,k) = pi_var_out(i,k,id_aero)           
575            ENDDO
576          ENDDO
577      ENDIF                  ! Was this the beginning of a new day?
578   ENDIF                     ! is_mpi_root
579     
580!$OMP END MASTER
581   CALL Scatter(pi_massvar,pi_massvar_p)           
582     
583 END SUBROUTINE readaerosol_preind
584     
585     
586     
587!-----------------------------------------------------------------------------
588! Routine for reading VAR data from files
589!-----------------------------------------------------------------------------
590           
591
592SUBROUTINE get_aero_fromfile (cyr, var, varname)
593  USE dimphy
594  IMPLICIT NONE
595     
596  INCLUDE "netcdf.inc"
597  INCLUDE "dimensions.h"     
598  INCLUDE "iniprint.h"
599
600  CHARACTER(len=4), INTENT(in)                       ::  cyr
601  REAL, DIMENSION(iim, jjm+1, klev, 12), INTENT(out) ::  var
602  CHARACTER(len=*), INTENT(in)                       ::  varname
603 
604
605  CHARACTER(len=30)     :: fname
606  CHARACTER(len=30)     :: cvar
607  INTEGER, DIMENSION(3) :: START, COUNT
608  INTEGER               :: STATUS, NCID, VARID
609  INTEGER               :: imth, i, j, k
610  INTEGER, PARAMETER    :: ny=jjm+1
611  REAL, DIMENSION(iim, ny, klev) :: varmth
612!
613!
614
615  fname = TRIM(varname)//'.run'//cyr//'.cdf'
616 
617  WRITE(lunout,*) 'reading ', TRIM(fname)
618  STATUS = NF_OPEN (TRIM(fname), NF_NOWRITE, NCID)
619  IF (STATUS .NE. NF_NOERR) THEN
620      WRITE(lunout,*) 'err in open ',status
621      CALL abort_gcm('get_aero_fromfile','Error in opening file',1)
622  ENDIF
623  DO imth=1, 12
624    IF (imth.EQ.1) THEN
625        cvar=TRIM(varname)//'JAN'
626    ELSEIF (imth.EQ.2) THEN
627        cvar=TRIM(varname)//'FEB'
628    ELSEIF (imth.EQ.3) THEN
629        cvar=TRIM(varname)//'MAR'
630    ELSEIF (imth.EQ.4) THEN
631        cvar=TRIM(varname)//'APR'
632    ELSEIF (imth.EQ.5) THEN
633        cvar=TRIM(varname)//'MAY'
634    ELSEIF (imth.EQ.6) THEN
635        cvar=TRIM(varname)//'JUN'
636    ELSEIF (imth.EQ.7) THEN
637        cvar=TRIM(varname)//'JUL'
638    ELSEIF (imth.EQ.8) THEN
639        cvar=TRIM(varname)//'AUG'
640    ELSEIF (imth.EQ.9) THEN
641        cvar=TRIM(varname)//'SEP'
642    ELSEIF (imth.EQ.10) THEN
643        cvar=TRIM(varname)//'OCT'
644    ELSEIF (imth.EQ.11) THEN
645        cvar=TRIM(varname)//'NOV'
646    ELSEIF (imth.EQ.12) THEN
647        cvar=TRIM(varname)//'DEC'
648    ENDIF
649    start(1)=1
650    start(2)=1
651    start(3)=1
652    COUNT(1)=iim
653    COUNT(2)=ny
654    COUNT(3)=klev
655    STATUS = NF_INQ_VARID (NCID, TRIM(cvar), VARID)
656    WRITE(lunout,*) ncid,imth,TRIM(cvar), varid
657   
658    IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in read ',status     
659
660#ifdef NC_DOUBLE
661    status = NF_GET_VARA_DOUBLE(NCID, VARID, START, COUNT,varmth)
662#else
663    status = NF_GET_VARA_REAL(NCID, VARID, START, COUNT, varmth)
664#endif
665    IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in read data',status
666         
667    DO k=1,klev
668      DO j=1,jjm+1
669        DO i=1,iim
670          IF (varmth(i,j,k).LT.0.) THEN
671              WRITE(lunout,*) 'this is shit'
672              WRITE(lunout,*) varname,'(',i,j,k,') =',varmth(i,j,k)
673          ENDIF
674          var(i,j,k,imth)=varmth(i,j,k)
675        ENDDO
676      ENDDO
677    ENDDO
678  ENDDO
679 
680  STATUS = NF_CLOSE(NCID)
681  IF (STATUS .NE. NF_NOERR) WRITE(lunout,*) 'err in closing file',status
682 
683
684END SUBROUTINE get_aero_fromfile
685     
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
Note: See TracBrowser for help on using the repository browser.