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

Last change on this file since 2666 was 2666, checked in by fhourdin, 8 years ago

Correction pour compilation gfortran
  • 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       first_call=.FALSE.
207    ENDIF
208 
209!****************************************************************************************
210! 1) Open the file limit.nc if it is the right moment to read, once a day.
211!    The file is read only by the master thread of the master mpi process(is_mpi_root)
212!    Check by the way if the number of records is correct.
213!
214!****************************************************************************************
215
216    is_modified = .FALSE.
217    IF (MOD(itime-1, lmt_pas) == 0 .OR. jour_lu /= jour ) THEN   ! time to read
218       jour_lu = jour
219       is_modified = .TRUE.
220!$OMP MASTER  ! Only master thread
221       IF (is_mpi_root) THEN ! Only master processus
222
223          ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
224          IF (ierr /= NF90_NOERR) CALL abort_physic(modname,&
225               'Pb d''ouverture du fichier de conditions aux limites',1)
226
227          !--- WARNING IF CALENDAR IS KNOWN AND DOES NOT MATCH THE ONE OF LMDZ
228          ierr=NF90_INQ_VARID(nid, 'TEMPS', nvarid)
229          ierr=NF90_GET_ATT(nid, nvarid, 'calendar', calendar)
230          IF(ierr==NF90_NOERR.AND.calendar/=calend.AND.prt_level>=1) THEN
231             WRITE(lunout,*)'BEWARE: gcm and limit.nc calendars differ: '
232             WRITE(lunout,*)'  '//TRIM(calend)//' for gcm'
233             WRITE(lunout,*)'  '//TRIM(calendar)//' for limit.nc file'
234          END IF
235
236          !--- ERROR IF FILE RECORDS NUMBER IS NOT EQUAL TO EXPECTED NUMBER OF DAYS
237          ierr=NF90_INQUIRE(nid, UnlimitedDimID=ndimid)
238          ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
239          str=""; DO WHILE(nn>0); str=TRIM(str)//CHAR(nn-10*(nn/10)-48); nn=nn/10; END DO
240          abort_message='limit.nc records number ('//TRIM(str)//') does'//&
241            ' not match year length ('
242          nn=year_len
243          str=""; DO WHILE(nn>0); str=TRIM(str)//CHAR(nn-10*(nn/10)-48); nn=nn/10; END DO
244          abort_message=TRIM(abort_message)//TRIM(str)//')'
245          IF(nn/=year_len) CALL abort_physic(modname,abort_message,1)
246
247          !--- ERROR IF FILES AND LMDZ HORIZONTAL RESOLUTIONS DO NOT MATCH
248          ierr=NF90_INQ_DIMID(nid, 'points_physiques', ndimid)
249          ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
250          str=""; DO WHILE(nn>0); str=TRIM(str)//CHAR(nn-10*(nn/10)-48); nn=nn/10; END DO
251          abort_message='limit.nc horizontal number of cells ('//TRIM(str)//') does'//&
252            ' not match LMDZ klon_glo ('
253          nn=klon_glo
254          str=""; DO WHILE(nn>0); str=TRIM(str)//CHAR(nn-10*(nn/10)-48); nn=nn/10; END DO
255          abort_message=TRIM(abort_message)//TRIM(str)//')'
256          IF(nn/=klon_glo) CALL abort_physic(modname,abort_message,1)
257
258          ! La tranche de donnees a lire:
259          start(1) = 1
260          start(2) = jour
261          epais(1) = klon_glo
262          epais(2) = 1
263
264
265!****************************************************************************************
266! 2) Read fraction if not type_ocean=couple
267!
268!****************************************************************************************
269
270          IF ( type_ocean /= 'couple') THEN
271!
272! Ocean fraction
273             ierr = NF90_INQ_VARID(nid, 'FOCE', nvarid)
274             IF (ierr /= NF90_NOERR) CALL abort_physic(modname, 'Le champ <FOCE> est absent',1)
275             
276             ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_oce),start,epais)
277             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FOCE>' ,1)
278!
279! Sea-ice fraction
280             ierr = NF90_INQ_VARID(nid, 'FSIC', nvarid)
281             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FSIC> est absent',1)
282
283             ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_sic),start,epais)
284             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FSIC>' ,1)
285
286
287! Read land and continentals fraction only if asked for
288             IF (read_continents .OR. itime == 1) THEN
289!
290! Land fraction
291                ierr = NF90_INQ_VARID(nid, 'FTER', nvarid)
292                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FTER> est absent',1)
293               
294                ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_ter),start,epais)
295                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FTER>',1)
296!
297! Continentale ice fraction
298                ierr = NF90_INQ_VARID(nid, 'FLIC', nvarid)
299                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FLIC> est absent',1)
300
301                ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_lic),start,epais)
302                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FLIC>',1)
303             END IF
304
305          END IF ! type_ocean /= couple
306
307!****************************************************************************************
308! 3) Read sea-surface temperature, if not coupled ocean
309!
310!****************************************************************************************
311          IF ( type_ocean /= 'couple') THEN
312
313             ierr = NF90_INQ_VARID(nid, 'SST', nvarid)
314             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <SST> est absent',1)
315
316             ierr = NF90_GET_VAR(nid,nvarid,sst_glo,start,epais)
317             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <SST>',1)
318         
319          END IF
320
321!****************************************************************************************
322! 4) Read albedo and rugosity for land surface, only in case of no vegetation model
323!
324!****************************************************************************************
325
326          IF (.NOT. ok_veget) THEN
327!
328! Read albedo
329             ierr = NF90_INQ_VARID(nid, 'ALB', nvarid)
330             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <ALB> est absent',1)
331
332             ierr = NF90_GET_VAR(nid,nvarid,alb_glo,start,epais)
333             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <ALB>',1)
334!
335! Read rugosity
336             ierr = NF90_INQ_VARID(nid, 'RUG', nvarid)
337             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <RUG> est absent',1)
338
339             ierr = NF90_GET_VAR(nid,nvarid,rug_glo,start,epais)
340             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <RUG>',1)
341
342          END IF
343
344!****************************************************************************************
345! 5) Close file and distribuate variables to all processus
346!
347!****************************************************************************************
348          ierr = NF90_CLOSE(nid)
349          IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Pb when closing file', 1)
350       ENDIF ! is_mpi_root
351
352!$OMP END MASTER
353!$OMP BARRIER
354
355       IF ( type_ocean /= 'couple') THEN
356          CALL Scatter(sst_glo,sst)
357          CALL Scatter(pct_glo(:,is_oce),pctsrf(:,is_oce))
358          CALL Scatter(pct_glo(:,is_sic),pctsrf(:,is_sic))
359          IF (read_continents .OR. itime == 1) THEN
360             CALL Scatter(pct_glo(:,is_ter),pctsrf(:,is_ter))
361             CALL Scatter(pct_glo(:,is_lic),pctsrf(:,is_lic))
362          END IF
363       END IF
364
365       IF (.NOT. ok_veget) THEN
366          CALL Scatter(alb_glo, albedo)
367          CALL Scatter(rug_glo, rugos)
368       END IF
369
370    ENDIF ! time to read
371
372  END SUBROUTINE limit_read_tot
373
374
375END MODULE limit_read_mod
Note: See TracBrowser for help on using the repository browser.