source: LMDZ6/branches/Ocean_skin/libf/phylmd/limit_read_mod.F90 @ 5454

Last change on this file since 5454 was 3605, checked in by lguez, 5 years ago

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