source: branches/publications/ORCHIDEE_gmd-2018-261/src_driver/getlandseamask.f90 @ 7844

Last change on this file since 7844 was 4333, checked in by josefine.ghattas, 8 years ago
  • Merge with trunk for revision r3623 - r3977
  • Some corrections for gfortran
File size: 38.1 KB
Line 
1MODULE getlandseamask
2  !
3  !================================================================================================================================
4  !! MODULE   : getlandseamask
5  !!
6  !> BRIEF      This module contains a few methods in order to build a land surface mask which can serve
7  !!            to test the routing and its ability to build a river network.
8  !!
9  !! DESCRIPTION  : None
10  !!
11  !! RECENT CHANGE(S): None
12  !!
13  !! MAIN OUTPUT VARIABLE(S):
14  !!
15  !! REFERENCES   : None
16  !!
17  !! FLOWCHART    : This module has 2 possible modes for building the land/sea mask :
18  !!                TOPOFILE : This configuration option contains a file name of the orography file
19  !!                           which will be used in order to build the land/sea mask.
20  !!                IIM_LON, JJM_lon : the number of points the grid should have in longitude and
21  !!                                   latitude. This is only used when we build the mask based on orography.
22  !!                                   Should these parameters not be provided, then the gird of the orography
23  !!                                   file is used.
24  !!                CONTFRACFILE :: File name in which we will find variable "Contfrac". This land/sea mask will be
25  !!                                used with the associated lat/lon grid.
26  !!                WEST_EAST,  SOUTH_NORTH : The zoom to be performed on the grid so as to retain only this window
27  !!                                          of the land-sea mask.
28  !! \n
29  !_================================================================================================================================
30
31  USE ioipsl_para
32  USE mod_orchidee_para
33  USE control
34  USE constantes_soil_var
35  USE constantes_var
36  USE constantes
37
38  USE netcdf
39
40  IMPLICIT NONE
41
42  PRIVATE
43  PUBLIC :: getlandseamask_init, getlandseamask_read
44
45  INTEGER(i_std), SAVE                           :: iim, jjm, nbland
46  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: lon, lat
47  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: contfrac, orog
48
49CONTAINS
50  !================================================================================================================================
51  !! SUBROUTINE   : getlandseamask_init
52  !!
53  !> BRIEF         Builds the land/sea mask and returns its size to the calling program.
54  !!
55  !! DESCRIPTION  : None
56  !!
57  !! RECENT CHANGE(S): None
58  !!
59  !! MAIN OUTPUT VARIABLE(S):
60  !!
61  !! REFERENCES   : None
62  !!
63  !! FLOWCHART    : None
64  !! \n
65  !_================================================================================================================================
66  SUBROUTINE getlandseamask_init(iim_out, jjm_out, nbland_out)
67    !
68    ! ARGUMENTS
69    !
70    INTEGER(i_std), INTENT(out) :: iim_out, jjm_out, nbland_out
71    !
72    ! LOCAL
73    !
74    CHARACTER(LEN=120)                                  :: filename                !! File name for topography (var: ROSE in "etopo20.nc")
75    CHARACTER(LEN=120)                                  :: contfracfile            !! filename for contfrac
76    INTEGER(i_std)                                      :: iim_full, jjm_full      !! Lon/Lat dimensions for ROSE in etopo20.nc
77    REAL(r_std), ALLOCATABLE, DIMENSION (:,:)           :: lon_full, lat_full      !! Lon/Lat array for ROSE
78    REAL(r_std), ALLOCATABLE, DIMENSION (:,:)           :: etopo, contfrac_full    !! Array for reading ROSE in etopo20.nc
79    INTEGER(i_std), DIMENSION(1)                        :: ibeg, iend
80    INTEGER(i_std), DIMENSION(1)                        :: jbeg, jend
81    INTEGER(i_std)                                      :: fid, i, j
82    REAL(r_std)                                         :: dx, dy
83    !
84    REAL(r_std), DIMENSION(2)                           :: zoom_lon, zoom_lat               !! For crop region when reading
85    !
86    LOGICAL                                             :: buildcontfrac=.FALSE.
87    !
88    !  Read the full topography which we will use to get the land-see mask for all resolutions.
89    !
90    !Config Key   = TOPOFILE
91    !Config Desc  = File containing high resolution topographic informations
92    !Config If    =
93    !Config Def   = etopo20.nc
94    !Config Help  = This needs to be a netCDF file.
95    !Config Units = [-]
96    !-
97    filename = 'NONE'
98    CALL getin('TOPOFILE', filename)
99    WRITE(*,*) "TOPOFILE = ", filename
100    !
101    !Config Key   = CONTFRACFILE
102    !Config Desc  = File containing a contfrac variable as defined in ORCHIDEE.
103    !Config If    =
104    !Config Def   = NONE
105    !Config Help  = This needs to be a netCDF file and can be any history file ORCHIDEE
106    !               has produced ina previous simulation.
107    !Config Units = [-]
108    !-
109    contfracfile = 'NONE'
110    CALL getin('CONTFRACFILE', contfracfile)
111    WRITE(*,*) "CONTFRACFILE = ", contfracfile
112    !
113    ! Initialize Orchidee parameter and calendar
114    !
115    IF ( INDEX(filename,"NONE") <= 0 ) THEN
116       IF ( is_root_prc) THEN
117          CALL topo_getsize(filename, iim_full, jjm_full)
118          !
119          ALLOCATE(lon_full(iim_full,jjm_full),lat_full(iim_full,jjm_full))
120          ALLOCATE(etopo(iim_full,jjm_full))
121          !
122          CALL topo_getvar(filename, iim_full, jjm_full, lon_full, lat_full, etopo)
123          !
124          WRITE(*,*) MINVAL(lon_full), ' < LON_FULL < ', MAXVAL(lon_full)
125          WRITE(*,*) MINVAL(lat_full), ' < LAT_FULL < ', MAXVAL(lat_full)
126          WRITE(*,*) 'ETOPO :', MINVAL(etopo), MAXVAL(etopo)
127          !
128          buildcontfrac=.TRUE.
129          !
130       ENDIF
131    ELSE IF ( INDEX( contfracfile,"NONE") <= 0 ) THEN
132       IF ( is_root_prc) THEN
133          !
134          ! Read lon, lat and contfrac from a netCDF file generated by ORCHIDEE
135          !
136          CALL opencontfrac(contfracfile, iim_full, jjm_full, fid)
137          !
138          ALLOCATE(lon_full(iim_full,jjm_full), lat_full(iim_full,jjm_full), contfrac_full(iim_full,jjm_full))
139          CALL readcontfrac(fid, contfracfile, iim_full, jjm_full, lon_full, lat_full, contfrac_full)
140          !
141          buildcontfrac=.FALSE.
142          !
143       ENDIF
144    ELSE
145       WRITE(*,*) "Neither a orography nore a contfrac file was provided"
146       WRITE(*,*) "CONTFRACFILE = ", contfracfile
147       write(*,*) "TOPOFILE = ", filename
148    ENDIF
149    !
150    !-
151    !- Define the zoom
152    !-
153    zoom_lon=(/-180,180/)
154    zoom_lat=(/-90,90/)
155    !
156    !Config Key   = LIMIT_WEST
157    !Config Desc  = Western limit of region
158    !Config If    = [-]
159    !Config Def   = -180.
160    !Config Help  = Western limit of the region we are
161    !Config         interested in. Between -180 and +180 degrees
162    !Config         The model will use the smalest regions from
163    !Config         region specified here and the one of the forcing file.
164    !Config Units = [Degrees]
165    !-
166    CALL getin_p('LIMIT_WEST',zoom_lon(1))
167    !-
168    !Config Key   = LIMIT_EAST
169    !Config Desc  = Eastern limit of region
170    !Config If    = [-]
171    !Config Def   = 180.
172    !Config Help  = Eastern limit of the region we are
173    !Config         interested in. Between -180 and +180 degrees
174    !Config         The model will use the smalest regions from
175    !Config         region specified here and the one of the forcing file.
176    !Config Units = [Degrees]
177    !-
178    CALL getin_p('LIMIT_EAST',zoom_lon(2))
179    !-
180    !Config Key   = LIMIT_NORTH
181    !Config Desc  = Northern limit of region
182    !Config If    = [-]
183    !Config Def   = 90.
184    !Config Help  = Northern limit of the region we are
185    !Config         interested in. Between +90 and -90 degrees
186    !Config         The model will use the smalest regions from
187    !Config         region specified here and the one of the forcing file.
188    !Config Units = [Degrees]
189    !-
190    CALL getin_p('LIMIT_NORTH',zoom_lat(2))
191    !-
192    !Config Key   = LIMIT_SOUTH
193    !Config Desc  = Southern limit of region
194    !Config If    = [-]
195    !Config Def   = -90.
196    !Config Help  = Southern limit of the region we are
197    !Config         interested in. Between 90 and -90 degrees
198    !Config         The model will use the smalest regions from
199    !Config         region specified here and the one of the forcing file.
200    !Config Units = [Degrees]
201    !-
202    CALL getin_p('LIMIT_SOUTH',zoom_lat(1))
203    IF ( (zoom_lon(1)+180 < EPSILON(zoom_lon(1))) .AND. (zoom_lon(2)-180 < EPSILON(zoom_lon(2))) .AND.&
204        &(zoom_lat(1)+90 < EPSILON(zoom_lat(1))) .AND. (zoom_lat(2)-90 < EPSILON(zoom_lat(2))) ) THEN
205       !
206       !Config Key   = WEST_EAST
207       !Config Desc  = Longitude interval to use from the forcing data
208       !Config If    = [-]
209       !Config Def   = -180, 180
210       !Config Help  = This function allows to zoom into the forcing data
211       !Config Units = [degrees east]
212       !-
213       CALL getin('WEST_EAST', zoom_lon)
214       !
215       !Config Key   = SOUTH_NORTH
216       !Config Desc  = Latitude interval to use from the forcing data
217       !Config If    = [-]
218       !Config Def   = -90, 90
219       !Config Help  = This function allows to zoom into the forcing data
220       !Config Units = [degrees north]
221       !-
222       CALL getin('SOUTH_NORTH', zoom_lat)
223    ENDIF
224    !
225    ! Get the finer coarser grid
226    !
227    !
228    ! We see if the user has specified the resolution, in number of
229    ! points he wishes to use.
230    !
231    IF ( buildcontfrac ) THEN
232       IF ( is_root_prc) THEN
233          iim=-1
234          jjm=-1
235          !Config Key   = IIM_LON
236          !Config Desc  = Number of points in longitude
237          !Config If    = [-]
238          !Config Def   = -1
239          !Config Help  = This allows to select the resolution at which the land/sea mask should be built.
240          !Config Units = -
241          CALL getin('IIM_LON', iim)
242          !Config Key   = IIM_LAT
243          !Config Desc  = Number of points in latitude
244          !Config If    = [-]
245          !Config Def   = -1
246          !Config Help  = This allows to select the resolution at which the land/sea mask should be built.
247          !Config Units = -
248          CALL getin('JJM_LAT', jjm)
249          !
250          IF (iim > 0 .AND. jjm > 0 ) THEN
251             !
252             ALLOCATE(lon(iim,jjm), lat(iim,jjm), orog(iim,jjm))
253             ! Trung: Calculate resolution
254             dx = (MAXVAL(zoom_lon)-MINVAL(zoom_lon))/REAL(iim, r_std)
255             dy = (MAXVAL(zoom_lat)-MINVAL(zoom_lat))/REAL(jjm, r_std)
256             WRITE (*,*) "Resolution in longitude [degrees] dx = ",dx
257             WRITE (*,*) "Resolution in laltiude [degrees] dy = ",dy
258             !
259             ! Generate lon and lat
260             !
261             DO i=1,iim
262                lon(i,:) = MINVAL(zoom_lon) + (i-1)*dx + dx/2.0 
263             ENDDO
264             DO j=1,jjm
265                lat(:,j) = MINVAL(zoom_lat) + (j-1)*dy + dy/2.0
266             ENDDO
267             !
268             ! Interpolate to new orography
269             !
270             CALL interpol(iim_full, jjm_full, lon_full, lat_full, etopo, iim, jjm, lon, lat, orog)
271             !
272             ! Write orog to netcdf for checking (deleted)
273             !
274          ELSE
275             ! Just transfer the data over the zoomed area
276             ibeg=MINLOC(ABS(lon_full(:,1)-MINVAL(zoom_lon)))
277             iend=MINLOC(ABS(lon_full(:,1)-MAXVAL(zoom_lon)))
278             jbeg=MINLOC(ABS(lat_full(1,:)-MINVAL(zoom_lat)))
279             jend=MINLOC(ABS(lat_full(1,:)-MAXVAL(zoom_lat)))
280             iim = (iend(1)-ibeg(1))+1
281             jjm = (jend(1)-jbeg(1))+1
282             !
283             !
284             ALLOCATE(lon(iim,jjm), lat(iim,jjm), orog(iim,jjm), contfrac(iim,jjm))
285             !
286             DO j=1,jjm
287                lon(:,j) = lon_full(ibeg(1):iend(1),j)
288             ENDDO
289             DO i=1,iim
290                lat(i,:) = lat_full(i,jbeg(1):jend(1))
291             ENDDO
292             ! Some treatment for the orogography
293             orog(:,:) = 0.0
294             DO i=ibeg(1),iend(1)
295                DO j=jbeg(1),jend(1)
296                   orog(i-ibeg(1)+1,j-jbeg(1)+1) = MAX(etopo(i,j), 0.0)
297                ENDDO
298             ENDDO
299             !
300             ! Generate the contfrac field
301             !
302             contfrac(:,:)= 0.0
303             DO i=1,iim
304                DO j=1,jjm
305                   IF ( orog(i,j) > 0 ) THEN
306                      contfrac(i,j) = 1.0
307                   ENDIF
308                ENDDO
309             ENDDO
310          ENDIF
311          !
312       ENDIF
313       !
314    ELSE
315       !
316       ! Just zoom into the contfrac if needed and build a dummy orog
317       !
318       IF ( is_root_prc) THEN
319          IF ( MINVAL(zoom_lon) > -180 .OR. MAXVAL(zoom_lon) < 180 .OR.&
320               & MINVAL(zoom_lat) > -90 .OR. MAXVAL(zoom_lat) < 90 ) THEN
321             ! Just transfer the data over the zoomed area
322             ibeg=MINLOC(ABS(lon_full(:,1)-MINVAL(zoom_lon)))
323             iend=MINLOC(ABS(lon_full(:,1)-MAXVAL(zoom_lon)))
324             jbeg=MINLOC(ABS(lat_full(1,:)-MINVAL(zoom_lat)))
325             jend=MINLOC(ABS(lat_full(1,:)-MAXVAL(zoom_lat)))
326             iim = (iend(1)-ibeg(1))+1
327             jjm = (jend(1)-jbeg(1))+1
328             !
329             !
330             ALLOCATE(lon(iim,jjm), lat(iim,jjm), orog(iim,jjm), contfrac(iim,jjm))
331             !
332             DO j=1,jjm
333                lon(:,j) = lon_full(ibeg(1):iend(1),j)
334             ENDDO
335             DO i=1,iim
336                lat(i,:) = lat_full(i,jbeg(1):jend(1))
337             ENDDO
338             ! Some treatment for the orogography
339             orog(:,:) = 0.0
340             contfrac(:,:) = 0.0
341             DO i=ibeg(1),iend(1)
342                DO j=jbeg(1),jend(1)
343                   orog(i-ibeg(1)+1,j-jbeg(1)+1) = 1.0
344                   contfrac(i-ibeg(1)+1,j-jbeg(1)+1) = contfrac_full(i,j)
345                ENDDO
346             ENDDO
347             !
348          ELSE
349             iim = iim_full
350             jjm = jjm_full
351             ALLOCATE(lon(iim,jjm), lat(iim,jjm), orog(iim,jjm), contfrac(iim,jjm))
352             lon(:,:) = lon_full(:,:)
353             lat(:,:) = lat_full(:,:)
354             contfrac(:,:) = contfrac_full(:,:)
355             orog(:,:) = contfrac_full(:,:)+100
356          ENDIF
357       ENDIF
358       !
359    ENDIF
360    !
361    IF ( is_root_prc) THEN
362       WRITE(*,*) "Dimensions : ", iim, jjm
363       WRITE(*,*) MINVAL(lon), ' < LON < ', MAXVAL(lon)
364       WRITE(*,*) MINVAL(lat), ' < LAT < ', MAXVAL(lat) 
365       WRITE(*,*) MINVAL(orog), " < OROG  < ", MAXVAL(orog)
366       WRITE(*,*) MINVAL(contfrac), " < contfrac  < ", MAXVAL(contfrac)
367    ENDIF
368    !
369    ! Free some memory
370    !
371    IF (ALLOCATED(lon_full)) DEALLOCATE(lon_full)
372    IF (ALLOCATED(lat_full)) DEALLOCATE(lat_full)
373    IF (ALLOCATED(etopo)) DEALLOCATE(etopo)
374    IF (ALLOCATED(contfrac_full)) DEALLOCATE(contfrac_full)
375    !
376    ! Distribute the information to all processors
377    !
378    CALL bcast(iim)
379    CALL bcast(jjm)
380    IF ( .NOT. ALLOCATED(lon)) ALLOCATE(lon(iim,jjm))
381    IF ( .NOT. ALLOCATED(lat)) ALLOCATE(lat(iim,jjm))
382    IF ( .NOT. ALLOCATED(orog)) ALLOCATE(orog(iim,jjm))
383    IF ( .NOT. ALLOCATED(contfrac)) ALLOCATE(contfrac(iim,jjm))
384    CALL bcast(lon)
385    CALL bcast(lat)
386    CALL bcast(orog)
387    CALL bcast(contfrac)
388    !
389    ! Count the number of land points
390    !
391    nbland = COUNT(contfrac > 0.0)
392    !
393    iim_out = iim
394    jjm_out = jjm
395    nbland_out = nbland
396    !
397  END SUBROUTINE getlandseamask_init
398  !!
399  !================================================================================================================================
400  !! SUBROUTINE   : getlandseamask_read
401  !!
402  !> BRIEF         Returns the land/sea mask
403  !!
404  !! DESCRIPTION  : None
405  !!
406  !! RECENT CHANGE(S): None
407  !!
408  !! MAIN OUTPUT VARIABLE(S):
409  !!
410  !! REFERENCES   : None
411  !!
412  !! FLOWCHART    : None
413  !! \n
414  !_================================================================================================================================
415  !!
416  SUBROUTINE getlandseamask_read(lon_out, lat_out, mask_out, orog_out)
417    !
418    ! ARGUMENTS
419    !
420    REAL(r_std), DIMENSION(:,:), INTENT(out) :: lon_out, lat_out
421    REAL(r_std), DIMENSION(:,:), INTENT(out) :: mask_out, orog_out
422    !
423    ! LOCAL
424    !
425    lon_out(:,:) = lon(:,:)
426    lat_out(:,:) = lat(:,:)
427    mask_out(:,:) = contfrac(:,:)
428    orog_out(:,:) = orog(:,:)
429    !
430  END SUBROUTINE getlandseamask_read
431  !!
432  !!
433  !
434  !================================================================================================================================
435  !! SUBROUTINE   : opencontfrac
436  !!
437  !> BRIEF         This subroutine gets the dimensions of the fields in the filename. Lat and lon in particular.
438  !!
439  !! DESCRIPTION  : This routines caters for ORCHIDEE as well as the WRF geo files.
440  !!
441  !! RECENT CHANGE(S): None
442  !!
443  !! MAIN OUTPUT VARIABLE(S):
444  !!
445  !! REFERENCES   : None
446  !!
447  !! FLOWCHART    : None
448  !! \n
449  !_================================================================================================================================
450  SUBROUTINE opencontfrac(filename, iim_full, jjm_full, fid)
451    !
452    ! ARGUMENTS
453    !
454    CHARACTER(LEN=*), INTENT(in) :: filename
455    INTEGER(i_std), INTENT(out) :: iim_full, jjm_full, fid
456    !
457    INTEGER(i_std) :: iret, i, lll, ndims, nvars
458    CHARACTER(LEN=20) :: dimname
459    !
460    iret = NF90_OPEN(filename, NF90_NOWRITE, fid)
461    IF (iret /= NF90_NOERR) THEN
462       WRITE(*,*) "==>",trim(nf90_strerror(iret))
463       WRITE(*,*) "Error opening ", filename
464       STOP
465    ENDIF
466    iret = NF90_INQUIRE (fid, nDimensions=ndims, nVariables=nvars)
467    DO i=1,ndims
468       iret = NF90_INQUIRE_DIMENSION (fid, i, name=dimname, len=lll)
469       CALL change_to_lowercase(dimname)
470       IF ( INDEX(dimname,"lon") > 0 ) THEN
471          iim_full = lll
472       ELSE IF (INDEX(dimname,"west_east") > 0 .AND. LEN_TRIM(dimname) == LEN_TRIM("west_east")) THEN
473          iim_full = lll
474       ELSE IF ( INDEX(dimname,"lat") > 0 ) THEN
475          jjm_full = lll
476       ELSE IF (INDEX(dimname,"south_north") > 0 .AND. LEN_TRIM(dimname) == LEN_TRIM("south_north")) THEN
477          jjm_full = lll
478       ENDIF
479    ENDDO
480    !
481    WRITE(*,*) "XXX iimf, jjmf = ", iim_full, jjm_full
482    !
483  END SUBROUTINE opencontfrac
484
485  !================================================================================================================================
486  !! SUBROUTINE   : readcontfrac
487  !!
488  !> BRIEF         This subroutine gets the topography out of the file
489  !!
490  !! DESCRIPTION  : This routines caters for ORCHIDEE as well as the WRF geo files.
491  !!
492  !! RECENT CHANGE(S): None
493  !!
494  !! MAIN OUTPUT VARIABLE(S):
495  !!
496  !! REFERENCES   : None
497  !!
498  !! FLOWCHART    : None
499  !! \n
500  !_================================================================================================================================
501  SUBROUTINE readcontfrac(fid, filename, iimf, jjmf, lon, lat, contf)
502    !
503    ! ARGUMENTS
504    !
505    INTEGER(i_std), INTENT(in)               :: fid
506    CHARACTER(LEN=*), INTENT(in)             :: filename
507    INTEGER(i_std), INTENT(in)               :: iimf, jjmf
508    REAL(r_std), DIMENSION(:,:), INTENT(out) :: lon, lat
509    REAL(r_std), DIMENSION(:,:), INTENT(out) :: contf
510    !
511    ! LOCAL
512    !
513    INTEGER(i_std) :: iret, ndims, nvars, xid, yid, vid, wid
514    INTEGER(i_std) :: xndims, yndims, vndims, wndims
515    INTEGER(i_std) :: iv, ndimsvar, i, j
516    INTEGER(i_std), DIMENSION(4) :: dimids
517    CHARACTER(LEN=60) :: varname, units
518    CHARACTER(LEN=3)  :: model='ORC'
519    INTEGER(i_std), DIMENSION(3) :: start, count
520    REAL(r_std), ALLOCATABLE, DIMENSION(:) :: vtmp
521    !
522    iret = NF90_INQUIRE (fid, nDimensions=ndims, nVariables=nvars)
523    !
524    xid = -1
525    yid = -1
526    vid = -1
527    wid = -1
528    !
529    DO iv=1,nvars
530       !
531       iret = NF90_INQUIRE_VARIABLE(fid, iv, name=varname, ndims=ndimsvar, dimids=dimids)
532       iret = NF90_GET_ATT(fid, iv, 'units', units)
533       CALL change_to_lowercase(units)
534       CALL change_to_lowercase(varname)
535       !
536       IF ( INDEX(units,"east") > 0 .OR. INDEX(varname, "xlong_m") > 0 ) THEN
537          xid = iv
538          xndims = ndimsvar
539          IF ( INDEX(varname, "xlong_m") > 0 ) THEN
540             model='WRF'
541          ENDIF
542       ELSE IF (INDEX(units,"north") > 0 .OR. INDEX(varname, "xlat_m") > 0 ) THEN
543          yid = iv
544          yndims = ndimsvar
545          IF ( INDEX(varname, "xlat_m") > 0 ) THEN
546             model='WRF'
547          ENDIF
548       ENDIF
549       !
550       IF ( INDEX(varname,"contfrac") > 0 ) THEN
551          vid = iv
552          vndims = ndimsvar
553       ENDIF
554       IF ( INDEX(varname,"landmask") > 0 ) THEN
555          wid = iv
556          wndims = ndimsvar
557       ENDIF
558       !
559    ENDDO
560    !
561    ! Get longitude
562    !
563    IF ( xid > 0 ) THEN
564       IF ( xndims == 1 ) THEN
565          ALLOCATE(vtmp(iimf))
566          iret = NF90_GET_VAR(fid, xid, vtmp)
567          DO j=1,jjmf
568             lon(:,j) = vtmp(:)
569          ENDDO
570          DEALLOCATE(vtmp)
571       ELSE IF ( xndims == 2 ) THEN
572          iret = NF90_GET_VAR(fid, xid, lon)
573       ELSE IF ( xndims == 3 ) THEN
574          start = (/1,1,1/)
575          count = (/iimf,jjmf,1/)
576          iret = NF90_GET_VAR(fid, xid, lon, start, count)
577       ELSE
578          WRITE(*,*) "Unforeseen number of dimensions for longitude ", xndims
579          STOP
580       ENDIF
581    ELSE
582       WRITE(*,*) "could not find the longitude in file ", filename
583       STOP
584    ENDIF
585    !
586    ! Get latitude
587    !
588    IF ( yid > 0 ) THEN
589       IF ( yndims == 1 ) THEN
590          ALLOCATE(vtmp(jjmf))
591          iret = NF90_GET_VAR(fid, xid, vtmp)
592          DO i=1,iimf
593             lat(i,:) = vtmp(:)
594          ENDDO
595          DEALLOCATE(vtmp)
596       ELSE IF ( yndims == 2 ) THEN
597          iret = NF90_GET_VAR(fid, yid, lat)
598       ELSE IF ( yndims == 3 ) THEN
599          start = (/1,1,1/)
600          count = (/iimf,jjmf,1/)
601          iret = NF90_GET_VAR(fid, yid, lat, start, count)
602       ELSE
603          WRITE(*,*) "Unforeseen number of dimensions for latitude ", yndims
604          STOP
605       ENDIF
606    ELSE
607       WRITE(*,*) "could not find the latitude in file ", filename
608       STOP
609    ENDIF
610    !
611    ! Get the variable
612    !
613    contf(:,:) = 0.0
614    IF ( model == "ORC" ) THEN
615       IF ( vid > 0 ) THEN
616          iret = NF90_GET_VAR(fid, vid, contf)
617          IF (iret /= NF90_NOERR) THEN
618             WRITE(*,*) "==>",trim(nf90_strerror(iret))
619             WRITE(*,*) "Cannot read variable contfrac"
620             STOP
621          ENDIF
622       ELSE
623          WRITE(*,*) "could not find the variable contfrac in file ", filename
624          STOP
625       ENDIF
626    ELSE IF ( model == "WRF" ) THEN
627       !
628       ! This is a time varying variable so start and count needs to be specified
629       !
630       start = (/1,1,1/)
631       count = (/iimf,jjmf,1/)
632       IF ( wid > 0 ) THEN
633          iret = NF90_GET_VAR(fid, wid, contf, start, count)
634          IF (iret /= NF90_NOERR) THEN
635             WRITE(*,*) "==>",trim(nf90_strerror(iret))
636             WRITE(*,*) "Cannot read variable landmask"
637             STOP
638          ENDIF
639       ELSE
640          WRITE(*,*) "could not find the variable contfrac in file ", filename
641          STOP
642       ENDIF
643       WRITE(*,*) "XXX", MINVAL(contf), " < contf < ", MAXVAL(contf)
644    ELSE
645       WRITE(*,*) "Unknown model = ", model
646       STOP
647    ENDIF
648    !
649    ! Delete undef values
650    !
651    DO i=1,iimf
652       DO j=1,jjmf
653          IF ( contf(i,j) > 1.0 .OR. contf(i,j) < 0.0 ) THEN
654             contf(i,j) = 0.0
655          ENDIF
656       ENDDO
657    ENDDO
658    !
659    iret = NF90_CLOSE(fid)
660    !
661  END SUBROUTINE readcontfrac
662  !!
663  !
664  !================================================================================================================================
665  !! SUBROUTINE   : topo_getsize
666  !!
667  !> BRIEF         This subroutine gets the size of the topography field
668  !!
669  !! DESCRIPTION  : None
670  !!
671  !! RECENT CHANGE(S): None
672  !!
673  !! MAIN OUTPUT VARIABLE(S):
674  !!
675  !! REFERENCES   : None
676  !!
677  !! FLOWCHART    : None
678  !! \n
679  !_================================================================================================================================
680  SUBROUTINE topo_getsize(filename, iim_full, jjm_full)
681    !
682    USE defprec
683    USE netcdf
684    !
685    ! ARGUMENTS
686    !
687    CHARACTER(LEN=*), INTENT(in) :: filename
688    INTEGER(i_std), INTENT(out)  :: iim_full, jjm_full
689    !
690    ! LOCAL
691    !
692    INTEGER(i_std) :: fid, iret, xid, yid, ndims, nvars, i, lll
693    CHARACTER(LEN=20) :: dimname
694    !
695    iret = NF90_OPEN(filename, NF90_NOWRITE, fid)
696    IF (iret /= NF90_NOERR) THEN
697       WRITE(*,*) "Error opening ", filename
698       STOP
699    ENDIF
700    iret = NF90_INQUIRE (fid, nDimensions=ndims, nVariables=nvars)
701    DO i=1,ndims
702       iret = NF90_INQUIRE_DIMENSION (fid, i, name=dimname, len=lll)
703       CALL change_to_lowercase(dimname)
704       IF ( INDEX(dimname,"x") > 0 ) THEN
705          iim_full = lll
706       ELSE IF ( INDEX(dimname,"y") > 0 ) THEN
707          jjm_full = lll
708       ENDIF
709    ENDDO
710    !
711    iret = NF90_CLOSE(fid)
712    !
713  END SUBROUTINE topo_getsize
714  !================================================================================================================================
715  !! SUBROUTINE   : topo_getvar
716  !!
717  !> BRIEF         This subroutine gets the topography out of the file
718  !!
719  !! DESCRIPTION  : None
720  !!
721  !! RECENT CHANGE(S): None
722  !!
723  !! MAIN OUTPUT VARIABLE(S):
724  !!
725  !! REFERENCES   : None
726  !!
727  !! FLOWCHART    : None
728  !! \n
729  !_================================================================================================================================
730  SUBROUTINE topo_getvar(filename, iim, jjm, lon, lat, etopo)
731    !
732    USE defprec
733    USE netcdf
734    !
735    ! ARGUMENTS
736    !
737    CHARACTER(LEN=*), INTENT(in)                 :: filename
738    INTEGER(i_std), INTENT(in)                   :: iim, jjm
739    REAL(r_std), DIMENSION(iim,jjm), INTENT(out) :: lon 
740    REAL(r_std), DIMENSION(iim,jjm), INTENT(out) :: lat
741    REAL(r_std), DIMENSION(iim,jjm), INTENT(out) :: etopo
742    !
743    ! LOCAL
744    !
745    INTEGER(i_std) :: fid, iret, xid, yid, vid
746    INTEGER(i_std) :: i, j, iv, ndimsvar, ndims, nvars
747    CHARACTER(LEN=20) :: varname, units
748    INTEGER(i_std), DIMENSION(1) :: il
749    REAL(r_std), ALLOCATABLE, DIMENSION(:)   :: lon_loc, lon_read, lat_read
750    REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: etopo_loc
751    INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: dimids, londim_ids, latdim_ids
752    !
753    ALLOCATE(lon_loc(iim))
754    ALLOCATE(lon_read(iim))
755    ALLOCATE(lat_read(jjm))
756    ALLOCATE(etopo_loc(iim,jjm))
757    ALLOCATE(dimids(NF90_MAX_VAR_DIMS), londim_ids(NF90_MAX_VAR_DIMS), latdim_ids(NF90_MAX_VAR_DIMS))
758    !
759    iret = NF90_OPEN(filename, NF90_NOWRITE, fid)
760    IF (iret /= NF90_NOERR) THEN
761       WRITE(*,*) "Error opening ", filename
762       STOP
763    ENDIF
764    !
765    iret = NF90_INQUIRE (fid, nDimensions=ndims, nVariables=nvars)
766    !
767    ! Identify the variables
768    !
769    xid = -1
770    yid = -1
771    vid = -1
772    !
773    DO iv=1,nvars
774       iret = NF90_INQUIRE_VARIABLE(fid, iv, name=varname, ndims=ndimsvar, dimids=dimids)
775       iret = NF90_GET_ATT(fid, iv, 'units', units)
776       CALL change_to_lowercase(units)
777       !
778       IF ( INDEX(units,"east") > 0 ) THEN
779          xid = iv
780       ELSE IF (INDEX(units,"north") > 0 ) THEN
781          yid = iv
782       ELSE IF ( INDEX(units,"meters") > 0 ) THEN
783          vid = iv
784       ENDIF
785    ENDDO
786    !
787    ! Get variables now that they have been identified
788    !
789    IF ( xid > 0 ) THEN
790       iret = NF90_GET_VAR(fid, xid, lon_read)
791    ELSE
792       WRITE(*,*) "cound not find the longitude in file ", filename
793       STOP
794    ENDIF
795    IF ( yid > 0 ) THEN
796       iret = NF90_GET_VAR(fid, yid, lat_read)
797    ELSE
798       WRITE(*,*) "cound not find the latitude in file ", filename
799       STOP
800    ENDIF
801    IF ( vid > 0 ) THEN
802       iret = NF90_GET_VAR(fid, vid, etopo_loc)
803    ELSE
804       WRITE(*,*) "cound not find the topography in file ", filename
805       STOP
806    ENDIF
807    iret = NF90_CLOSE(fid)
808    !
809    ! Shift the topography around the longitude to get a -180:180 range
810    ! Only if needed !
811    !
812    IF ( MAXVAL(lon_read) > 180. ) THEN
813       !
814       lon_loc(:) = lon_read(:)
815       !
816       DO i=1,iim
817          IF (lon_loc(i) > 180. ) THEN
818             lon_loc(i) = lon_loc(i)-360.
819          ENDIF
820       ENDDO
821       !
822       DO i=1,iim
823          il = MINLOC(lon_loc)
824          lon_read(i) = lon_loc(il(1))
825          etopo(i,:) = etopo_loc(il(1),:)
826          lon_loc(il(1)) = 9999999.999999
827       ENDDO
828       !
829       DO j=1,jjm
830          lon(:,j) = lon_read(:)
831       ENDDO
832       DO i=1,iim
833          lat(i,:) = lat_read(:)
834       ENDDO
835    ELSE
836       DO j=1,jjm
837          lon(:,j) = lon_read(:)
838       ENDDO
839       DO i=1,iim
840          lat(i,:) = lat_read(:)
841       ENDDO
842       etopo(:,:) = etopo_loc(:,:)
843    ENDIF
844    !
845    !
846    DEALLOCATE(lon_read)
847    DEALLOCATE(lon_loc)
848    DEALLOCATE(lat_read)
849    DEALLOCATE(etopo_loc)
850    !
851  END SUBROUTINE topo_getvar
852  !
853  !================================================================================================================================
854  !! SUBROUTINE   : interpol
855  !!
856  !> BRIEF         This subroutine interpolates topography to new resolution.
857  !!
858  !! DESCRIPTION  : None
859  !!
860  !! RECENT CHANGE(S): None
861  !!
862  !! MAIN OUTPUT VARIABLE(S):
863  !!
864  !! REFERENCES   : None
865  !!
866  !! FLOWCHART    : None
867  !! \n
868  !_================================================================================================================================
869
870  SUBROUTINE interpol(iim_fine, jjm_fine, lon_fine, lat_fine, orog, iim, jjm, lon, lat, neworog)
871    !
872    USE constantes_var
873    USE ioipsl_para
874    !
875    IMPLICIT NONE
876    ! Parameter
877    !INTEGER(i_std), PARAMETER                                  :: nbvmax = 800
878    INTEGER(i_std), PARAMETER                                  :: nbvmax = 30000
879    !
880    INTEGER(i_std)                                             :: iim_fine, jjm_fine, iim, jjm
881    REAL(r_std), DIMENSION (iim_fine,jjm_fine)                 :: lon_fine
882    REAL(r_std), DIMENSION (iim_fine,jjm_fine)                 :: lat_fine
883    REAL(r_std), DIMENSION (iim_fine,jjm_fine)                 :: orog
884    REAL(r_std), DIMENSION (iim,jjm)                           :: lon
885    REAL(r_std), DIMENSION (jjm,iim)                           :: lat
886    REAL(r_std), DIMENSION (iim,jjm)                           :: neworog
887    !
888    REAL(r_std), DIMENSION (nbvmax)                            :: area, topo
889    !
890    INTEGER(i_std)                                             :: ip, jp, ig, jg, i, fopt, lastjp
891    REAL(r_std)                                                :: dx, dy, ax, ay, ox, lon_up, lon_low, lat_up, lat_low, sgn
892    REAL(r_std)                                                :: totarea, landarea, height
893    REAL(r_std), ALLOCATABLE, DIMENSION (:,:)                  :: lon_ful, lat_ful, laup_rel, loup_rel, lalow_rel, lolow_rel
894    REAL(r_std)                                                :: lonrel, lolowrel, louprel, coslat
895    !
896    ! Allocate arrays
897    !
898    ALLOCATE (laup_rel(iim_fine,jjm_fine))
899    ALLOCATE (loup_rel(iim_fine,jjm_fine))
900    ALLOCATE (lalow_rel(iim_fine,jjm_fine))
901    ALLOCATE (lolow_rel(iim_fine,jjm_fine))
902    ALLOCATE (lat_ful(iim_fine+2,jjm_fine+2))
903    ALLOCATE (lon_ful(iim_fine+2,jjm_fine+2))
904    !
905    ! Duplicate the border assuming we have a global grid going from west to east
906    !
907    DO jp=2,jjm_fine+1
908       lon_ful(2:iim_fine+1,jp) = lon_fine(1:iim_fine,jp)
909    ENDDO
910    DO ip=2,iim_fine+1
911       lat_ful(ip,2:jjm_fine+1) = lat_fine(ip,1:jjm_fine)
912    ENDDO
913    !
914    IF ( lon_fine(iim_fine,1) .LT. lon_ful(2,2)) THEN
915       lon_ful(1,2:jjm_fine+1) = lon_fine(iim_fine,:)
916       lat_ful(1,2:jjm_fine+1) = lat_fine(1,1:jjm_fine)
917    ELSE
918       lon_ful(1,2:jjm_fine+1) = lon_fine(iim_fine,:)-360.
919       lat_ful(1,2:jjm_fine+1) = lat_fine(1,1:jjm_fine)
920    ENDIF
921
922    IF ( lon_fine(1,1) .GT. lon_ful(iim_fine+1,2)) THEN
923       lon_ful(iim_fine+2,2:jjm_fine+1) = lon_fine(1,1:jjm_fine)
924       lat_ful(iim_fine+2,2:jjm_fine+1) = lat_fine(iim_fine,1:jjm_fine)
925    ELSE
926       lon_ful(iim_fine+2,2:jjm_fine+1) = lon_fine(1,1:jjm_fine)+360
927       lat_ful(iim_fine+2,2:jjm_fine+1) = lat_fine(iim_fine,1:jjm_fine)
928    ENDIF
929    !
930    sgn = lat_fine(1,1)/ABS(lat_fine(1,1))
931    lat_ful(2:iim_fine+1,1) = sgn*180 - lat_fine(1,1)
932    sgn = lat_fine(iim_fine,jjm_fine)/ABS(lat_fine(iim_fine,jjm_fine))
933    lat_ful(2:iim_fine+1,jjm_fine+2) = sgn*180 - lat_fine(iim_fine,jjm_fine)
934    lat_ful(1,1) = lat_ful(iim_fine+1,1)
935    lat_ful(iim_fine+2,1) = lat_ful(2,1)
936    lat_ful(1,jjm_fine+2) = lat_ful(iim_fine+1,jjm_fine+2)
937    lat_ful(iim_fine+2,jjm_fine+2) = lat_ful(2,jjm_fine+2)
938    !
939    ! Add the longitude lines to the top and bottom
940    !
941    lon_ful(:,1) = lon_ful(:,2) 
942    lon_ful(:,jjm_fine+2) = lon_ful(:,jjm_fine+1) 
943    !
944    ! Get the upper and lower limits of each grid box
945    !
946    DO ip=1,iim_fine
947       DO jp=1,jjm_fine
948          loup_rel(ip,jp) =MAX(0.5*(lon_ful(ip,jp+1)+lon_ful(ip+1,jp+1)), 0.5*(lon_ful(ip+1,jp+1)+lon_ful(ip+2,jp+1)))
949          lolow_rel(ip,jp) =MIN(0.5*(lon_ful(ip,jp+1)+lon_ful(ip+1,jp+1)), 0.5*(lon_ful(ip+1,jp+1)+lon_ful(ip+2,jp+1)))
950          laup_rel(ip,jp) =MAX(0.5*(lat_ful(ip+1,jp)+lat_ful(ip+1,jp+1)), 0.5*(lat_ful(ip+1,jp+1)+lat_ful(ip+1,jp+2)))
951          lalow_rel(ip,jp) =MIN(0.5*(lat_ful(ip+1,jp)+lat_ful(ip+1,jp+1)), 0.5*(lat_ful(ip+1,jp+1)+lat_ful(ip+1,jp+2)))
952       ENDDO
953    ENDDO
954    !
955    ! Now we take each grid point and find out which values from the forcing we need to average
956    !
957    dx = ABS(lon(2,1) - lon(1,1))
958    dy = ABS(lat(1,2) - lat(1,1))
959    !
960    DO ig = 1,iim
961       DO jg = 1,jjm
962          !
963          !
964          lon_up = lon(ig,jg) + dx/2.0
965          lon_low = lon(ig,jg) - dx/2.0
966          !
967          !
968          lat_up = lat(ig,jg) + dy/2.0
969          lat_low = lat(ig,jg) - dy/2.0
970          !
971          !  Find the grid boxes from the data that go into the model's boxes
972          !  We still work as if we had a regular grid ! Well it needs to be localy regular so
973          !  so that the longitude at the latitude of the last found point is close to the one of the next point.
974          !
975          totarea = 0.0
976          landarea = 0.0
977          height = 0.0
978          !
979          fopt = 0
980          lastjp = 1
981          DO ip=1,iim_fine
982             !
983             DO jp = 1, jjm_fine
984                !  Either the center of the data grid point is in the interval of the model grid or
985                !  the East and West limits of the data grid point are on either sides of the border of
986                !  the data grid.
987                !
988                !  To do that correctly we have to check if the grid box sits on the date-line.
989                !
990                IF ( lon_low < -180.0 ) THEN
991                   lonrel = MOD( lon_fine(ip,jp) - 360.0, 360.0)
992                   lolowrel = MOD( lolow_rel(ip,jp) - 360.0, 360.0)
993                   louprel = MOD( loup_rel(ip,jp) - 360.0, 360.0)
994                   !
995                ELSE IF ( lon_up > 180.0 ) THEN
996                   lonrel = MOD( 360. - lon_fine(ip,jp), 360.0)
997                   lolowrel = MOD( 360. - lolow_rel(ip,jp), 360.0)
998                   louprel = MOD( 360. - loup_rel(ip,jp), 360.0)
999                ELSE
1000                   lonrel = lon_fine(ip,jp)
1001                   lolowrel = lolow_rel(ip,jp)
1002                   louprel = loup_rel(ip,jp)
1003                ENDIF
1004                !
1005                !
1006                IF ( lonrel > lon_low .AND. lonrel < lon_up .OR. &
1007                     & lolowrel < lon_low .AND.  louprel > lon_low .OR. &
1008                     & lolowrel < lon_up  .AND.  louprel > lon_up ) THEN
1009                   ! 
1010                   !
1011                   ! Now that we have the longitude let us find the latitude
1012                   !
1013                   IF ( lat_fine(ip,jp) > lat_low .AND. lat_fine(ip,jp) < lat_up .OR. &
1014                        & lalow_rel(ip,jp) < lat_low .AND. laup_rel(ip,jp) > lat_low .OR.&
1015                        & lalow_rel(ip,jp) < lat_up .AND. laup_rel(ip,jp) > lat_up) THEN
1016                      !
1017                      lastjp = jp
1018                      !
1019                      fopt = fopt + 1
1020                      IF ( fopt .GT. nbvmax) THEN
1021                         WRITE(*,*) 'Please increase nbvmax in subroutine interpol', ig, jg
1022                         WRITE(*,*) "nbvmax = ", nbvmax
1023                         STOP
1024                      ELSE
1025                         !
1026                         ! If we sit on the date line we need to do the same transformations as above.
1027                         !
1028                         IF ( lon_low < -180.0 ) THEN
1029                            lolowrel = MOD( lolow_rel(ip,jp) - 360.0, 360.0)
1030                            louprel = MOD( loup_rel(ip,jp) - 360.0, 360.0)
1031                            !
1032                         ELSE IF ( lon_up > 180.0 ) THEN
1033                            lolowrel = MOD( 360. - lolow_rel(ip,jp), 360.0)
1034                            louprel = MOD( 360. - loup_rel(ip,jp), 360.0)
1035                         ELSE
1036                            lolowrel = lolow_rel(ip,jp)
1037                            louprel = loup_rel(ip,jp)
1038                         ENDIF
1039                         !
1040                         ! Get the area of the fine grid in the model grid
1041                         !
1042                         coslat = MAX( COS( lat_fine(ip,jp) * pi/180. ), 0.001 )
1043                         ax = (MIN(lon_up,louprel)-MAX(lon_low, lolowrel))*pi/180. * R_Earth * coslat
1044                         ay = (MIN(lat_up, laup_rel(ip,jp))-MAX(lat_low,lalow_rel(ip,jp)))*pi/180. * R_Earth
1045                         !
1046                         totarea = totarea + ax*ay
1047                         IF ( orog(ip, jp) > 0.0 ) THEN
1048                            landarea = landarea + ax*ay
1049                            height = height + ax*ay*orog(ip,jp)
1050                         ENDIF
1051                         !
1052                      ENDIF
1053                   ENDIF
1054                   !
1055                ENDIF
1056                !
1057             ENDDO
1058             !
1059          ENDDO
1060          !
1061          IF ( landarea > 0 .AND. landarea > 0.1*totarea ) THEN
1062             neworog(ig,jg) = height/landarea
1063          ELSE
1064             neworog(ig,jg) = 0.0
1065          ENDIF
1066          !
1067       ENDDO
1068       !
1069    ENDDO
1070  END SUBROUTINE interpol
1071
1072  SUBROUTINE change_to_lowercase (str)
1073    !---------------------------------------------------------------------
1074    !- Converts a string into lower case letters
1075    !---------------------------------------------------------------------
1076    IMPLICIT NONE
1077    !-
1078    CHARACTER(LEN=*) :: str
1079    !-
1080    INTEGER :: i,ic
1081    !---------------------------------------------------------------------
1082    DO i=1,LEN_TRIM(str)
1083       ic = IACHAR(str(i:i))
1084       IF ( (ic >= 65).AND.(ic <= 90) )  str(i:i) = ACHAR(ic+32)
1085    ENDDO
1086    !--------------------------
1087  END SUBROUTINE change_to_lowercase
1088
1089END MODULE getlandseamask
Note: See TracBrowser for help on using the repository browser.