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

Last change on this file since 2887 was 2884, checked in by romain.vande, 3 years ago

Mars PCM:
Watercaptag is now outputed in the starfi.nc.
There is a new mode : icelocationmode=5 to read watercaptag from startfi and to do as iceloctaionmode=4 if the variable is not present.
RV
File size: 24.7 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!           watercaptag(:) = .false.
145           longwatercaptag(:) = .false.
146         endif
147         
148         write(*,*) "surfini: Ice dryness ?"
149         icedryness=1. ! default value
150         call getin_p("icedryness",icedryness)
151         write(*,*) "surfini: icedryness = ",icedryness
152         dryness (:) = icedryness
153         
154      ! To be able to run in parallel, we work on the full grid
155      ! and dispatch results afterwards
156
157      ! start by geting latitudes and logitudes on full grid
158      ! (in serial mode, this is just a copy)
159      call gather(latitude,lati_glo)
160      call gather(longitude,long_glo)
161      call gather(watercaptag,watercaptag_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 IF (icelocationmode .eq. 5) THEN
451     
452         print*,'icelocationmode = 5'
453         print*,'Surfini: ice caps defined using startfi.nc data'
454         do ig=1,ngrid
455           if(any(watercaptag_glo)) then
456           else
457              call locate_watercaptag(klon_glo,lati_glo,
458     &            long_glo,watercaptag_glo)
459           endif
460         enddo
461
462!         print*,'watercaptag_glo(:), ',watercaptag_glo(:)
463
464        ELSE
465     
466         print*, 'In surfini.F, icelocationmode is ', icelocationmode
467         print*, 'It should be 1, 2, 3, 4 or 5 (default is 5)'
468         call abort_physic("surfini","wrong icelocationmode",1)
469
470        ENDIF ! of if (icelocation)
471       
472       
473        ! print caps locations - useful for plots too
474        print*,'surfini: latitude | longitude | ig'
475        do ig=1,klon_glo
476          dryness_glo(ig) = icedryness
477
478          if (watercaptag_glo(ig)) then
479             print*,'surfini: ice water cap', lati_glo(ig)*180./pi,
480     &              long_glo(ig)*180./pi, ig
481!             write(1,*),ig, lati_glo(ig)*180./pi,
482!     &              cell_area(ig)
483!             write(2,*), lati_glo(ig)*180./pi,
484!     &              long_glo(ig)*180./pi,cell_area(ig)
485!             write(3,*), ig, lati_glo(ig)*180./pi,
486!     &              long_glo(ig)*180./pi,cell_area(ig)
487          endif
488        enddo
489       
490       endif !of if (is_master)
491       
492       if (ngrid.gt.1) then
493        ! Now scatter fields watercaptag and dryness from master to all
494        ! (is just a plain copy in serial mode)
495        call scatter(dryness_glo,dryness)
496        call scatter(watercaptag_glo,watercaptag)
497       endif
498       ELSE
499         watercaptag(:) = .false.
500       ENDIF ! water
501! end of #else of #ifndef MESOSCALE
502#endif       
503!      END SUBROUTINE surfini(ngrid,piceco2,qsurf)
504      END !SUBROUTINE surfini(ngrid,piceco2,qsurf)
505
506      SUBROUTINE locate_watercaptag(klon_glo,lati_glo,
507     &            long_glo,watercaptag_glo)
508
509      USE comcstfi_h, ONLY: pi
510
511      integer, intent(in) :: klon_glo
512      real, intent(in) :: lati_glo(klon_glo)
513      real, intent(in) :: long_glo(klon_glo)
514      logical, intent(out) :: watercaptag_glo(klon_glo)
515      integer :: ig,i
516!      real, dimension(klon_glo,120) :: wcap
517      real, dimension(120,2) :: latedge
518      real, dimension(120,2) :: lonedge
519
520! In icelocationmode=2 there are 120 manually predefined grid points where
521! watercaptag is true (for the 64x48 resolution). The grid cells corners
522! coordinates in latitude and longitude are written below. With this
523! routine, we check if the grid cell center is in between any of those
524! points. If so, watercaptag = true.
525
526
527
528
529      latedge(:,1)=(/
530     & 88.125, 84.375, 84.375, 84.375, 84.375, 84.375,84.375, 84.375,
531     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
532     & 84.375, 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, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
539     & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
540     & 80.625, 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, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875,
543     & 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 73.125,
544     & 73.125, 73.125, 73.125, 73.125, 73.125, 73.125, 73.125, 73.125/)
545
546
547      latedge(:,2)=(/
548     & 90. , 88.125, 88.125, 88.125, 88.125, 88.125,88.125, 88.125,   
549     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
550     & 88.125, 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, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
557     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
558     & 84.375, 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, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
561     & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 76.875,
562     & 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875/)
563
564
565      lonedge(:,1)=(/
566     &-180.    , -180.    , -177.1875, -171.5625,-165.9375, -160.3125,
567     &-154.6875, -149.0625, -143.4375, -137.8125, -132.1875,-126.5625,
568     &-120.9375, -115.3125, -109.6875, -104.0625,  -98.4375, -92.8125,
569     & -87.1875,  -81.5625,  -75.9375,  -70.3125,  -64.6875, -59.0625,
570     & -53.4375,  -47.8125,  -42.1875,  -36.5625,  -30.9375, -25.3125,
571     & -19.6875,  -14.0625,   -8.4375,   -2.8125,    2.8125,   8.4375,
572     &  14.0625,   19.6875,   25.3125,   30.9375,   36.5625,  42.1875,
573     &  47.8125,   53.4375,   59.0625,   64.6875,   70.3125,  75.9375,
574     &  81.5625,   87.1875,   92.8125,   98.4375,  104.0625, 109.6875,
575     & 115.3125,  120.9375,  126.5625,  132.1875,  137.8125, 143.4375,
576     & 149.0625,  154.6875,  160.3125,  165.9375,  171.5625,-132.1875,
577     &-126.5625, -120.9375, -115.3125, -109.6875, -104.0625, -98.4375,
578     & -92.8125,  -87.1875,  -81.5625,  -75.9375,  -30.9375, -25.3125,
579     & -19.6875,  -14.0625,   -8.4375,   -2.8125,    2.8125,   8.4375,
580     &  14.0625,   19.6875,   25.3125,   30.9375,   36.5625,  42.1875,
581     &  47.8125,   53.4375,   59.0625,   64.6875,   70.3125,  75.9375,
582     &  81.5625,   87.1875, -149.0625, -137.8125, -126.5625,-115.3125,
583     &  -8.4375,    2.8125,   14.0625,  115.3125,  126.5625, 137.8125,
584     & 149.0625,  160.3125,  171.5625, -180.    , -132.1875,-109.6875,
585     &  98.4375,  109.6875,  132.1875,  143.4375,  154.6875,165.9375/)
586
587      lonedge(:,2)=(/
588     & 180.    , -180.    , -171.5625, -165.9375,-160.3125, -154.6875,
589     &-149.0625,-143.4375, -137.8125, -132.1875, -126.5625, -120.9375,
590     &-115.3125,-109.6875, -104.0625,  -98.4375,  -92.8125,  -87.1875,
591     & -81.5625, -75.9375,  -70.3125,  -64.6875,  -59.0625,  -53.4375,
592     & -47.8125, -42.1875,  -36.5625,  -30.9375,  -25.3125,  -19.6875,
593     & -14.0625,  -8.4375,   -2.8125,    2.8125,    8.4375,   14.0625,
594     &  19.6875,  25.3125,   30.9375,   36.5625,   42.1875,   47.8125,
595     &  53.4375,  59.0625,   64.6875,   70.3125,   75.9375,   81.5625,
596     &  87.1875,  92.8125,   98.4375,  104.0625,  109.6875,  115.3125,
597     & 120.9375, 126.5625,  132.1875,  137.8125,  143.4375,  149.0625,
598     & 154.6875, 160.3125,  165.9375,  171.5625,  177.1875, -126.5625,
599     &-120.9375,-115.3125, -109.6875, -104.0625,  -98.4375,  -92.8125,
600     & -87.1875, -81.5625,  -75.9375,  -70.3125,  -25.3125,  -19.6875,
601     & -14.0625,  -8.4375,   -2.8125,    2.8125,    8.4375,   14.0625,
602     &  19.6875,  25.3125,   30.9375,   36.5625,   42.1875,   47.8125,
603     &  53.4375,  59.0625,   64.6875,   70.3125,   75.9375,   81.5625,
604     &  87.1875,  92.8125, -143.4375, -132.1875, -120.9375, -109.6875,
605     &  -2.8125,   8.4375,   19.6875,  120.9375,  132.1875,  143.4375,
606     & 154.6875, 165.9375,  177.1875, -177.1875, -126.5625, -104.0625,
607     & 104.0625, 115.3125,  137.8125,  149.0625,  160.3125,171.5625/)
608
609
610      watercaptag_glo(:) = .false.
611      DO ig=1, klon_glo
612        DO i=1, 120
613           if ((long_glo(ig)*180./pi.ge.lonedge(i,1))
614     &         .and.(long_glo(ig)*180./pi.le.lonedge(i,2))
615     &         .and.(lati_glo(ig)*180./pi.ge.latedge(i,1))
616     &         .and.(lati_glo(ig)*180./pi.le.latedge(i,2))) then
617             watercaptag_glo(ig) = .true.
618           endif
619        ENDDO !i=1, 120
620      ENDDO ! ig=1, klon_glo
621
622      END SUBROUTINE locate_watercaptag
Note: See TracBrowser for help on using the repository browser.