source: LMDZ6/branches/Amaury_dev/libf/phylmd/limit_read_mod.F90 @ 5218

Last change on this file since 5218 was 5123, checked in by abarral, 5 months ago

Correct various minor mistakes from previous commits

  • 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: 18.0 KB
Line 
1
2! $Id: limit_read_mod.F90 5123 2024-07-25 06:45:50Z abarral $
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
34  SUBROUTINE init_limit_read(first_day)
35  USE lmdz_grid_phy
36  USE surface_data
37  USE lmdz_phys_para
38  USE lmdz_xios
39
40  IMPLICIT NONE
41    INTEGER, INTENT(IN) :: first_day
42   
43   
44    IF ( type_ocean /= 'couple') THEN
45      IF (grid_type==unstructured) THEN
46        IF (is_omp_master) CALL xios_set_file_attr("limit_read",enabled=.TRUE.,record_offset=first_day)
47      ENDIF 
48    ENDIF
49
50  END SUBROUTINE init_limit_read
51 
52  SUBROUTINE limit_read_frac(itime, dtime, jour, pctsrf_new, is_modified)
53
54! This SUBROUTINE is called from "change_srf_frac" for case of
55! ocean=force or from ocean_slab_frac for ocean=slab.
56! The fraction for all sub-surfaces at actual time step is returned.
57
58    USE dimphy
59    USE indice_sol_mod
60
61! Input arguments
62!****************************************************************************************
63    INTEGER, INTENT(IN) :: itime   ! time step
64    INTEGER, INTENT(IN) :: jour    ! current day
65    REAL   , INTENT(IN) :: dtime   ! length of time step
66 
67! Output arguments
68!****************************************************************************************
69    REAL, DIMENSION(klon,nbsrf), INTENT(OUT) :: pctsrf_new  ! sub surface fractions
70    LOGICAL, INTENT(OUT)                     :: is_modified ! true if pctsrf is modified at this time step
71
72! End declaration
73!****************************************************************************************
74
75! 1) Read file limit.nc
76    CALL limit_read_tot(itime, dtime, jour, is_modified)
77
78! 2) Return the fraction read in limit_read_tot
79    pctsrf_new(:,:) = pctsrf(:,:)
80   
81  END SUBROUTINE limit_read_frac
82
83!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
84
85  SUBROUTINE limit_read_rug_alb(itime, dtime, jour, &
86       knon, knindex, &
87       rugos_out, alb_out)
88
89! This SUBROUTINE is called from surf_land_bucket.
90! The flag "ok_veget" must can not be true. If coupled run, "ocean=couple"
91! then this routine will CALL limit_read_tot.
92
93    USE dimphy
94    USE surface_data
95
96! Input arguments
97!****************************************************************************************
98    INTEGER, INTENT(IN) :: itime                     ! numero du pas de temps courant
99    INTEGER, INTENT(IN) :: jour                      ! jour a lire dans l'annee
100    REAL   , INTENT(IN) :: dtime                     ! pas de temps de la physique (en s)
101    INTEGER, INTENT(IN) :: knon                      ! nomber of points on compressed grid
102    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex  ! grid point number for compressed grid
103! Output arguments
104!****************************************************************************************
105    REAL, DIMENSION(klon), INTENT(OUT) :: rugos_out
106    REAL, DIMENSION(klon), INTENT(OUT) :: alb_out
107   
108! Local variables
109!****************************************************************************************
110    INTEGER :: i
111    LOGICAL :: is_modified
112!****************************************************************************************
113
114IF (type_ocean == 'couple'.OR. &
115         (type_ocean == 'slab' .AND. version_ocean == 'sicINT')) THEN
116       ! limit.nc has not yet been read. Do it now!
117       CALL limit_read_tot(itime, dtime, jour, is_modified)
118    END IF
119
120    DO i=1,knon
121       rugos_out(i) = rugos(knindex(i))
122       alb_out(i)  = albedo(knindex(i))
123    END DO
124
125  END SUBROUTINE limit_read_rug_alb
126
127!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
128
129  SUBROUTINE limit_read_sst(knon, knindex, sst_out)
130
131! This SUBROUTINE returns the sea surface temperature already read from limit.nc.
132
133    USE dimphy, ONLY: klon
134
135    INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
136    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex  ! grid point number for compressed grid
137    REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out
138
139    INTEGER :: i
140
141    DO i = 1, knon
142       sst_out(i) = sst(knindex(i))
143    END DO
144
145  END SUBROUTINE limit_read_sst
146
147!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
148!!
149!! Private SUBROUTINE :
150!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
151
152  SUBROUTINE limit_read_tot(itime, dtime, jour, is_modified)
153
154! Read everything needed from limit.nc
155
156! 0) Initialize
157! 1) Open the file limit.nc, if it is time
158! 2) Read fraction, if not type_ocean=couple
159! 3) Read sea surface temperature, if not type_ocean=couple
160! 4) Read albedo and rugosity for land surface, ONLY in case of no vegetation model
161! 5) Close file and distribuate variables to all processus
162
163    USE dimphy
164    USE lmdz_grid_phy
165    USE lmdz_phys_para
166    USE surface_data, ONLY: type_ocean, ok_veget
167    USE netcdf, ONLY:nf90_get_var,nf90_inq_varid,nf90_close,nf90_inquire_dimension,&
168            nf90_inquire,nf90_get_att,nf90_inq_dimid,nf90_nowrite,nf90_noerr,nf90_open
169    USE indice_sol_mod
170    USE phys_cal_mod, ONLY: calend, year_len
171    USE lmdz_print_control, ONLY: lunout, prt_level
172    USE lmdz_xios, ONLY: xios_recv_field, using_xios
173    USE lmdz_abort_physic, ONLY: abort_physic
174   
175    IMPLICIT NONE
176   
177! In- and ouput arguments
178!****************************************************************************************
179    INTEGER, INTENT(IN) :: itime   ! numero du pas de temps courant
180    INTEGER, INTENT(IN) :: jour    ! jour a lire dans l'annee
181    REAL   , INTENT(IN) :: dtime   ! pas de temps de la physique (en s)
182
183    LOGICAL, INTENT(OUT) :: is_modified  ! true if pctsrf is modified at this time step
184
185! Locals variables with attribute SAVE
186!****************************************************************************************
187! frequence de lecture des conditions limites (en pas de physique)
188    INTEGER,SAVE                              :: lmt_pas
189!$OMP THREADPRIVATE(lmt_pas)
190    LOGICAL, SAVE                             :: first_call=.TRUE.
191!$OMP THREADPRIVATE(first_call) 
192    INTEGER, SAVE                             :: jour_lu = -1
193!$OMP THREADPRIVATE(jour_lu) 
194! Locals variables
195!****************************************************************************************
196    INTEGER                                   :: nid, nvarid, ndimid, nn
197    INTEGER                                   :: ii, ierr
198    INTEGER, DIMENSION(2)                     :: start, epais
199    REAL, DIMENSION(klon_glo,nbsrf)           :: pct_glo  ! fraction at global grid
200    REAL, DIMENSION(klon_glo)                 :: sst_glo  ! sea-surface temperature at global grid
201    REAL, DIMENSION(klon_glo)                 :: rug_glo  ! rugosity at global grid
202    REAL, DIMENSION(klon_glo)                 :: alb_glo  ! albedo at global grid
203
204    REAL, DIMENSION(klon_mpi,nbsrf)           :: pct_mpi  ! fraction at global grid
205    REAL, DIMENSION(klon_mpi)                 :: sst_mpi  ! sea-surface temperature at global grid
206    REAL, DIMENSION(klon_mpi)                 :: rug_mpi  ! rugosity at global grid
207    REAL, DIMENSION(klon_mpi)                 :: alb_mpi  ! albedo at global grid
208
209    CHARACTER(len=20)                         :: modname='limit_read_mod'     
210    CHARACTER(LEN=99)                         :: abort_message, calendar, str
211
212! End declaration
213!****************************************************************************************
214
215!****************************************************************************************
216! 0) Initialization
217
218!****************************************************************************************
219    IF (first_call) THEN
220       first_call=.FALSE.
221       ! calculate number of time steps for one day
222       lmt_pas = NINT(86400./dtime * 1.0)
223       
224       ! Allocate module save variables
225       IF ( type_ocean /= 'couple' ) THEN
226          ALLOCATE(pctsrf(klon,nbsrf), sst(klon), stat=ierr)
227          IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating pctsrf and sst',1)
228       END IF
229
230       IF ( .NOT. ok_veget ) THEN
231          ALLOCATE(rugos(klon), albedo(klon), stat=ierr)
232          IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating rugos and albedo',1)
233       END IF
234
235!$OMP MASTER  ! Only master thread
236       IF (is_mpi_root) THEN ! Only master processus
237          ierr = nf90_open ('limit.nc', nf90_nowrite, nid)
238          IF (ierr /= nf90_noerr) CALL abort_physic(modname,&
239               'Pb d''ouverture du fichier de conditions aux limites',1)
240
241          !--- WARNING IF CALENDAR IS KNOWN AND DOES NOT MATCH THE ONE OF LMDZ
242          ierr=nf90_inq_varid(nid, 'TEMPS', nvarid)
243          ierr=nf90_get_att(nid, nvarid, 'calendar', calendar)
244          IF(ierr==nf90_noerr.AND.calendar/=calend.AND.prt_level>=1) THEN
245             WRITE(lunout,*)'BEWARE: gcm and limit.nc calendars differ: '
246             WRITE(lunout,*)'  '//TRIM(calend)//' for gcm'
247             WRITE(lunout,*)'  '//TRIM(calendar)//' for limit.nc file'
248          END IF
249
250          !--- ERROR IF FILE RECORDS NUMBER IS NOT EQUAL TO EXPECTED NUMBER OF DAYS         
251          IF (grid_type==unstructured) THEN
252            ierr=nf90_inq_dimid(nid,"time_year",ndimid)
253          ELSE
254            ierr=nf90_inquire(nid, UnlimitedDimID=ndimid)
255          ENDIF
256          ierr=nf90_inquire_dimension(nid, ndimid, len=nn)
257          WRITE(abort_message,'(a,2(i0,a))')'limit.nc records number (',nn,') does no'//&
258            't match year length (',year_len,')'
259          IF(nn/=year_len) CALL abort_physic(modname,abort_message,1)
260
261          !--- ERROR IF FILES AND LMDZ HORIZONTAL RESOLUTIONS DO NOT MATCH
262          IF (grid_type==unstructured) THEN
263            ierr=nf90_inq_dimid(nid, 'cell', ndimid)
264          ELSE
265            ierr=nf90_inq_dimid(nid, 'points_physiques', ndimid)
266          ENDIF
267          ierr=nf90_inquire_dimension(nid, ndimid, len=nn)
268          WRITE(abort_message,'(a,2(i0,a))')'limit.nc horizontal number of cells (',nn, &
269            ') does not match LMDZ klon_glo (',klon_glo,')'
270          IF(nn/=klon_glo) CALL abort_physic(modname,abort_message,1)
271
272          ierr = nf90_close(nid)
273          IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Pb when closing file', 1)
274       END IF ! is_mpi_root
275!$OMP END MASTER
276!$OMP BARRIER
277    END IF
278
279!****************************************************************************************
280! 1) Open the file limit.nc if it is the right moment to read, once a day.
281!    The file is read only by the master thread of the master mpi process(is_mpi_root)
282!    Check by the way if the number of records is correct.
283
284!****************************************************************************************
285
286    is_modified = .FALSE.
287!ym    IF (MOD(itime-1, lmt_pas) == 0 .OR. jour_lu /= jour ) THEN   ! time to read
288!  not REALLY PERIODIC
289    IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN   ! time to read
290!    IF (MOD(itime-1, lmt_pas) == 0) THEN   ! time to read
291       jour_lu = jour
292       is_modified = .TRUE.
293
294      IF (grid_type==unstructured) THEN
295
296        IF ( type_ocean /= 'couple') THEN
297
298           IF (is_omp_master) CALL xios_recv_field("foce_limin",pct_mpi(:,is_oce))
299           IF (is_omp_master) CALL xios_recv_field("fsic_limin",pct_mpi(:,is_sic))
300  !         IF (read_continents .OR. itime == 1) THEN
301           IF (is_omp_master) CALL xios_recv_field("fter_limin",pct_mpi(:,is_ter))
302           IF (is_omp_master) CALL xios_recv_field("flic_limin",pct_mpi(:,is_lic))
303  !         ENDIF
304         ENDIF! type_ocean /= couple
305         
306         IF ( type_ocean /= 'couple') THEN                   
307             IF (is_omp_master) CALL xios_recv_field("sst_limin",sst_mpi)
308         ENDIF
309       
310         IF (.NOT. ok_veget) THEN
311           IF (is_omp_master) CALL xios_recv_field("alb_limin",alb_mpi)
312           IF (is_omp_master) CALL xios_recv_field("rug_limin",rug_mpi)
313         ENDIF
314
315       IF ( type_ocean /= 'couple') THEN
316          CALL Scatter_omp(sst_mpi,sst)
317          CALL Scatter_omp(pct_mpi(:,is_oce),pctsrf(:,is_oce))
318          CALL Scatter_omp(pct_mpi(:,is_sic),pctsrf(:,is_sic))
319!          IF (read_continents .OR. itime == 1) THEN
320             CALL Scatter_omp(pct_mpi(:,is_ter),pctsrf(:,is_ter))
321             CALL Scatter_omp(pct_mpi(:,is_lic),pctsrf(:,is_lic))
322!          END IF
323       END IF
324
325       IF (.NOT. ok_veget) THEN
326          CALL Scatter_omp(alb_mpi, albedo)
327          CALL Scatter_omp(rug_mpi, rugos)
328       END IF
329 
330     ELSE      ! grid_type==regular
331
332!$OMP MASTER  ! Only master thread
333       IF (is_mpi_root) THEN ! Only master processus!
334
335          ierr = nf90_open ('limit.nc', nf90_nowrite, nid)
336          IF (ierr /= nf90_noerr) CALL abort_physic(modname,&
337               'Pb d''ouverture du fichier de conditions aux limites',1)
338
339          ! La tranche de donnees a lire:
340          start(1) = 1
341          start(2) = jour
342          epais(1) = klon_glo
343          epais(2) = 1
344
345
346!****************************************************************************************
347! 2) Read fraction if not type_ocean=couple
348
349!****************************************************************************************
350
351          IF ( type_ocean /= 'couple') THEN
352
353! Ocean fraction
354             ierr = nf90_inq_varid(nid, 'FOCE', nvarid)
355             IF (ierr /= nf90_noerr) CALL abort_physic(modname, 'Le champ <FOCE> est absent',1)
356             
357             ierr = nf90_get_var(nid,nvarid,pct_glo(:,is_oce),start,epais)
358             IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Lecture echouee pour <FOCE>' ,1)
359
360! Sea-ice fraction
361             ierr = nf90_inq_varid(nid, 'FSIC', nvarid)
362             IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Le champ <FSIC> est absent',1)
363
364             ierr = nf90_get_var(nid,nvarid,pct_glo(:,is_sic),start,epais)
365             IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Lecture echouee pour <FSIC>' ,1)
366
367
368! Read land and continentals fraction only if asked for
369             IF (read_continents .OR. itime == 1) THEN
370
371! Land fraction
372                ierr = nf90_inq_varid(nid, 'FTER', nvarid)
373                IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Le champ <FTER> est absent',1)
374               
375                ierr = nf90_get_var(nid,nvarid,pct_glo(:,is_ter),start,epais)
376                IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Lecture echouee pour <FTER>',1)
377
378! Continentale ice fraction
379                ierr = nf90_inq_varid(nid, 'FLIC', nvarid)
380                IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Le champ <FLIC> est absent',1)
381
382                ierr = nf90_get_var(nid,nvarid,pct_glo(:,is_lic),start,epais)
383                IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Lecture echouee pour <FLIC>',1)
384             END IF
385
386          END IF ! type_ocean /= couple
387
388!****************************************************************************************
389! 3) Read sea-surface temperature, if not coupled ocean
390
391!****************************************************************************************
392          IF ( type_ocean /= 'couple') THEN
393
394             ierr = nf90_inq_varid(nid, 'SST', nvarid)
395             IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Le champ <SST> est absent',1)
396
397             ierr = nf90_get_var(nid,nvarid,sst_glo,start,epais)
398             IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Lecture echouee pour <SST>',1)
399         
400          END IF
401
402!****************************************************************************************
403! 4) Read albedo and rugosity for land surface, ONLY in case of no vegetation model
404
405!****************************************************************************************
406
407          IF (.NOT. ok_veget) THEN
408
409! Read albedo
410             ierr = nf90_inq_varid(nid, 'ALB', nvarid)
411             IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Le champ <ALB> est absent',1)
412
413             ierr = nf90_get_var(nid,nvarid,alb_glo,start,epais)
414             IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Lecture echouee pour <ALB>',1)
415
416! Read rugosity
417             ierr = nf90_inq_varid(nid, 'RUG', nvarid)
418             IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Le champ <RUG> est absent',1)
419
420             ierr = nf90_get_var(nid,nvarid,rug_glo,start,epais)
421             IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Lecture echouee pour <RUG>',1)
422
423          END IF
424
425!****************************************************************************************
426! 5) Close file and distribuate variables to all processus
427
428!****************************************************************************************
429          ierr = nf90_close(nid)
430          IF (ierr /= nf90_noerr) CALL abort_physic(modname,'Pb when closing file', 1)
431       ENDIF ! is_mpi_root
432
433!$OMP END MASTER
434!$OMP BARRIER
435
436       IF ( type_ocean /= 'couple') THEN
437          CALL Scatter(sst_glo,sst)
438          CALL Scatter(pct_glo(:,is_oce),pctsrf(:,is_oce))
439          CALL Scatter(pct_glo(:,is_sic),pctsrf(:,is_sic))
440          IF (read_continents .OR. itime == 1) THEN
441             CALL Scatter(pct_glo(:,is_ter),pctsrf(:,is_ter))
442             CALL Scatter(pct_glo(:,is_lic),pctsrf(:,is_lic))
443          END IF
444       END IF
445
446       IF (.NOT. ok_veget) THEN
447          CALL Scatter(alb_glo, albedo)
448          CALL Scatter(rug_glo, rugos)
449       END IF
450
451      ENDIF ! Grid type
452
453    ENDIF ! time to read
454
455  END SUBROUTINE limit_read_tot
456
457END MODULE limit_read_mod
Note: See TracBrowser for help on using the repository browser.