source: LMDZ6/branches/DYNAMICO-conv-GC/libf/phylmd/limit_read_mod.F90 @ 5185

Last change on this file since 5185 was 3406, checked in by jghattas, 6 years ago

Added all modifications in the model code that were used for the simulations with DYANMICO during the Grand Challeng 2018. Modifications done by Y. Meurdesoif, L. Fairhead and A.K. Traore

  • 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 3406 2018-10-23 09:16:41Z 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#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) THEN   ! time to read
293       jour_lu = jour
294       is_modified = .TRUE.
295
296      IF (grid_type==unstructured) THEN
297
298#ifdef CPP_XIOS
299        IF ( type_ocean /= 'couple') THEN
300
301           IF (is_omp_master) CALL xios_recv_field("foce_limin",pct_mpi(:,is_oce))
302           IF (is_omp_master) CALL xios_recv_field("fsic_limin",pct_mpi(:,is_sic))
303  !         IF (read_continents .OR. itime == 1) THEN
304           IF (is_omp_master) CALL xios_recv_field("fter_limin",pct_mpi(:,is_ter))
305           IF (is_omp_master) CALL xios_recv_field("flic_limin",pct_mpi(:,is_lic))
306  !         ENDIF
307         ENDIF! type_ocean /= couple
308         
309         IF ( type_ocean /= 'couple') THEN                   
310             IF (is_omp_master) CALL xios_recv_field("sst_limin",sst_mpi)
311         ENDIF
312       
313         IF (.NOT. ok_veget) THEN
314           IF (is_omp_master) CALL xios_recv_field("alb_limin",alb_mpi)
315           IF (is_omp_master) CALL xios_recv_field("rug_limin",rug_mpi)
316         ENDIF
317
318       IF ( type_ocean /= 'couple') THEN
319          CALL Scatter_omp(sst_mpi,sst)
320          CALL Scatter_omp(pct_mpi(:,is_oce),pctsrf(:,is_oce))
321          CALL Scatter_omp(pct_mpi(:,is_sic),pctsrf(:,is_sic))
322!          IF (read_continents .OR. itime == 1) THEN
323             CALL Scatter_omp(pct_mpi(:,is_ter),pctsrf(:,is_ter))
324             CALL Scatter_omp(pct_mpi(:,is_lic),pctsrf(:,is_lic))
325!          END IF
326       END IF
327
328       IF (.NOT. ok_veget) THEN
329          CALL Scatter_omp(alb_mpi, albedo)
330          CALL Scatter_omp(rug_mpi, rugos)
331       END IF
332#endif
333
334 
335     ELSE      ! grid_type==regular
336
337!$OMP MASTER  ! Only master thread
338       IF (is_mpi_root) THEN ! Only master processus!
339
340          ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
341          IF (ierr /= NF90_NOERR) CALL abort_physic(modname,&
342               'Pb d''ouverture du fichier de conditions aux limites',1)
343
344          ! La tranche de donnees a lire:
345          start(1) = 1
346          start(2) = jour
347          epais(1) = klon_glo
348          epais(2) = 1
349
350
351!****************************************************************************************
352! 2) Read fraction if not type_ocean=couple
353!
354!****************************************************************************************
355
356          IF ( type_ocean /= 'couple') THEN
357!
358! Ocean fraction
359             ierr = NF90_INQ_VARID(nid, 'FOCE', nvarid)
360             IF (ierr /= NF90_NOERR) CALL abort_physic(modname, 'Le champ <FOCE> est absent',1)
361             
362             ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_oce),start,epais)
363             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FOCE>' ,1)
364!
365! Sea-ice fraction
366             ierr = NF90_INQ_VARID(nid, 'FSIC', nvarid)
367             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FSIC> est absent',1)
368
369             ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_sic),start,epais)
370             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FSIC>' ,1)
371
372
373! Read land and continentals fraction only if asked for
374             IF (read_continents .OR. itime == 1) THEN
375!
376! Land fraction
377                ierr = NF90_INQ_VARID(nid, 'FTER', nvarid)
378                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FTER> est absent',1)
379               
380                ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_ter),start,epais)
381                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FTER>',1)
382!
383! Continentale ice fraction
384                ierr = NF90_INQ_VARID(nid, 'FLIC', nvarid)
385                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FLIC> est absent',1)
386
387                ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_lic),start,epais)
388                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FLIC>',1)
389             END IF
390
391          END IF ! type_ocean /= couple
392
393!****************************************************************************************
394! 3) Read sea-surface temperature, if not coupled ocean
395!
396!****************************************************************************************
397          IF ( type_ocean /= 'couple') THEN
398
399             ierr = NF90_INQ_VARID(nid, 'SST', nvarid)
400             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <SST> est absent',1)
401
402             ierr = NF90_GET_VAR(nid,nvarid,sst_glo,start,epais)
403             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <SST>',1)
404         
405          END IF
406
407!****************************************************************************************
408! 4) Read albedo and rugosity for land surface, only in case of no vegetation model
409!
410!****************************************************************************************
411
412          IF (.NOT. ok_veget) THEN
413!
414! Read albedo
415             ierr = NF90_INQ_VARID(nid, 'ALB', nvarid)
416             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <ALB> est absent',1)
417
418             ierr = NF90_GET_VAR(nid,nvarid,alb_glo,start,epais)
419             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <ALB>',1)
420!
421! Read rugosity
422             ierr = NF90_INQ_VARID(nid, 'RUG', nvarid)
423             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <RUG> est absent',1)
424
425             ierr = NF90_GET_VAR(nid,nvarid,rug_glo,start,epais)
426             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <RUG>',1)
427
428          END IF
429
430!****************************************************************************************
431! 5) Close file and distribuate variables to all processus
432!
433!****************************************************************************************
434          ierr = NF90_CLOSE(nid)
435          IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Pb when closing file', 1)
436       ENDIF ! is_mpi_root
437
438!$OMP END MASTER
439!$OMP BARRIER
440
441       IF ( type_ocean /= 'couple') THEN
442          CALL Scatter(sst_glo,sst)
443          CALL Scatter(pct_glo(:,is_oce),pctsrf(:,is_oce))
444          CALL Scatter(pct_glo(:,is_sic),pctsrf(:,is_sic))
445          IF (read_continents .OR. itime == 1) THEN
446             CALL Scatter(pct_glo(:,is_ter),pctsrf(:,is_ter))
447             CALL Scatter(pct_glo(:,is_lic),pctsrf(:,is_lic))
448          END IF
449       END IF
450
451       IF (.NOT. ok_veget) THEN
452          CALL Scatter(alb_glo, albedo)
453          CALL Scatter(rug_glo, rugos)
454       END IF
455
456      ENDIF ! Grid type
457
458    ENDIF ! time to read
459
460  END SUBROUTINE limit_read_tot
461
462END MODULE limit_read_mod
Note: See TracBrowser for help on using the repository browser.