GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/limit_read_mod.F90 Lines: 107 115 93.0 %
Date: 2023-06-30 12:51:15 Branches: 75 138 54.3 %

Line Branch Exec Source
1
!
2
! $Id: limit_read_mod.F90 3435 2019-01-22 15:21:59Z fairhead $
3
!
4
MODULE 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
27
CONTAINS
28
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
29
!!
30
!! Public subroutines :
31
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
32
33
34
1
  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
1
  END SUBROUTINE init_limit_read
54
55
288
  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
288
    CALL limit_read_tot(itime, dtime, jour, is_modified)
80
81
! 2) Return the fraction read in limit_read_tot
82

1146528
    pctsrf_new(:,:) = pctsrf(:,:)
83
84
288
  END SUBROUTINE limit_read_frac
85
86
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
87
88
288
  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
117

288
IF (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
148896
    DO i=1,knon
124
148608
       rugos_out(i) = rugos(knindex(i))
125
148896
       alb_out(i)  = albedo(knindex(i))
126
    END DO
127
128
288
  END SUBROUTINE limit_read_rug_alb
129
130
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
131
132
288
  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
218016
    DO i = 1, knon
145
218016
       sst_out(i) = sst(knindex(i))
146
    END DO
147
148
288
  END SUBROUTINE limit_read_sst
149
150
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
151
!!
152
!! Private subroutine :
153
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
154
155
288
  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
576
    REAL, DIMENSION(klon_glo,nbsrf)           :: pct_glo  ! fraction at global grid
203
576
    REAL, DIMENSION(klon_glo)                 :: sst_glo  ! sea-surface temperature at global grid
204
576
    REAL, DIMENSION(klon_glo)                 :: rug_glo  ! rugosity at global grid
205
576
    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
288
    IF (first_call) THEN
223
1
       first_call=.FALSE.
224
       ! calculate number of time steps for one day
225
1
       lmt_pas = NINT(86400./dtime * 1.0)
226
227
       ! Allocate module save variables
228
1
       IF ( type_ocean /= 'couple' ) THEN
229




1
          ALLOCATE(pctsrf(klon,nbsrf), sst(klon), stat=ierr)
230
1
          IF (ierr /= 0) CALL abort_physic(modname, 'PB in allocating pctsrf and sst',1)
231
       END IF
232
233
1
       IF ( .NOT. ok_veget ) THEN
234




1
          ALLOCATE(rugos(klon), albedo(klon), stat=ierr)
235
1
          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
1
       IF (is_mpi_root) THEN ! Only master processus
240
1
          ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
241
1
          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
1
          ierr=NF90_INQ_VARID(nid, 'TEMPS', nvarid)
246
1
          ierr=NF90_GET_ATT(nid, nvarid, 'calendar', calendar)
247

1
          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
1
          IF (grid_type==unstructured) THEN
255
            ierr=NF90_INQ_DIMID(nid,"time_year",ndimid)
256
          ELSE
257
1
            ierr=NF90_INQUIRE(nid, UnlimitedDimID=ndimid)
258
          ENDIF
259
1
          ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
260
1
          WRITE(abort_message,'(a,2(i0,a))')'limit.nc records number (',nn,') does no'//&
261
2
            't match year length (',year_len,')'
262
1
          IF(nn/=year_len) CALL abort_physic(modname,abort_message,1)
263
264
          !--- ERROR IF FILES AND LMDZ HORIZONTAL RESOLUTIONS DO NOT MATCH
265
1
          IF (grid_type==unstructured) THEN
266
            ierr=NF90_INQ_DIMID(nid, 'cell', ndimid)
267
          ELSE
268
1
            ierr=NF90_INQ_DIMID(nid, 'points_physiques', ndimid)
269
          ENDIF
270
1
          ierr=NF90_INQUIRE_DIMENSION(nid, ndimid, len=nn)
271
1
          WRITE(abort_message,'(a,2(i0,a))')'limit.nc horizontal number of cells (',nn, &
272
2
            ') does not match LMDZ klon_glo (',klon_glo,')'
273
1
          IF(nn/=klon_glo) CALL abort_physic(modname,abort_message,1)
274
275
1
          ierr = NF90_CLOSE(nid)
276
1
          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
288
    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

288
    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
6
       jour_lu = jour
295
6
       is_modified = .TRUE.
296
297
6
      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
6
       IF (is_mpi_root) THEN ! Only master processus!
340
341
6
          ierr = NF90_OPEN ('limit.nc', NF90_NOWRITE, nid)
342
6
          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
6
          start(1) = 1
347
6
          start(2) = jour
348
6
          epais(1) = klon_glo
349
6
          epais(2) = 1
350
351
352
!****************************************************************************************
353
! 2) Read fraction if not type_ocean=couple
354
!
355
!****************************************************************************************
356
357
6
          IF ( type_ocean /= 'couple') THEN
