source: trunk/LMDZ.MARS/libf/phymars/surfini.F @ 2947

Last change on this file since 2947 was 2892, checked in by romain.vande, 23 months ago

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