source: LMDZ5/trunk/libf/phylmd/limit_read_mod.F90 @ 2769

Last change on this file since 2769 was 2768, checked in by dcugnet, 8 years ago

Fix a bug about limit.nc file checking.
The limit.nc number of days does not need to be checked every day ; it
is now checked only at the start of the run ; this avoids the gcm to
stop because of a mismatch between the file and model number of days
(in the case of an irregular calendar) on last day due to the fact
that the calendar is updated before the end of the day.
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.1 KB
Line 
1!
2! $Header$
3!
4MODULE limit_read_mod
5!
6! This module reads the fichier "limit.nc" containing fields for surface forcing.
7!
8! Module subroutines :
9!  limit_read_frac    : call limit_read_tot and return the fractions
10!  limit_read_rug_alb : return rugosity and albedo, if coupled ocean call limit_read_tot first
11!  limit_read_sst     : return sea ice temperature   
12!  limit_read_tot     : read limit.nc and store the fields in local modules variables
13!
14  IMPLICIT NONE
15
16  REAL, ALLOCATABLE, DIMENSION(:,:), SAVE, PRIVATE :: pctsrf
17!$OMP THREADPRIVATE(pctsrf)
18  REAL, ALLOCATABLE, DIMENSION(:),   SAVE, PRIVATE :: rugos
19!$OMP THREADPRIVATE(rugos)
20  REAL, ALLOCATABLE, DIMENSION(:),   SAVE, PRIVATE :: albedo
21!$OMP THREADPRIVATE(albedo) 
22  REAL, ALLOCATABLE, DIMENSION(:),   SAVE, PRIVATE :: sst
23!$OMP THREADPRIVATE(sst) 
24  LOGICAL,SAVE :: read_continents=.FALSE.
25!$OMP THREADPRIVATE(read_continents)
26
27CONTAINS
28!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
29!!
30!! Public subroutines :
31!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
32
33  SUBROUTINE limit_read_frac(itime, dtime, jour, pctsrf_new, is_modified)
34!
35! This subroutine is called from "change_srf_frac" for case of
36! ocean=force or from ocean_slab_frac for ocean=slab.
37! The fraction for all sub-surfaces at actual time step is returned.
38
39    USE dimphy
40    USE indice_sol_mod
41
42! Input arguments
43!****************************************************************************************
44    INTEGER, INTENT(IN) :: itime   ! time step
45    INTEGER, INTENT(IN) :: jour    ! current day
46    REAL   , INTENT(IN) :: dtime   ! length of time step
47 
48! Output arguments
49!****************************************************************************************
50    REAL, DIMENSION(klon,nbsrf), INTENT(OUT) :: pctsrf_new  ! sub surface fractions
51    LOGICAL, INTENT(OUT)                     :: is_modified ! true if pctsrf is modified at this time step
52
53! End declaration
54!****************************************************************************************
55
56! 1) Read file limit.nc
57    CALL limit_read_tot(itime, dtime, jour, is_modified)
58
59! 2) Return the fraction read in limit_read_tot
60    pctsrf_new(:,:) = pctsrf(:,:)
61   
62  END SUBROUTINE limit_read_frac
63
64!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
65
66  SUBROUTINE limit_read_rug_alb(itime, dtime, jour, &
67       knon, knindex, &
68       rugos_out, alb_out)
69!
70! This subroutine is called from surf_land_bucket.
71! The flag "ok_veget" must can not be true. If coupled run, "ocean=couple"
72! then this routine will call limit_read_tot.
73!
74    USE dimphy
75    USE surface_data
76
77! Input arguments
78!****************************************************************************************
79    INTEGER, INTENT(IN) :: itime                     ! numero du pas de temps courant
80    INTEGER, INTENT(IN) :: jour                      ! jour a lire dans l'annee
81    REAL   , INTENT(IN) :: dtime                     ! pas de temps de la physique (en s)
82    INTEGER, INTENT(IN) :: knon                      ! nomber of points on compressed grid
83    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex  ! grid point number for compressed grid
84! Output arguments
85!****************************************************************************************
86    REAL, DIMENSION(klon), INTENT(OUT) :: rugos_out
87    REAL, DIMENSION(klon), INTENT(OUT) :: alb_out
88   
89! Local variables
90!****************************************************************************************
91    INTEGER :: i
92    LOGICAL :: is_modified
93!****************************************************************************************
94
95IF (type_ocean == 'couple'.OR. &
96         (type_ocean == 'slab' .AND. version_ocean == 'sicINT')) THEN
97       ! limit.nc has not yet been read. Do it now!
98       CALL limit_read_tot(itime, dtime, jour, is_modified)
99    END IF
100
101    DO i=1,knon
102       rugos_out(i) = rugos(knindex(i))
103       alb_out(i)  = albedo(knindex(i))
104    END DO
105
106  END SUBROUTINE limit_read_rug_alb
107
108!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
109
110  SUBROUTINE limit_read_sst(knon, knindex, sst_out)
111!
112! This subroutine returns the sea surface temperature already read from limit.nc.
113!
114    USE dimphy, ONLY : klon
115
116    INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
117    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex  ! grid point number for compressed grid
118    REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out
119
120    INTEGER :: i
121
122    DO i = 1, knon
123       sst_out(i) = sst(knindex(i))
124    END DO
125
126  END SUBROUTINE limit_read_sst
127
128!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
129!!
130!! Private subroutine :
131!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
132
133  SUBROUTINE limit_read_tot(itime, dtime, jour, is_modified)
134!
135! Read everything needed from limit.nc
136!
137! 0) Initialize
138! 1) Open the file limit.nc, if it is time
139! 2) Read fraction, if not type_ocean=couple
140! 3) Read sea surface temperature, if not type_ocean=couple
141! 4) Read albedo and rugosity for land surface, only in case of no vegetation model
142! 5) Close file and distribuate variables to all processus
143
144    USE dimphy
145    USE mod_grid_phy_lmdz
146    USE mod_phys_lmdz_para
147    USE surface_data, ONLY : type_ocean, ok_veget
148    USE netcdf
149    USE indice_sol_mod
150    USE phys_cal_mod, ONLY : calend, year_len
151    USE print_control_mod, ONLY: lunout, prt_level
152
153    IMPLICIT NONE
154   
155! In- and ouput arguments
156!****************************************************************************************
157    INTEGER, INTENT(IN) :: itime   ! numero du pas de temps courant
158    INTEGER, INTENT(IN) :: jour    ! jour a lire dans l'annee
159    REAL   , INTENT(IN) :: dtime   ! pas de temps de la physique (en s)
160
161    LOGICAL, INTENT(OUT) :: is_modified  ! true if pctsrf is modified at this time step
162
163! Locals variables with attribute SAVE
164!****************************************************************************************
165! frequence de lecture des conditions limites (en pas de physique)
166    INTEGER,SAVE                              :: lmt_pas
167!$OMP THREADPRIVATE(lmt_pas)
168    LOGICAL, SAVE                             :: first_call=.TRUE.
169!$OMP THREADPRIVATE(first_call) 
170    INTEGER, SAVE                             :: jour_lu = -1
171!$OMP THREADPRIVATE(jour_lu) 
172! Locals variables
173!****************************************************************************************
174    INTEGER                                   :: nid, nvarid, ndimid, nn
175    INTEGER                                   :: ii, ierr
176    INTEGER, DIMENSION(2)                     :: start, epais
177    REAL, DIMENSION(klon_glo,nbsrf)           :: pct_glo  ! fraction at global grid
178    REAL, DIMENSION(klon_glo)                 :: sst_glo  ! sea-surface temperature at global grid
179    REAL, DIMENSION(klon_glo)                 :: rug_glo  ! rugosity at global grid
180    REAL, DIMENSION(klon_glo)                 :: alb_glo  ! albedo at global grid
181    CHARACTER(len=20)                         :: modname='limit_read_mod'     
182    CHARACTER(LEN=99)                         :: abort_message, calendar, str
183
184! End declaration
185!****************************************************************************************
186
187!****************************************************************************************
188! 0) Initialization
189!
190!****************************************************************************************
191    IF (first_call) THEN
192       ! calculate number of time steps for one day
193       lmt_pas = NINT(86400./dtime * 1.0)
194       
195       ! Allocate module save variables
196       IF ( type_ocean /= 'couple' ) THEN
197          ALLOCATE(pctsrf(klon,nbsrf), sst(klon), stat=ierr)
198          IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating pctsrf and sst',1)
199       END IF
200
201       IF ( .NOT. ok_veget ) THEN
202          ALLOCATE(rugos(klon), albedo(klon), stat=ierr)
203          IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating rugos and albedo',1)
204       END IF
205
206!$OMP MASTER  ! Only master thread
207       IF (is_mpi_root) THEN ! Only master processus
208          ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
209          IF (ierr /= NF90_NOERR) CALL abort_physic(modname,&
210               'Pb d''ouverture du fichier de conditions aux limites',1)
211
212          !--- WARNING IF CALENDAR IS KNOWN AND DOES NOT MATCH THE ONE OF LMDZ
213          ierr=NF90_INQ_VARID(nid, 'TEMPS', nvarid)
214          ierr=NF90_GET_ATT(nid, nvarid, 'calendar', calendar)
215          IF(ierr==NF90_NOERR.AND.calendar/=calend.AND.prt_level>=1) THEN
216             WRITE(lunout,*)'BEWARE: gcm and limit.nc calendars differ: '
217             WRITE(lunout,*)'  '//TRIM(calend)//' for gcm'
218             WRITE(lunout,*)'  '//TRIM(calendar)//' for limit.nc file'
219          END IF
220
221          !--- ERROR IF FILE RECORDS NUMBER IS NOT EQUAL TO EXPECTED NUMBER OF DAYS
222          ierr=NF90_INQUIRE(nid, UnlimitedDimID=ndimid)
223          ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
224          WRITE(abort_message,'(a,2(i0,a))')'limit.nc records number (',nn,') does no'//&
225            't match year length (',year_len,')'
226          IF(nn/=year_len) CALL abort_physic(modname,abort_message,1)
227
228          !--- ERROR IF FILES AND LMDZ HORIZONTAL RESOLUTIONS DO NOT MATCH
229          ierr=NF90_INQ_DIMID(nid, 'points_physiques', ndimid)
230          ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
231          WRITE(abort_message,'(a,2(i0,a))')'limit.nc horizontal number of cells (',nn, &
232            ') does not match LMDZ klon_glo (',klon_glo,')'
233          IF(nn/=klon_glo) CALL abort_physic(modname,abort_message,1)
234
235          ierr = NF90_CLOSE(nid)
236          IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Pb when closing file', 1)
237          first_call=.FALSE.
238       END IF ! is_mpi_root
239!$OMP END MASTER
240!$OMP BARRIER
241    END IF
242
243!****************************************************************************************
244! 1) Open the file limit.nc if it is the right moment to read, once a day.
245!    The file is read only by the master thread of the master mpi process(is_mpi_root)
246!    Check by the way if the number of records is correct.
247!
248!****************************************************************************************
249
250    is_modified = .FALSE.
251    IF (MOD(itime-1, lmt_pas) == 0 .OR. jour_lu /= jour ) THEN   ! time to read
252       jour_lu = jour
253       is_modified = .TRUE.
254!$OMP MASTER  ! Only master thread
255       IF (is_mpi_root) THEN ! Only master processus
256
257          ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
258          IF (ierr /= NF90_NOERR) CALL abort_physic(modname,&
259               'Pb d''ouverture du fichier de conditions aux limites',1)
260
261          ! La tranche de donnees a lire:
262          start(1) = 1
263          start(2) = jour
264          epais(1) = klon_glo
265          epais(2) = 1
266
267
268!****************************************************************************************
269! 2) Read fraction if not type_ocean=couple
270!
271!****************************************************************************************
272
273          IF ( type_ocean /= 'couple') THEN
274!
275! Ocean fraction
276             ierr = NF90_INQ_VARID(nid, 'FOCE', nvarid)
277             IF (ierr /= NF90_NOERR) CALL abort_physic(modname, 'Le champ <FOCE> est absent',1)
278             
279             ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_oce),start,epais)
280             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FOCE>' ,1)
281!
282! Sea-ice fraction
283             ierr = NF90_INQ_VARID(nid, 'FSIC', nvarid)
284             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FSIC> est absent',1)
285
286             ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_sic),start,epais)
287             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FSIC>' ,1)
288
289
290! Read land and continentals fraction only if asked for
291             IF (read_continents .OR. itime == 1) THEN
292!
293! Land fraction
294                ierr = NF90_INQ_VARID(nid, 'FTER', nvarid)
295                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FTER> est absent',1)
296               
297                ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_ter),start,epais)
298                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FTER>',1)
299!
300! Continentale ice fraction
301                ierr = NF90_INQ_VARID(nid, 'FLIC', nvarid)
302                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FLIC> est absent',1)
303
304                ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_lic),start,epais)
305                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FLIC>',1)
306             END IF
307
308          END IF ! type_ocean /= couple
309
310!****************************************************************************************
311! 3) Read sea-surface temperature, if not coupled ocean
312!
313!****************************************************************************************
314          IF ( type_ocean /= 'couple') THEN
315
316             ierr = NF90_INQ_VARID(nid, 'SST', nvarid)
317             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <SST> est absent',1)
318
319             ierr = NF90_GET_VAR(nid,nvarid,sst_glo,start,epais)
320             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <SST>',1)
321         
322          END IF
323
324!****************************************************************************************
325! 4) Read albedo and rugosity for land surface, only in case of no vegetation model
326!
327!****************************************************************************************
328
329          IF (.NOT. ok_veget) THEN
330!
331! Read albedo
332             ierr = NF90_INQ_VARID(nid, 'ALB', nvarid)
333             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <ALB> est absent',1)
334
335             ierr = NF90_GET_VAR(nid,nvarid,alb_glo,start,epais)
336             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <ALB>',1)
337!
338! Read rugosity
339             ierr = NF90_INQ_VARID(nid, 'RUG', nvarid)
340             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <RUG> est absent',1)
341
342             ierr = NF90_GET_VAR(nid,nvarid,rug_glo,start,epais)
343             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <RUG>',1)
344
345          END IF
346
347!****************************************************************************************
348! 5) Close file and distribuate variables to all processus
349!
350!****************************************************************************************
351          ierr = NF90_CLOSE(nid)
352          IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Pb when closing file', 1)
353       ENDIF ! is_mpi_root
354
355!$OMP END MASTER
356!$OMP BARRIER
357
358       IF ( type_ocean /= 'couple') THEN
359          CALL Scatter(sst_glo,sst)
360          CALL Scatter(pct_glo(:,is_oce),pctsrf(:,is_oce))
361          CALL Scatter(pct_glo(:,is_sic),pctsrf(:,is_sic))
362          IF (read_continents .OR. itime == 1) THEN
363             CALL Scatter(pct_glo(:,is_ter),pctsrf(:,is_ter))
364             CALL Scatter(pct_glo(:,is_lic),pctsrf(:,is_lic))
365          END IF
366       END IF
367
368       IF (.NOT. ok_veget) THEN
369          CALL Scatter(alb_glo, albedo)
370          CALL Scatter(rug_glo, rugos)
371       END IF
372
373    ENDIF ! time to read
374
375  END SUBROUTINE limit_read_tot
376
377END MODULE limit_read_mod
Note: See TracBrowser for help on using the repository browser.