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

Last change on this file since 3026 was 3013, checked in by llange, 17 months ago

MARS PCM
CO2 condens: Fixing bug introduced in r2977: as the cooling of the air before condensing is now done in the surface condensation loop, the boundary conditions for the VL scheme must be changed in order to avoid duplicating the cooling.

Also fixed a bug in surfini.F (2892): qsurf is now ngrid x nq x nslope: I added the dimension nslope.

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