source: LMDZ6/trunk/libf/phylmd/limit_read_mod.F90 @ 5130

Last change on this file since 5130 was 5084, checked in by Laurent Fairhead, 4 months ago

Reverting to r4065. Updating fortran standard broke too much stuff. Will do it by smaller chunks
AB, LF
  • 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: 17.8 KB
Line 
1!
2! $Id: limit_read_mod.F90 5084 2024-07-19 16:40:44Z 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 mod_grid_phy_lmdz
36  USE surface_data
37  USE mod_phys_lmdz_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 mod_grid_phy_lmdz
165    USE mod_phys_lmdz_para
166    USE surface_data, ONLY : type_ocean, ok_veget
167    USE netcdf
168    USE indice_sol_mod
169    USE phys_cal_mod, ONLY : calend, year_len
170    USE print_control_mod, ONLY: lunout, prt_level
171    USE lmdz_xios, ONLY: xios_recv_field, using_xios
172   
173    IMPLICIT NONE
174   
175! In- and ouput arguments
176!****************************************************************************************
177    INTEGER, INTENT(IN) :: itime   ! numero du pas de temps courant
178    INTEGER, INTENT(IN) :: jour    ! jour a lire dans l'annee
179    REAL   , INTENT(IN) :: dtime   ! pas de temps de la physique (en s)
180
181    LOGICAL, INTENT(OUT) :: is_modified  ! true if pctsrf is modified at this time step
182
183! Locals variables with attribute SAVE
184!****************************************************************************************
185! frequence de lecture des conditions limites (en pas de physique)
186    INTEGER,SAVE                              :: lmt_pas
187!$OMP THREADPRIVATE(lmt_pas)
188    LOGICAL, SAVE                             :: first_call=.TRUE.
189!$OMP THREADPRIVATE(first_call) 
190    INTEGER, SAVE                             :: jour_lu = -1
191!$OMP THREADPRIVATE(jour_lu) 
192! Locals variables
193!****************************************************************************************
194    INTEGER                                   :: nid, nvarid, ndimid, nn
195    INTEGER                                   :: ii, ierr
196    INTEGER, DIMENSION(2)                     :: start, epais
197    REAL, DIMENSION(klon_glo,nbsrf)           :: pct_glo  ! fraction at global grid
198    REAL, DIMENSION(klon_glo)                 :: sst_glo  ! sea-surface temperature at global grid
199    REAL, DIMENSION(klon_glo)                 :: rug_glo  ! rugosity at global grid
200    REAL, DIMENSION(klon_glo)                 :: alb_glo  ! albedo at global grid
201
202    REAL, DIMENSION(klon_mpi,nbsrf)           :: pct_mpi  ! fraction at global grid
203    REAL, DIMENSION(klon_mpi)                 :: sst_mpi  ! sea-surface temperature at global grid
204    REAL, DIMENSION(klon_mpi)                 :: rug_mpi  ! rugosity at global grid
205    REAL, DIMENSION(klon_mpi)                 :: alb_mpi  ! albedo at global grid
206
207    CHARACTER(len=20)                         :: modname='limit_read_mod'     
208    CHARACTER(LEN=99)                         :: abort_message, calendar, str
209
210! End declaration
211!****************************************************************************************
212
213!****************************************************************************************
214! 0) Initialization
215!
216!****************************************************************************************
217    IF (first_call) THEN
218       first_call=.FALSE.
219       ! calculate number of time steps for one day
220       lmt_pas = NINT(86400./dtime * 1.0)
221       
222       ! Allocate module save variables
223       IF ( type_ocean /= 'couple' ) THEN
224          ALLOCATE(pctsrf(klon,nbsrf), sst(klon), stat=ierr)
225          IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating pctsrf and sst',1)
226       END IF
227
228       IF ( .NOT. ok_veget ) THEN
229          ALLOCATE(rugos(klon), albedo(klon), stat=ierr)
230          IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating rugos and albedo',1)
231       END IF
232
233!$OMP MASTER  ! Only master thread
234       IF (is_mpi_root) THEN ! Only master processus
235          ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
236          IF (ierr /= NF90_NOERR) CALL abort_physic(modname,&
237               'Pb d''ouverture du fichier de conditions aux limites',1)
238
239          !--- WARNING IF CALENDAR IS KNOWN AND DOES NOT MATCH THE ONE OF LMDZ
240          ierr=NF90_INQ_VARID(nid, 'TEMPS', nvarid)
241          ierr=NF90_GET_ATT(nid, nvarid, 'calendar', calendar)
242          IF(ierr==NF90_NOERR.AND.calendar/=calend.AND.prt_level>=1) THEN
243             WRITE(lunout,*)'BEWARE: gcm and limit.nc calendars differ: '
244             WRITE(lunout,*)'  '//TRIM(calend)//' for gcm'
245             WRITE(lunout,*)'  '//TRIM(calendar)//' for limit.nc file'
246          END IF
247
248          !--- ERROR IF FILE RECORDS NUMBER IS NOT EQUAL TO EXPECTED NUMBER OF DAYS         
249          IF (grid_type==unstructured) THEN
250            ierr=NF90_INQ_DIMID(nid,"time_year",ndimid)
251          ELSE
252            ierr=NF90_INQUIRE(nid, UnlimitedDimID=ndimid)
253          ENDIF
254          ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
255          WRITE(abort_message,'(a,2(i0,a))')'limit.nc records number (',nn,') does no'//&
256            't match year length (',year_len,')'
257          IF(nn/=year_len) CALL abort_physic(modname,abort_message,1)
258
259          !--- ERROR IF FILES AND LMDZ HORIZONTAL RESOLUTIONS DO NOT MATCH
260          IF (grid_type==unstructured) THEN
261            ierr=NF90_INQ_DIMID(nid, 'cell', ndimid)
262          ELSE
263            ierr=NF90_INQ_DIMID(nid, 'points_physiques', ndimid)
264          ENDIF
265          ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
266          WRITE(abort_message,'(a,2(i0,a))')'limit.nc horizontal number of cells (',nn, &
267            ') does not match LMDZ klon_glo (',klon_glo,')'
268          IF(nn/=klon_glo) CALL abort_physic(modname,abort_message,1)
269
270          ierr = NF90_CLOSE(nid)
271          IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Pb when closing file', 1)
272       END IF ! is_mpi_root
273!$OMP END MASTER
274!$OMP BARRIER
275    END IF
276
277!****************************************************************************************
278! 1) Open the file limit.nc if it is the right moment to read, once a day.
279!    The file is read only by the master thread of the master mpi process(is_mpi_root)
280!    Check by the way if the number of records is correct.
281!
282!****************************************************************************************
283
284    is_modified = .FALSE.
285!ym    IF (MOD(itime-1, lmt_pas) == 0 .OR. jour_lu /= jour ) THEN   ! time to read
286!  not REALLY PERIODIC
287    IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN   ! time to read
288!    IF (MOD(itime-1, lmt_pas) == 0) THEN   ! time to read
289       jour_lu = jour
290       is_modified = .TRUE.
291
292      IF (grid_type==unstructured) THEN
293
294        IF ( type_ocean /= 'couple') THEN
295
296           IF (is_omp_master) CALL xios_recv_field("foce_limin",pct_mpi(:,is_oce))
297           IF (is_omp_master) CALL xios_recv_field("fsic_limin",pct_mpi(:,is_sic))
298  !         IF (read_continents .OR. itime == 1) THEN
299           IF (is_omp_master) CALL xios_recv_field("fter_limin",pct_mpi(:,is_ter))
300           IF (is_omp_master) CALL xios_recv_field("flic_limin",pct_mpi(:,is_lic))
301  !         ENDIF
302         ENDIF! type_ocean /= couple
303         
304         IF ( type_ocean /= 'couple') THEN                   
305             IF (is_omp_master) CALL xios_recv_field("sst_limin",sst_mpi)
306         ENDIF
307       
308         IF (.NOT. ok_veget) THEN
309           IF (is_omp_master) CALL xios_recv_field("alb_limin",alb_mpi)
310           IF (is_omp_master) CALL xios_recv_field("rug_limin",rug_mpi)
311         ENDIF
312
313       IF ( type_ocean /= 'couple') THEN
314          CALL Scatter_omp(sst_mpi,sst)
315          CALL Scatter_omp(pct_mpi(:,is_oce),pctsrf(:,is_oce))
316          CALL Scatter_omp(pct_mpi(:,is_sic),pctsrf(:,is_sic))
317!          IF (read_continents .OR. itime == 1) THEN
318             CALL Scatter_omp(pct_mpi(:,is_ter),pctsrf(:,is_ter))
319             CALL Scatter_omp(pct_mpi(:,is_lic),pctsrf(:,is_lic))
320!          END IF
321       END IF
322
323       IF (.NOT. ok_veget) THEN
324          CALL Scatter_omp(alb_mpi, albedo)
325          CALL Scatter_omp(rug_mpi, rugos)
326       END IF
327 
328     ELSE      ! grid_type==regular
329
330!$OMP MASTER  ! Only master thread
331       IF (is_mpi_root) THEN ! Only master processus!
332
333          ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
334          IF (ierr /= NF90_NOERR) CALL abort_physic(modname,&
335               'Pb d''ouverture du fichier de conditions aux limites',1)
336
337          ! La tranche de donnees a lire:
338          start(1) = 1
339          start(2) = jour
340          epais(1) = klon_glo
341          epais(2) = 1
342
343
344!****************************************************************************************
345! 2) Read fraction if not type_ocean=couple
346!
347!****************************************************************************************
348
349          IF ( type_ocean /= 'couple') THEN
350!
351! Ocean fraction
352             ierr = NF90_INQ_VARID(nid, 'FOCE', nvarid)
353             IF (ierr /= NF90_NOERR) CALL abort_physic(modname, 'Le champ <FOCE> est absent',1)
354             
355             ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_oce),start,epais)
356             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FOCE>' ,1)
357!
358! Sea-ice fraction
359             ierr = NF90_INQ_VARID(nid, 'FSIC', nvarid)
360             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FSIC> est absent',1)
361
362             ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_sic),start,epais)
363             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FSIC>' ,1)
364
365
366! Read land and continentals fraction only if asked for
367             IF (read_continents .OR. itime == 1) THEN
368!
369! Land fraction
370                ierr = NF90_INQ_VARID(nid, 'FTER', nvarid)
371                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FTER> est absent',1)
372               
373                ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_ter),start,epais)
374                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FTER>',1)
375!
376! Continentale ice fraction
377                ierr = NF90_INQ_VARID(nid, 'FLIC', nvarid)
378                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FLIC> est absent',1)
379
380                ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_lic),start,epais)
381                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FLIC>',1)
382             END IF
383
384          END IF ! type_ocean /= couple
385
386!****************************************************************************************
387! 3) Read sea-surface temperature, if not coupled ocean
388!
389!****************************************************************************************
390          IF ( type_ocean /= 'couple') THEN
391
392             ierr = NF90_INQ_VARID(nid, 'SST', nvarid)
393             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <SST> est absent',1)
394
395             ierr = NF90_GET_VAR(nid,nvarid,sst_glo,start,epais)
396             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <SST>',1)
397         
398          END IF
399
400!****************************************************************************************
401! 4) Read albedo and rugosity for land surface, only in case of no vegetation model
402!
403!****************************************************************************************
404
405          IF (.NOT. ok_veget) THEN
406!
407! Read albedo
408             ierr = NF90_INQ_VARID(nid, 'ALB', nvarid)
409             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <ALB> est absent',1)
410
411             ierr = NF90_GET_VAR(nid,nvarid,alb_glo,start,epais)
412             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <ALB>',1)
413!
414! Read rugosity
415             ierr = NF90_INQ_VARID(nid, 'RUG', nvarid)
416             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <RUG> est absent',1)
417
418             ierr = NF90_GET_VAR(nid,nvarid,rug_glo,start,epais)
419             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <RUG>',1)
420
421          END IF
422
423!****************************************************************************************
424! 5) Close file and distribuate variables to all processus
425!
426!****************************************************************************************
427          ierr = NF90_CLOSE(nid)
428          IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Pb when closing file', 1)
429       ENDIF ! is_mpi_root
430
431!$OMP END MASTER
432!$OMP BARRIER
433
434       IF ( type_ocean /= 'couple') THEN
435          CALL Scatter(sst_glo,sst)
436          CALL Scatter(pct_glo(:,is_oce),pctsrf(:,is_oce))
437          CALL Scatter(pct_glo(:,is_sic),pctsrf(:,is_sic))
438          IF (read_continents .OR. itime == 1) THEN
439             CALL Scatter(pct_glo(:,is_ter),pctsrf(:,is_ter))
440             CALL Scatter(pct_glo(:,is_lic),pctsrf(:,is_lic))
441          END IF
442       END IF
443
444       IF (.NOT. ok_veget) THEN
445          CALL Scatter(alb_glo, albedo)
446          CALL Scatter(rug_glo, rugos)
447       END IF
448
449      ENDIF ! Grid type
450
451    ENDIF ! time to read
452
453  END SUBROUTINE limit_read_tot
454
455END MODULE limit_read_mod
Note: See TracBrowser for help on using the repository browser.