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

Last change on this file since 2800 was 2741, checked in by jnaar, 2 years ago

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