source: trunk/LMDZ.MARS/libf/phymars/surfini_mod.F90 @ 4063

Last change on this file since 4063 was 3903, checked in by jbclement, 6 months ago

Mars PCM:
Case 'icelocationmode = 5' (default) is corrected compared to r2884. Now it reads 'watercaptag' from "startfi.nc" in any circumstances, unless it is missing, in which case 'icelocationmode' is set to 4.
JBC
File size: 24.6 KB
Line 
1MODULE surfini_mod
2
3implicit none
4
5! There are 4 different modes for ice distribution:
6! icelocationmode = 1 ---> based on data from surface.nc
7! icelocationmode = 2 ---> directly predefined for GCM resolutions 32x24 or 64x48
8! icelocationmode = 3 ---> based on logical relations for latitude and longitude
9! icelocationmode = 4 ---> predefined 64x48 but usable with every
10! resolution, and easily adaptable for dynamico
11! For visualisation : > /u/tnalmd/bin/watercaps gcm_txt_output_file
12integer :: icelocationmode = 5
13
14!$OMP THREADPRIVATE(icelocationmode)
15
16!=======================================================================
17contains
18!=======================================================================
19
20      SUBROUTINE surfini(ngrid,nslope,qsurf)
21
22      USE ioipsl_getin_p_mod, ONLY : getin_p
23      use netcdf
24      use tracer_mod, only: nqmx, noms
25      use geometry_mod, only: longitude, latitude, & ! in radians
26                          cell_area ! for watercaptag diagnosis
27      use surfdat_h, only: watercaptag, frost_albedo_threshold, &
28                          albedo_h2o_cap, inert_h2o_ice, albedodat, &
29                          albedice, dryness
30#ifndef MESOSCALE
31      use mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid
32      use mod_phys_lmdz_para, only : is_master, gather, scatter
33#endif
34      USE comcstfi_h, ONLY: pi
35      use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
36      use datafile_mod, only: datadir
37      use callkeys_mod, only: water
38
39      IMPLICIT NONE
40!=======================================================================
41!
42!   creation of caps for initial state
43!
44!=======================================================================
45!-----------------------------------------------------------------------
46!   Declarations:
47!   -------------
48
49      integer,intent(in) :: ngrid ! number of atmospheric columns
50      integer,intent(in) :: nslope ! number of sub-grid surfaces
51      real,intent(inout) :: qsurf(ngrid,nqmx,nslope) ! tracer on surface (kg/m2)
52
53      INTEGER ig,icap,iq,alternate,islope
54      REAL icedryness ! ice dryness
55
56      ! longwatercaptag is watercaptag. Trick for some compilers
57      LOGICAL, DIMENSION(100000) :: longwatercaptag
58
59      !in case icelocationmode == 1
60      INTEGER i,j
61      INTEGER     imd,jmd
62      PARAMETER   (imd=360,jmd=180)
63      REAL        zdata(imd,jmd)
64      REAL        zelat,zelon
65
66#ifndef MESOSCALE
67      INTEGER nb_ice(klon_glo,2)   ! number of counts | detected ice for GCM grid
68#endif
69      INTEGER latice(nbp_lat-1,2),lonice (nbp_lon,2) ! number of counts | detected ice along lat & lon axis
70
71      REAL step,count,ratiolat
72
73      INTEGER   ierr,nid,nvarid
74
75      REAL,SAVE :: min_icevalue = 500.
76
77!$OMP THREADPRIVATE(min_icevalue)
78
79      character(len=50) :: string = 'thermal'
80
81      character (len=100) :: zedatafile
82
83#ifdef MESOSCALE
84
85      do ig=1,ngrid
86
87         !write(*,*) "all qsurf to zero. dirty."
88         do iq=1,nqmx
89            do islope = 1,nslope
90              qsurf(ig,iq,islope)=0.  !! on jette les inputs GCM
91                          !! on regle juste watercaptag
92                          !! il faudrait garder les inputs GCM
93                          !! si elles sont consequentes
94            enddo
95         enddo
96         if ( ( latitude(ig)*180./pi .gt. 70. ) .and. ( albedodat(ig) .ge. 0.26   ) )  then
97                 write(*,*)"outlier ",ig
98                 watercaptag(ig)  = .true.
99                 dryness(ig)      = 1.
100                 albedodat(ig)    = albedo_h2o_cap  !! pour output
101         else
102                 watercaptag(ig)  = .false.
103                 dryness(ig)      = 1.
104         endif
105
106      enddo
107#endif
108! problem with nested precompiling flags
109
110#ifndef MESOSCALE
111      ! to handle parallel cases
112#if CPP_PARA
113      logical watercaptag_glo(klon_glo)
114      real dryness_glo(klon_glo)
115      real lati_glo(klon_glo)
116      real long_glo(klon_glo)
117#else
118      logical watercaptag_glo(ngrid)
119      real dryness_glo(ngrid)
120      real lati_glo(ngrid)
121      real long_glo(ngrid)
122#endif
123#endif
124
125#ifndef MESOSCALE
126
127!
128!=======================================================================
129! Initialize watercaptag (default is false)
130!      watercaptag_glo(:)=.false. !Already done in phyetat0 if needed
131
132!     water ice outliers
133!     ------------------------------------------
134
135      IF (water) THEN
136
137! Perennial H20 north cap defined by watercaptag=true (allows surface to be
138! hollowed by sublimation in vdifc).
139
140! We might not want albedodat to be modified because it is used to write
141! restart files. Instead, albedo is directly modified when needed (i.e.
142! if we have watercaptag and no co2 ice), below and in albedocaps.F90
143
144!       "Dryness coefficient" controlling the evaporation and
145!        sublimation from the ground water ice (close to 1)
146!        HERE, the goal is to correct for the fact
147!        that the simulated permanent water ice polar caps
148!        is larger than the actual cap and the atmospheric
149!        opacity not always realistic.
150
151         alternate = 0
152
153         if (ngrid .ne. 1) then
154           IF (icelocationmode .ne. 5) then
155             watercaptag(:) = .false.
156           ENDIF
157           longwatercaptag(:) = .false.
158         endif
159
160         write(*,*) "surfini: Ice dryness ?"
161         icedryness=1. ! default value
162         call getin_p("icedryness",icedryness)
163         write(*,*) "surfini: icedryness = ",icedryness
164         dryness (:) = icedryness
165
166      ! To be able to run in parallel, we work on the full grid
167      ! and dispatch results afterwards
168
169      ! start by geting latitudes and logitudes on full grid
170      ! (in serial mode, this is just a copy)
171      call gather(latitude,lati_glo)
172      call gather(longitude,long_glo)
173      call gather(watercaptag,watercaptag_glo)
174
175      if (is_master) then
176
177        IF (ngrid == 1) THEN ! special case for 1d --> do nothing
178
179         write(*,*) 'ngrid = 1, do no put ice caps in surfini.F'
180
181        ELSE
182
183        select case (icelocationmode)
184            case(1) ! icelocationmode == 1
185
186         write(*,*)'Surfini: ice caps defined from surface.nc'
187
188! This method detects ice as gridded value above min_icevalue in the field "string" from surface.nc
189! Typically, it is for thermal inertia above 500 tiu.
190! Two conditions are verified:
191! 1. GCM ice caps are defined such as area is conserved for a given latitude
192! (the approximation is that all points within the GCM latitude resolution have the same area).
193! 2. caps are placed to fill the GCM points with the most detected ice first.
194
195         zedatafile = trim(datadir)
196
197         ierr=nf90_open(trim(zedatafile)//'/surface.nc',NF90_NOWRITE,nid)
198
199         IF (ierr.NE.nf90_noerr) THEN
200       write(*,*)'Error : cannot open file surface.nc '
201       write(*,*)'(in phymars/surfini.F)'
202       write(*,*)'It should be in :',trim(zedatafile),'/'
203       write(*,*)'1) You can set this path in the callphys.def file:'
204       write(*,*)'   datadir=/path/to/the/datafiles'
205       write(*,*)'2) If necessary, surface.nc (and other datafiles)'
206       write(*,*)'   can be obtained online on:'
207       write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
208       call abort_physic("surfini","missing surface.nc file",1)
209         ENDIF
210
211
212         ierr=nf90_inq_varid(nid, string, nvarid)
213         if (ierr.ne.nf90_noerr) then
214          write(*,*) 'surfini error, cannot find ',trim(string)
215          write(*,*) ' in file ',trim(zedatafile),'/surface.nc'
216          write(*,*)trim(nf90_strerror(ierr))
217          call abort_physic("surfini","missing "//trim(string),1)
218         endif
219
220         ierr=nf90_get_var(nid, nvarid, zdata)
221
222         if (ierr.ne.nf90_noerr) then
223          write(*,*) 'surfini: error failed loading ',trim(string)
224          write(*,*)trim(nf90_strerror(ierr))
225          call abort_physic("surfini","failed loading "//trim(string),1)
226         endif
227
228
229         ierr=nf90_close(nid)
230
231
232         nb_ice(:,1) = 1 ! default: there is no ice
233         latice(:,1) = 1
234         lonice(:,1) = 1
235         nb_ice(:,2) = 0
236         latice(:,2) = 0
237         lonice(:,2) = 0
238         !write(*,*)'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi
239
240         ! loop over the GCM grid - except for poles (ig=1 and ngrid)
241         do ig=2,klon_glo-1
242
243          ! loop over the surface file grid
244          do i=1,imd
245           do j=1,jmd
246             zelon = i - 180.
247             zelat = 90. - j
248            if ((abs(lati_glo(ig)*180./pi-zelat) .le. 90./real(nbp_lat-1)) .and. &
249                (abs(long_glo(ig)*180./pi-zelon) .le. 180./real(nbp_lon))) then
250              ! count all points in that GCM grid point
251              nb_ice(ig,1) = nb_ice(ig,1) + 1
252              ! count all detected points in that GCM grid point
253              if (zdata(i,j) > min_icevalue) nb_ice(ig,2) = nb_ice(ig,2) + 1
254             endif
255           enddo
256          enddo
257
258        ! projection of nb_ice on GCM lat and lon axes
259          latice(1+(ig-2)/nbp_lon,:) = latice(1+(ig-2)/nbp_lon,:) + nb_ice(ig,:)
260          lonice(1+mod(ig-2,nbp_lon),:) = lonice(1+mod(ig-2,nbp_lon),:) + nb_ice(ig,:) ! lonice is USELESS ...
261
262         enddo ! of do ig=2,klon_glo-1
263
264
265
266         ! special case for poles
267         nb_ice(1,2)   = 1  ! ice prescribed on north pole
268         latice(1,:)   = nb_ice(1,:)
269         lonice(1,:)   = nb_ice(1,:)
270         latice(nbp_lat-1,:) = nb_ice(ngrid,:)
271         lonice(nbp_lon,:) = nb_ice(ngrid,:)
272
273
274!      write(*,*) 'latice TOT', latice(:,1)
275!      write(*,*) 'latice FOUND', latice(:,2)
276!      write(*,*) 'lonice TOT', lonice(:,1)
277!      write(*,*) 'lonice FOUND', lonice(:,2)
278
279!      write(*,*) 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
280!      write(*,*) 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
281
282!      write(*,*)''
283!      write(*,*)'sum lat', sum(latice(:,1)), sum(lonice(:,1))
284!      write(*,*)'sum lon', sum(latice(:,2)), sum(lonice(:,2))
285
286
287         ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
288         do i=1,(nbp_lat-1)/2
289          step  = 1. ! threshold to add ice cap
290          count = 0. ! number of ice GCM caps at this latitude
291          ! ratiolat is the ratio of area covered by ice within this GCM latitude range
292          ratiolat  = real(latice(i,2))/real(latice(i,1))
293          !write(*,*)'i',i,(i-1)*iim+2,i*iim+1
294
295          ! put ice caps while there is not enough ice,
296          ! as long as the threshold is above 20%
297          do while ((count.le.ratiolat*nbp_lon).and.(step.ge.0.2))
298           count = 0.
299           ! loop over GCM longitudes
300           do j=1,nbp_lon
301            ! if the detected ice ratio in the GCM grid point
302            ! is more than 'step', then add ice
303            if (real(nb_ice((i-1)*nbp_lon+1+j,2)) / real(nb_ice((i-1)*nbp_lon+1+j,1)) .ge. step) then
304                  watercaptag_glo((i-1)*nbp_lon+1+j) = .true.
305                  count = count + 1
306            endif
307           enddo ! of do j=1,nbp_lon
308           !write(*,*) 'step',step,count,ratiolat*nbp_lon
309           step = step - 0.01
310          enddo ! of do while
311          !write(*,*) 'step',step,count,ratiolat*nbp_lon
312
313         enddo ! of do i=1,jjm/2
314
315
316            case(2) ! icelocationmode == 2
317
318         write(*,*)'Surfini: predefined ice caps'
319
320         if ((nbp_lon.eq.32).and.((nbp_lat-1).eq.24)) then ! 32x24
321
322          write(*,*)'water ice caps distribution for 32x24 resolution'
323          longwatercaptag(1:9)    = .true. ! central cap - core
324          longwatercaptag(26:33)  = .true. ! central cap
325          longwatercaptag(1:33)  = .true. ! central cap
326          longwatercaptag(56)  = .true. ! central cap
327          longwatercaptag(58)  = .true. ! central cap
328          longwatercaptag(60)  = .true. ! central cap
329          longwatercaptag(62)  = .true. ! central cap
330          longwatercaptag(64)  = .true. ! central cap
331!---------------------   OUTLIERS  ----------------------------
332
333         else if ((nbp_lon.eq.64).and.((nbp_lat-1).eq.48)) then ! 64x48
334
335          write(*,*)'water ice caps distribution for 64x48 resolution'
336          longwatercaptag(1:65)   = .true. ! central cap - core
337          longwatercaptag(75:85)  = .true. ! central cap
338          longwatercaptag(93:114) = .true. ! central cap
339!---------------------   OUTLIERS  ----------------------------
340          if (.true.) then
341          longwatercaptag(136)    = .true. ! outlier, lat = 78.75
342          longwatercaptag(138)    = .true. ! outlier, lat = 78.75
343          longwatercaptag(140)    = .true. ! outlier, lat = 78.75
344          longwatercaptag(142)    = .true. ! outlier, lat = 78.75
345          longwatercaptag(161)    = .true. ! outlier, lat = 78.75
346          longwatercaptag(163)    = .true. ! outlier, lat = 78.75
347          longwatercaptag(165)    = .true. ! outlier, lat = 78.75
348          longwatercaptag(183)    = .true. ! outlier, lat = 78.75
349          longwatercaptag(185)    = .true. ! outlier, lat = 78.75
350          longwatercaptag(187)    = .true. ! outlier, lat = 78.75
351          longwatercaptag(189)    = .true. ! outlier, lat = 78.75
352          longwatercaptag(191)    = .true. ! outlier, lat = 78.75
353          longwatercaptag(193)    = .true. ! outlier, lat = 78.75
354          longwatercaptag(194)    = .true. ! outlier, lat = 75
355          longwatercaptag(203)    = .true. ! outlier, lat = 75
356          longwatercaptag(207)    = .true. ! outlier, lat = 75
357          longwatercaptag(244)    = .true. ! outlier, lat = 75
358          longwatercaptag(246)    = .true. ! outlier, lat = 75
359          longwatercaptag(250)    = .true. ! outlier, lat = 75
360          longwatercaptag(252)    = .true. ! outlier, lat = 75
361          longwatercaptag(254)    = .true. ! outlier, lat = 75
362          longwatercaptag(256)    = .true. ! outlier, lat = 75
363          endif
364!--------------------------------------------------------------
365
366
367
368         else if (klon_glo .ne. 1) then
369
370          write(*,*)'No predefined ice location for this resolution :',nbp_lon,nbp_lat-1
371          write(*,*)'Please change icelocationmode in surfini.F'
372          write(*,*)'Or add some new definitions ...'
373          call abort_physic("surfini","no pre-definitions for this resolution",1)
374
375         endif
376
377         do ig=1,klon_glo
378          if (longwatercaptag(ig)) watercaptag_glo(ig) = .true.
379         enddo
380
381
382            case(3) ! icelocationmode == 3
383
384         write(*,*)'Surfini: ice caps defined by lat and lon values'
385
386         do ig=1,klon_glo
387
388!-------- Towards olympia planitia water caps -----------
389!--------------------------------------------------------
390
391          if ( ( ( lati_glo(ig)*180./pi .ge. 77.  ) .and.  & ! cap #2
392                 ( lati_glo(ig)*180./pi .le. 80.  ) .and.  &
393                 ( long_glo(ig)*180./pi .ge. 110. ) .and.  &
394                 ( long_glo(ig)*180./pi .le. 181. ) ) .or. &
395               ( ( lati_glo(ig)*180./pi .ge. 75.  ) .and.  & ! cap #4 (Korolev crater)
396                 ( lati_glo(ig)*180./pi .le. 76.  ) .and.  &
397                 ( long_glo(ig)*180./pi .ge. 150. ) .and.  &
398                 ( long_glo(ig)*180./pi .le. 168. ) ) .or. &
399               ( ( lati_glo(ig)*180./pi .ge. 77 ) .and.    & ! cap #5
400                 ( lati_glo(ig)*180./pi .le. 80.  ) .and.  &
401                 ( long_glo(ig)*180./pi .ge. -150.) .and.  &
402                 ( long_glo(ig)*180./pi .le. -110.) ) ) then
403
404               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
405              !    watercaptag(ig)=.true.
406                  alternate = 1
407               else
408                  alternate = 0
409               endif !end if alternate = 0
410
411          endif
412
413!----------- Opposite olympia planitia water cap --------
414!--------------------------------------------------------
415
416          if ( ( ( lati_glo(ig)*180./pi .ge. 80 ) .and.   &
417               ( lati_glo(ig)*180./pi   .le. 84 ) ) .and. &
418             ( ( long_glo(ig)*180./pi .lt. -95. ) .or.    &  !!! 32x24
419               ( long_glo(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
420!           ( ( ( long_glo(ig)*180./pi .ge. -29. ) .and.  &   !!! 64x48
421!               ( long_glo(ig)*180./pi .le.  90. ) ) .or. &   !!! 64x48
422!             ( ( long_glo(ig)*180./pi .ge. -77. ) .and.  &   !!! 64x48
423!               ( long_glo(ig)*180./pi .le. -70. ) ) ) ) then !!! 64x48
424        !   watercaptag_glo(ig)=.true.
425          endif
426
427
428! -------------------- Central cap ----------------------
429!--------------------------------------------------------
430
431          if (abs(lati_glo(ig)*180./pi).gt.80) watercaptag_glo(ig)=.true.
432
433!--------------------------------------------------------
434!--------------------------------------------------------
435         end do ! of (klon_glo)
436
437            case(4) ! icelocationmode == 4
438
439         write(*,*)'icelocationmode = 4'
440         write(*,*)'Surfini: ice caps defined using manual 64x48 settings'
441         write(*,*)'(although, it should work with any resolution)'
442         call locate_watercaptag(klon_glo,lati_glo,long_glo,watercaptag_glo)
443
444!         write(*,*)'watercaptag_glo(:), ',watercaptag_glo(:)
445
446            case(5) ! icelocationmode == 5
447
448         write(*,*)'icelocationmode = 5'
449         write(*,*)'Surfini: ice caps defined using startfi.nc data'
450
451            case default
452
453         write(*,*) 'In surfini.F, icelocationmode is ', icelocationmode
454         write(*,*) 'It should be 1, 2, 3, 4 or 5 (default is 5)'
455         call abort_physic("surfini","wrong icelocationmode",1)
456
457        end select
458
459        ENDIF ! of if (ngrid)
460
461
462        ! print caps locations - useful for plots too
463        write(*,*)'surfini: latitude | longitude | ig'
464        do ig=1,klon_glo
465          dryness_glo(ig) = icedryness
466
467          if (watercaptag_glo(ig)) then
468             write(*,*)'surfini: ice water cap', lati_glo(ig)*180./pi,long_glo(ig)*180./pi, ig
469!             write(1,*),ig, lati_glo(ig)*180./pi,cell_area(ig)
470!             write(2,*), lati_glo(ig)*180./pi,long_glo(ig)*180./pi,cell_area(ig)
471!             write(3,*), ig, lati_glo(ig)*180./pi,long_glo(ig)*180./pi,cell_area(ig)
472          endif
473        enddo
474
475       endif !of if (is_master)
476
477       if (ngrid.gt.1) then
478        ! Now scatter fields watercaptag and dryness from master to all
479        ! (is just a plain copy in serial mode)
480        call scatter(dryness_glo,dryness)
481        call scatter(watercaptag_glo,watercaptag)
482       endif
483       ELSE
484         watercaptag(:) = .false.
485       ENDIF ! water
486! end of #else of #ifndef MESOSCALE
487#endif
488      END SUBROUTINE surfini
489
490!=======================================================================
491
492      SUBROUTINE locate_watercaptag(klon_glo,lati_glo,long_glo,watercaptag_glo)
493
494      USE comcstfi_h, ONLY: pi
495
496      integer, intent(in) :: klon_glo
497      real, intent(in) :: lati_glo(klon_glo)
498      real, intent(in) :: long_glo(klon_glo)
499      logical, intent(out) :: watercaptag_glo(klon_glo)
500      integer :: ig,i
501!      real, dimension(klon_glo,120) :: wcap
502      real, dimension(120,2) :: latedge
503      real, dimension(120,2) :: lonedge
504
505! In icelocationmode=2 there are 120 manually predefined grid points where
506! watercaptag is true (for the 64x48 resolution). The grid cells corners
507! coordinates in latitude and longitude are written below. With this
508! routine, we check if the grid cell center is in between any of those
509! points. If so, watercaptag = true.
510
511      latedge(:,1)=(/                                                 &
512      88.125, 84.375, 84.375, 84.375, 84.375, 84.375,84.375, 84.375,  &
513      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
514      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
515      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
516      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
517      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
518      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
519      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
520      84.375, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, &
521      80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, &
522      80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, &
523      80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, &
524      80.625, 80.625, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, &
525      76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 73.125, &
526      73.125, 73.125, 73.125, 73.125, 73.125, 73.125, 73.125, 73.125/)
527
528      latedge(:,2)=(/                                                 &
529      90. , 88.125, 88.125, 88.125, 88.125, 88.125,88.125, 88.125,    &
530      88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, &
531      88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, &
532      88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, &
533      88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, &
534      88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, &
535      88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, &
536      88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, &
537      88.125, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
538      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
539      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
540      84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, &
541      84.375, 84.375, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, &
542      80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 76.875, &
543      76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875/)
544
545      lonedge(:,1)=(/                                                 &
546     -180.    , -180.    , -177.1875, -171.5625,-165.9375, -160.3125, &
547     -154.6875, -149.0625, -143.4375, -137.8125, -132.1875,-126.5625, &
548     -120.9375, -115.3125, -109.6875, -104.0625,  -98.4375, -92.8125, &
549      -87.1875,  -81.5625,  -75.9375,  -70.3125,  -64.6875, -59.0625, &
550      -53.4375,  -47.8125,  -42.1875,  -36.5625,  -30.9375, -25.3125, &
551      -19.6875,  -14.0625,   -8.4375,   -2.8125,    2.8125,   8.4375, &
552       14.0625,   19.6875,   25.3125,   30.9375,   36.5625,  42.1875, &
553       47.8125,   53.4375,   59.0625,   64.6875,   70.3125,  75.9375, &
554       81.5625,   87.1875,   92.8125,   98.4375,  104.0625, 109.6875, &
555      115.3125,  120.9375,  126.5625,  132.1875,  137.8125, 143.4375, &
556      149.0625,  154.6875,  160.3125,  165.9375,  171.5625,-132.1875, &
557     -126.5625, -120.9375, -115.3125, -109.6875, -104.0625, -98.4375, &
558      -92.8125,  -87.1875,  -81.5625,  -75.9375,  -30.9375, -25.3125, &
559      -19.6875,  -14.0625,   -8.4375,   -2.8125,    2.8125,   8.4375, &
560       14.0625,   19.6875,   25.3125,   30.9375,   36.5625,  42.1875, &
561       47.8125,   53.4375,   59.0625,   64.6875,   70.3125,  75.9375, &
562       81.5625,   87.1875, -149.0625, -137.8125, -126.5625,-115.3125, &
563       -8.4375,    2.8125,   14.0625,  115.3125,  126.5625, 137.8125, &
564      149.0625,  160.3125,  171.5625, -180.    , -132.1875,-109.6875, &
565       98.4375,  109.6875,  132.1875,  143.4375,  154.6875,165.9375/)
566
567      lonedge(:,2)=(/                                                  &
568       180.    , -180.    , -171.5625, -165.9375,-160.3125, -154.6875, &
569      -149.0625,-143.4375, -137.8125, -132.1875, -126.5625, -120.9375, &
570      -115.3125,-109.6875, -104.0625,  -98.4375,  -92.8125,  -87.1875, &
571       -81.5625, -75.9375,  -70.3125,  -64.6875,  -59.0625,  -53.4375, &
572       -47.8125, -42.1875,  -36.5625,  -30.9375,  -25.3125,  -19.6875, &
573       -14.0625,  -8.4375,   -2.8125,    2.8125,    8.4375,   14.0625, &
574        19.6875,  25.3125,   30.9375,   36.5625,   42.1875,   47.8125, &
575        53.4375,  59.0625,   64.6875,   70.3125,   75.9375,   81.5625, &
576        87.1875,  92.8125,   98.4375,  104.0625,  109.6875,  115.3125, &
577       120.9375, 126.5625,  132.1875,  137.8125,  143.4375,  149.0625, &
578       154.6875, 160.3125,  165.9375,  171.5625,  177.1875, -126.5625, &
579      -120.9375,-115.3125, -109.6875, -104.0625,  -98.4375,  -92.8125, &
580       -87.1875, -81.5625,  -75.9375,  -70.3125,  -25.3125,  -19.6875, &
581       -14.0625,  -8.4375,   -2.8125,    2.8125,    8.4375,   14.0625, &
582        19.6875,  25.3125,   30.9375,   36.5625,   42.1875,   47.8125, &
583        53.4375,  59.0625,   64.6875,   70.3125,   75.9375,   81.5625, &
584        87.1875,  92.8125, -143.4375, -132.1875, -120.9375, -109.6875, &
585        -2.8125,   8.4375,   19.6875,  120.9375,  132.1875,  143.4375, &
586       154.6875, 165.9375,  177.1875, -177.1875, -126.5625, -104.0625, &
587       104.0625, 115.3125,  137.8125,  149.0625,  160.3125,171.5625/)
588
589      watercaptag_glo(:) = .false.
590      DO ig=1, klon_glo
591        DO i=1, 120
592           if ((long_glo(ig)*180./pi.ge.lonedge(i,1))      &
593               .and.(long_glo(ig)*180./pi.le.lonedge(i,2)) &
594               .and.(lati_glo(ig)*180./pi.ge.latedge(i,1)) &
595               .and.(lati_glo(ig)*180./pi.le.latedge(i,2))) then
596             watercaptag_glo(ig) = .true.
597           endif
598        ENDDO !i=1, 120
599      ENDDO ! ig=1, klon_glo
600
601      END SUBROUTINE locate_watercaptag
602
603END MODULE surfini_mod
Note: See TracBrowser for help on using the repository browser.