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

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

Mars GCM:
The variable co2ice is deleted. All the co2 ice on surface is now in qsurf(:,igcm_co2).
CO2 tracer is now mandatory. diagfi output is unchanged.
RV
File size: 24.2 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 = 4
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.
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
162      if (is_master) then
163
164        IF (ngrid .eq. 1) THEN ! special case for 1d --> do nothing
165     
166         print*, 'ngrid = 1, do no put ice caps in surfini.F'
167
168        ELSE IF (icelocationmode .eq. 1) THEN
169     
170         print*,'Surfini: ice caps defined from surface.nc'
171           
172! This method detects ice as gridded value above min_icevalue in the field "string" from surface.nc
173! Typically, it is for thermal inertia above 500 tiu.
174! Two conditions are verified:
175! 1. GCM ice caps are defined such as area is conserved for a given latitude
176! (the approximation is that all points within the GCM latitude resolution have the same area).
177! 2. caps are placed to fill the GCM points with the most detected ice first.
178     
179
180           
181         zedatafile = trim(datadir)
182 
183       
184         ierr=nf90_open(trim(zedatafile)//'/surface.nc',
185     &   NF90_NOWRITE,nid)
186     
187         IF (ierr.NE.nf90_noerr) THEN
188       write(*,*)'Error : cannot open file surface.nc '
189       write(*,*)'(in phymars/surfini.F)'
190       write(*,*)'It should be in :',trim(zedatafile),'/'
191       write(*,*)'1) You can set this path in the callphys.def file:'
192       write(*,*)'   datadir=/path/to/the/datafiles'
193       write(*,*)'2) If necessary, surface.nc (and other datafiles)'
194       write(*,*)'   can be obtained online on:'
195       write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
196       call abort_physic("surfini","missing surface.nc file",1)
197         ENDIF
198     
199     
200         ierr=nf90_inq_varid(nid, string, nvarid)
201         if (ierr.ne.nf90_noerr) then
202          write(*,*) 'surfini error, cannot find ',trim(string)
203          write(*,*) ' in file ',trim(zedatafile),'/surface.nc'
204          write(*,*)trim(nf90_strerror(ierr))
205          call abort_physic("surfini","missing "//trim(string),1)
206         endif
207
208         ierr=nf90_get_var(nid, nvarid, zdata)
209
210         if (ierr.ne.nf90_noerr) then
211          write(*,*) 'surfini: error failed loading ',trim(string)
212          write(*,*)trim(nf90_strerror(ierr))
213          call abort_physic("surfini","failed loading "//trim(string),1)
214         endif
215 
216                     
217         ierr=nf90_close(nid)
218 
219
220         nb_ice(:,1) = 1 ! default: there is no ice
221         latice(:,1) = 1
222         lonice(:,1) = 1
223         nb_ice(:,2) = 0
224         latice(:,2) = 0
225         lonice(:,2) = 0
226         !print*,'jjm,iim',jjm,iim ! jjm =  nb lati , iim = nb longi
227
228         ! loop over the GCM grid - except for poles (ig=1 and ngrid)
229         do ig=2,klon_glo-1
230     
231          ! loop over the surface file grid     
232          do i=1,imd
233           do j=1,jmd
234             zelon = i - 180.
235             zelat = 90. - j
236            if ((abs(lati_glo(ig)*180./pi-zelat).le.
237     &           90./real(nbp_lat-1)) .and.
238     &          (abs(long_glo(ig)*180./pi-zelon).le.
239     &           180./real(nbp_lon))) then
240              ! count all points in that GCM grid point
241              nb_ice(ig,1) = nb_ice(ig,1) + 1
242              if (zdata(i,j) > min_icevalue)
243                 ! count all detected points in that GCM grid point
244     &           nb_ice(ig,2) = nb_ice(ig,2) + 1
245             endif
246           enddo
247          enddo 
248
249        ! projection of nb_ice on GCM lat and lon axes
250          latice(1+(ig-2)/nbp_lon,:) =
251     &     latice(1+(ig-2)/nbp_lon,:) + nb_ice(ig,:)
252          lonice(1+mod(ig-2,nbp_lon),:) =
253     &     lonice(1+mod(ig-2,nbp_lon),:) + nb_ice(ig,:) ! lonice is USELESS ...
254
255         enddo ! of do ig=2,klon_glo-1
256     
257
258     
259         ! special case for poles
260         nb_ice(1,2)   = 1  ! ice prescribed on north pole
261         latice(1,:)   = nb_ice(1,:)
262         lonice(1,:)   = nb_ice(1,:)
263         latice(nbp_lat-1,:) = nb_ice(ngrid,:)
264         lonice(nbp_lon,:) = nb_ice(ngrid,:)
265     
266     
267!      print*, 'latice TOT', latice(:,1)
268!      print*, 'latice FOUND', latice(:,2)
269!      print*, 'lonice TOT', lonice(:,1)
270!      print*, 'lonice FOUND', lonice(:,2)
271     
272!      print*, 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
273!      print*, 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
274     
275!      print*,''
276!      print*,'sum lat', sum(latice(:,1)), sum(lonice(:,1))
277!      print*,'sum lon', sum(latice(:,2)), sum(lonice(:,2))
278     
279   
280         ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
281         do i=1,(nbp_lat-1)/2
282          step  = 1. ! threshold to add ice cap
283          count = 0. ! number of ice GCM caps at this latitude
284          ! ratiolat is the ratio of area covered by ice within this GCM latitude range
285          ratiolat  = real(latice(i,2))/real(latice(i,1))
286          !print*,'i',i,(i-1)*iim+2,i*iim+1
287     
288          ! put ice caps while there is not enough ice,
289          ! as long as the threshold is above 20%
290          do while ((count.le.ratiolat*nbp_lon).and.(step.ge.0.2))
291           count = 0.
292           ! loop over GCM longitudes
293           do j=1,nbp_lon
294            ! if the detected ice ratio in the GCM grid point
295            ! is more than 'step', then add ice
296            if (real(nb_ice((i-1)*nbp_lon+1+j,2))
297     &        / real(nb_ice((i-1)*nbp_lon+1+j,1)) .ge. step) then
298                  watercaptag_glo((i-1)*nbp_lon+1+j) = .true.
299                  count = count + 1
300            endif
301           enddo ! of do j=1,nbp_lon
302           !print*, 'step',step,count,ratiolat*nbp_lon
303           step = step - 0.01
304          enddo ! of do while
305          !print*, 'step',step,count,ratiolat*nbp_lon
306
307         enddo ! of do i=1,jjm/2
308           
309
310        ELSE IF (icelocationmode .eq. 2) THEN
311     
312         print*,'Surfini: predefined ice caps'
313     
314         if ((nbp_lon.eq.32).and.((nbp_lat-1).eq.24)) then ! 32x24
315           
316          print*,'water ice caps distribution for 32x24 resolution'
317          longwatercaptag(1:9)    = .true. ! central cap - core
318          longwatercaptag(26:33)  = .true. ! central cap
319          longwatercaptag(1:33)  = .true. ! central cap
320          longwatercaptag(56)  = .true. ! central cap
321          longwatercaptag(58)  = .true. ! central cap
322          longwatercaptag(60)  = .true. ! central cap
323          longwatercaptag(62)  = .true. ! central cap
324          longwatercaptag(64)  = .true. ! central cap
325!---------------------   OUTLIERS  ----------------------------
326
327         else if ((nbp_lon.eq.64).and.((nbp_lat-1).eq.48)) then ! 64x48
328
329          print*,'water ice caps distribution for 64x48 resolution'
330          longwatercaptag(1:65)   = .true. ! central cap - core
331          longwatercaptag(75:85)  = .true. ! central cap
332          longwatercaptag(93:114) = .true. ! central cap
333!---------------------   OUTLIERS  ----------------------------
334          if (.true.) then
335          longwatercaptag(136)    = .true. ! outlier, lat = 78.75
336          longwatercaptag(138)    = .true. ! outlier, lat = 78.75
337          longwatercaptag(140)    = .true. ! outlier, lat = 78.75
338          longwatercaptag(142)    = .true. ! outlier, lat = 78.75
339          longwatercaptag(161)    = .true. ! outlier, lat = 78.75
340          longwatercaptag(163)    = .true. ! outlier, lat = 78.75
341          longwatercaptag(165)    = .true. ! outlier, lat = 78.75
342          longwatercaptag(183)    = .true. ! outlier, lat = 78.75
343          longwatercaptag(185)    = .true. ! outlier, lat = 78.75
344          longwatercaptag(187)    = .true. ! outlier, lat = 78.75
345          longwatercaptag(189)    = .true. ! outlier, lat = 78.75
346          longwatercaptag(191)    = .true. ! outlier, lat = 78.75
347          longwatercaptag(193)    = .true. ! outlier, lat = 78.75
348          longwatercaptag(194)    = .true. ! outlier, lat = 75
349          longwatercaptag(203)    = .true. ! outlier, lat = 75
350          longwatercaptag(207)    = .true. ! outlier, lat = 75
351          longwatercaptag(244)    = .true. ! outlier, lat = 75
352          longwatercaptag(246)    = .true. ! outlier, lat = 75
353          longwatercaptag(250)    = .true. ! outlier, lat = 75
354          longwatercaptag(252)    = .true. ! outlier, lat = 75
355          longwatercaptag(254)    = .true. ! outlier, lat = 75
356          longwatercaptag(256)    = .true. ! outlier, lat = 75
357          endif
358!--------------------------------------------------------------       
359
360
361           
362         else if (klon_glo .ne. 1) then
363       
364          print*,'No predefined ice location for this resolution :',
365     &           nbp_lon,nbp_lat-1
366          print*,'Please change icelocationmode in surfini.F'
367          print*,'Or add some new definitions ...'
368          call abort_physic("surfini",
369     &         "no pre-definitions for this resolution",1)
370         
371         endif
372
373         do ig=1,klon_glo
374          if (longwatercaptag(ig)) watercaptag_glo(ig) = .true.
375         enddo
376
377
378        ELSE IF (icelocationmode .eq. 3) THEN
379     
380         print*,'Surfini: ice caps defined by lat and lon values'
381
382         do ig=1,klon_glo
383         
384c-------- Towards olympia planitia water caps -----------
385c--------------------------------------------------------
386
387          if ( ( ( lati_glo(ig)*180./pi .ge. 77.  ) .and. ! cap #2
388     .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
389     .           ( long_glo(ig)*180./pi .ge. 110. ) .and.
390     .           ( long_glo(ig)*180./pi .le. 181. ) )
391     .         .or.
392
393     .         ( ( lati_glo(ig)*180./pi .ge. 75.  ) .and. ! cap #4 (Korolev crater)
394     .           ( lati_glo(ig)*180./pi .le. 76.  ) .and.
395     .           ( long_glo(ig)*180./pi .ge. 150. ) .and.
396     .           ( long_glo(ig)*180./pi .le. 168. ) )
397     .         .or.
398     .         ( ( lati_glo(ig)*180./pi .ge. 77 ) .and. ! cap #5
399     .           ( lati_glo(ig)*180./pi .le. 80.  ) .and.
400     .           ( long_glo(ig)*180./pi .ge. -150.) .and.
401     .           ( long_glo(ig)*180./pi .le. -110.) ) )
402     .         then
403             
404               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
405              !    watercaptag(ig)=.true.
406                  alternate = 1
407               else
408                  alternate = 0
409               endif !end if alternate = 0
410               
411          endif
412
413c----------- Opposite olympia planitia water cap --------
414c--------------------------------------------------------
415
416          if ( ( ( lati_glo(ig)*180./pi     .ge.  80 ) .and.
417     .         ( lati_glo(ig)*180./pi     .le.  84 ) )
418     .         .and.
419     .       ( ( long_glo(ig)*180./pi .lt. -95. ) .or.       !!! 32x24
420     .         ( long_glo(ig)*180./pi .gt.  85. ) ) ) then   !!! 32x24
421!     .     ( ( ( long_glo(ig)*180./pi .ge. -29. ) .and.       !!! 64x48
422!     .         ( long_glo(ig)*180./pi .le.  90. ) ) .or.      !!! 64x48
423!     .       ( ( long_glo(ig)*180./pi .ge. -77. ) .and.       !!! 64x48
424!     .         ( long_glo(ig)*180./pi .le. -70. ) ) ) ) then  !!! 64x48
425        !   watercaptag_glo(ig)=.true.
426          endif
427
428
429c -------------------- Central cap ----------------------
430c--------------------------------------------------------
431
432          if (abs(lati_glo(ig)*180./pi).gt.80)
433     .          watercaptag_glo(ig)=.true.
434           
435c--------------------------------------------------------
436c--------------------------------------------------------
437         end do ! of (klon_glo)
438
439        ELSE IF (icelocationmode .eq. 4) THEN
440     
441         print*,'icelocationmode = 4'
442         print*,'Surfini: ice caps defined using manual 64x48 settings'
443         print*,'(although, it should work with any resolution)'
444         call locate_watercaptag(klon_glo,lati_glo,
445     &            long_glo,watercaptag_glo)
446
447!         print*,'watercaptag_glo(:), ',watercaptag_glo(:)
448
449        ELSE
450     
451         print*, 'In surfini.F, icelocationmode is ', icelocationmode
452         print*, 'It should be 1, 2, 3 or 4 (default is 4)'
453         call abort_physic("surfini","wrong icelocationmode",1)
454
455        ENDIF ! of if (icelocation)
456       
457       
458        ! print caps locations - useful for plots too
459        print*,'surfini: latitude | longitude | ig'
460        do ig=1,klon_glo
461          dryness_glo(ig) = icedryness
462
463          if (watercaptag_glo(ig)) then
464             print*,'surfini: ice water cap', lati_glo(ig)*180./pi,
465     &              long_glo(ig)*180./pi, ig
466!             write(1,*),ig, lati_glo(ig)*180./pi,
467!     &              cell_area(ig)
468!             write(2,*), lati_glo(ig)*180./pi,
469!     &              long_glo(ig)*180./pi,cell_area(ig)
470!             write(3,*), ig, lati_glo(ig)*180./pi,
471!     &              long_glo(ig)*180./pi,cell_area(ig)
472          endif
473        enddo
474       
475       endif !of if (is_master)
476       
477       if (ngrid.gt.1) then
478        ! Now scatter fields watercaptag and dryness from master to all
479        ! (is just a plain copy in serial mode)
480        call scatter(dryness_glo,dryness)
481        call scatter(watercaptag_glo,watercaptag)
482       endif
483       ELSE
484         watercaptag(:) = .false.
485       ENDIF ! water
486! end of #else of #ifndef MESOSCALE
487#endif       
488!      END SUBROUTINE surfini(ngrid,piceco2,qsurf)
489      END !SUBROUTINE surfini(ngrid,piceco2,qsurf)
490
491      SUBROUTINE locate_watercaptag(klon_glo,lati_glo,
492     &            long_glo,watercaptag_glo)
493
494      USE comcstfi_h, ONLY: pi
495
496      integer, intent(in) :: klon_glo
497      real, intent(in) :: lati_glo(klon_glo)
498      real, intent(in) :: long_glo(klon_glo)
499      logical, intent(out) :: watercaptag_glo(klon_glo)
500      integer :: ig,i
501!      real, dimension(klon_glo,120) :: wcap
502      real, dimension(120,2) :: latedge
503      real, dimension(120,2) :: lonedge
504
505! In icelocationmode=2 there are 120 manually predefined grid points where
506! watercaptag is true (for the 64x48 resolution). The grid cells corners
507! coordinates in latitude and longitude are written below. With this
508! routine, we check if the grid cell center is in between any of those
509! points. If so, watercaptag = true.
510
511
512
513
514      latedge(:,1)=(/
515     & 88.125, 84.375, 84.375, 84.375, 84.375, 84.375,84.375, 84.375,
516     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
517     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
518     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
519     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
520     & 84.375, 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, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
524     & 80.625, 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, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875,
528     & 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 73.125,
529     & 73.125, 73.125, 73.125, 73.125, 73.125, 73.125, 73.125, 73.125/)
530
531
532      latedge(:,2)=(/
533     & 90. , 88.125, 88.125, 88.125, 88.125, 88.125,88.125, 88.125,   
534     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
535     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
536     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
537     & 88.125, 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, 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, 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, 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, 76.875,
547     & 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875/)
548
549
550      lonedge(:,1)=(/
551     &-180.    , -180.    , -177.1875, -171.5625,-165.9375, -160.3125,
552     &-154.6875, -149.0625, -143.4375, -137.8125, -132.1875,-126.5625,
553     &-120.9375, -115.3125, -109.6875, -104.0625,  -98.4375, -92.8125,
554     & -87.1875,  -81.5625,  -75.9375,  -70.3125,  -64.6875, -59.0625,
555     & -53.4375,  -47.8125,  -42.1875,  -36.5625,  -30.9375, -25.3125,
556     & -19.6875,  -14.0625,   -8.4375,   -2.8125,    2.8125,   8.4375,
557     &  14.0625,   19.6875,   25.3125,   30.9375,   36.5625,  42.1875,
558     &  47.8125,   53.4375,   59.0625,   64.6875,   70.3125,  75.9375,
559     &  81.5625,   87.1875,   92.8125,   98.4375,  104.0625, 109.6875,
560     & 115.3125,  120.9375,  126.5625,  132.1875,  137.8125, 143.4375,
561     & 149.0625,  154.6875,  160.3125,  165.9375,  171.5625,-132.1875,
562     &-126.5625, -120.9375, -115.3125, -109.6875, -104.0625, -98.4375,
563     & -92.8125,  -87.1875,  -81.5625,  -75.9375,  -30.9375, -25.3125,
564     & -19.6875,  -14.0625,   -8.4375,   -2.8125,    2.8125,   8.4375,
565     &  14.0625,   19.6875,   25.3125,   30.9375,   36.5625,  42.1875,
566     &  47.8125,   53.4375,   59.0625,   64.6875,   70.3125,  75.9375,
567     &  81.5625,   87.1875, -149.0625, -137.8125, -126.5625,-115.3125,
568     &  -8.4375,    2.8125,   14.0625,  115.3125,  126.5625, 137.8125,
569     & 149.0625,  160.3125,  171.5625, -180.    , -132.1875,-109.6875,
570     &  98.4375,  109.6875,  132.1875,  143.4375,  154.6875,165.9375/)
571
572      lonedge(:,2)=(/
573     & 180.    , -180.    , -171.5625, -165.9375,-160.3125, -154.6875,
574     &-149.0625,-143.4375, -137.8125, -132.1875, -126.5625, -120.9375,
575     &-115.3125,-109.6875, -104.0625,  -98.4375,  -92.8125,  -87.1875,
576     & -81.5625, -75.9375,  -70.3125,  -64.6875,  -59.0625,  -53.4375,
577     & -47.8125, -42.1875,  -36.5625,  -30.9375,  -25.3125,  -19.6875,
578     & -14.0625,  -8.4375,   -2.8125,    2.8125,    8.4375,   14.0625,
579     &  19.6875,  25.3125,   30.9375,   36.5625,   42.1875,   47.8125,
580     &  53.4375,  59.0625,   64.6875,   70.3125,   75.9375,   81.5625,
581     &  87.1875,  92.8125,   98.4375,  104.0625,  109.6875,  115.3125,
582     & 120.9375, 126.5625,  132.1875,  137.8125,  143.4375,  149.0625,
583     & 154.6875, 160.3125,  165.9375,  171.5625,  177.1875, -126.5625,
584     &-120.9375,-115.3125, -109.6875, -104.0625,  -98.4375,  -92.8125,
585     & -87.1875, -81.5625,  -75.9375,  -70.3125,  -25.3125,  -19.6875,
586     & -14.0625,  -8.4375,   -2.8125,    2.8125,    8.4375,   14.0625,
587     &  19.6875,  25.3125,   30.9375,   36.5625,   42.1875,   47.8125,
588     &  53.4375,  59.0625,   64.6875,   70.3125,   75.9375,   81.5625,
589     &  87.1875,  92.8125, -143.4375, -132.1875, -120.9375, -109.6875,
590     &  -2.8125,   8.4375,   19.6875,  120.9375,  132.1875,  143.4375,
591     & 154.6875, 165.9375,  177.1875, -177.1875, -126.5625, -104.0625,
592     & 104.0625, 115.3125,  137.8125,  149.0625,  160.3125,171.5625/)
593
594
595      watercaptag_glo(:) = .false.
596      DO ig=1, klon_glo
597        DO i=1, 120
598           if ((long_glo(ig)*180./pi.ge.lonedge(i,1))
599     &         .and.(long_glo(ig)*180./pi.le.lonedge(i,2))
600     &         .and.(lati_glo(ig)*180./pi.ge.latedge(i,1))
601     &         .and.(lati_glo(ig)*180./pi.le.latedge(i,2))) then
602             watercaptag_glo(ig) = .true.
603           endif
604        ENDDO !i=1, 120
605      ENDDO ! ig=1, klon_glo
606
607      END SUBROUTINE locate_watercaptag
Note: See TracBrowser for help on using the repository browser.