358
!
359
! Ocean fraction
360
6
             ierr = NF90_INQ_VARID(nid, 'FOCE', nvarid)
361
6
             IF (ierr /= NF90_NOERR) CALL abort_physic(modname, 'Le champ <FOCE> est absent',1)
362
363
6
             ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_oce),start,epais)
364
6
             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FOCE>' ,1)
365
!
366
! Sea-ice fraction
367
6
             ierr = NF90_INQ_VARID(nid, 'FSIC', nvarid)
368
6
             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FSIC> est absent',1)
369
370
6
             ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_sic),start,epais)
371
6
             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

6
             IF (read_continents .OR. itime == 1) THEN
376
!
377
! Land fraction
378
1
                ierr = NF90_INQ_VARID(nid, 'FTER', nvarid)
379
1
                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FTER> est absent',1)
380
381
1
                ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_ter),start,epais)
382
1
                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <FTER>',1)
383
!
384
! Continentale ice fraction
385
1
                ierr = NF90_INQ_VARID(nid, 'FLIC', nvarid)
386
1
                IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <FLIC> est absent',1)
387
388
1
                ierr = NF90_GET_VAR(nid,nvarid,pct_glo(:,is_lic),start,epais)
389
1
                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
6
          IF ( type_ocean /= 'couple') THEN
399
400
6
             ierr = NF90_INQ_VARID(nid, 'SST', nvarid)
401
6
             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <SST> est absent',1)
402
403
6
             ierr = NF90_GET_VAR(nid,nvarid,sst_glo,start,epais)
404
6
             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
6
          IF (.NOT. ok_veget) THEN
414
!
415
! Read albedo
416
6
             ierr = NF90_INQ_VARID(nid, 'ALB', nvarid)
417
6
             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <ALB> est absent',1)
418
419
6
             ierr = NF90_GET_VAR(nid,nvarid,alb_glo,start,epais)
420
6
             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Lecture echouee pour <ALB>',1)
421
!
422
! Read rugosity
423
6
             ierr = NF90_INQ_VARID(nid, 'RUG', nvarid)
424
6
             IF (ierr /= NF90_NOERR) CALL abort_physic(modname,'Le champ <RUG> est absent',1)
425
426
6
             ierr = NF90_GET_VAR(nid,nvarid,rug_glo,start,epais)
427
6
             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
6
          ierr = NF90_CLOSE(nid)
436
6
          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
6
       IF ( type_ocean /= 'couple') THEN
443
6
          CALL Scatter(sst_glo,sst)
444
6
          CALL Scatter(pct_glo(:,is_oce),pctsrf(:,is_oce))
445
6
          CALL Scatter(pct_glo(:,is_sic),pctsrf(:,is_sic))
446

6
          IF (read_continents .OR. itime == 1) THEN
447
1
             CALL Scatter(pct_glo(:,is_ter),pctsrf(:,is_ter))
448
1
             CALL Scatter(pct_glo(:,is_lic),pctsrf(:,is_lic))
449
          END IF
450
       END IF
451
452
6
       IF (.NOT. ok_veget) THEN
453
6
          CALL Scatter(alb_glo, albedo)
454
6
          CALL Scatter(rug_glo, rugos)
455
       END IF
456
457
      ENDIF ! Grid type
458
459
    ENDIF ! time to read
460
461
288
  END SUBROUTINE limit_read_tot
462
463
END MODULE limit_read_mod