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

Last change on this file since 2596 was 2508, checked in by jnaar, 4 years ago

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