source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/limit_read_mod.F90 @ 3983

Last change on this file since 3983 was 3895, checked in by ymipsl, 9 years ago

Make LMDZ5 be compliant to generate initiale state and compute in openMP mode suing unstructured grid.

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