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

Last change on this file since 2814 was 2741, checked in by jnaar, 3 years ago

Initialization of watercaptag no longer checks "caps" boolean option in callphys.def. JN
File size: 24.3 KB
RevLine 
[1944]1      SUBROUTINE surfini(ngrid,piceco2,qsurf)
[1918]2
[2311]3      USE ioipsl_getin_p_mod, ONLY : getin_p
[697]4      use netcdf
[1224]5      use tracer_mod, only: nqmx, noms
[2502]6      use geometry_mod, only: longitude, latitude, ! in radians
7     &                     cell_area ! for watercaptag diagnosis
[1047]8      use surfdat_h, only: watercaptag, frost_albedo_threshold,
[2508]9     &                     albedo_h2o_cap, inert_h2o_ice, albedodat,
[1224]10     &                     albedice, dryness
[1212]11#ifndef MESOSCALE
[1130]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
[1212]14#endif
[1944]15      USE comcstfi_h, ONLY: pi
[1528]16      use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
[1918]17      use datafile_mod, only: datadir
[38]18      IMPLICIT NONE
19c=======================================================================
20c
21c   creation des calottes pour l'etat initial
22c
23c=======================================================================
24c-----------------------------------------------------------------------
25c   Declarations:
26c   -------------
[1918]27      include "callkeys.h"
[669]28
[1130]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
[520]34      REAL icedryness ! ice dryness
[633]35     
36      ! longwatercaptag is watercaptag. Trick for some compilers
37      LOGICAL, DIMENSION(100000) :: longwatercaptag
[520]38     
[2502]39! There are 4 different modes for ice distribution:
[669]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
[2502]43! icelocationmode = 4 ---> predefined 64x48 but usable with every
44! resolution, and easily adaptable for dynamico
[669]45! For visualisation : > /u/tnalmd/bin/watercaps gcm_txt_output_file
[2502]46      INTEGER,SAVE :: icelocationmode = 4
[2616]47     
48!$OMP THREADPRIVATE(icelocationmode)
[669]49       
50       
51      !in case icelocationmode == 1
52      INTEGER i,j
53      INTEGER     imd,jmd
54      PARAMETER   (imd=360,jmd=180)
[712]55      REAL        zdata(imd,jmd)
[669]56      REAL        zelat,zelon
57
[1212]58#ifndef MESOSCALE
[1207]59      INTEGER nb_ice(klon_glo,2)   ! number of counts | detected ice for GCM grid
[1212]60#endif
[1528]61      INTEGER latice(nbp_lat-1,2),lonice (nbp_lon,2) ! number of counts | detected ice along lat & lon axis
[669]62
63      REAL step,count,ratiolat
64
65      INTEGER   ierr,nid,nvarid
66     
67      REAL,SAVE :: min_icevalue = 500.
[2616]68     
69!$OMP THREADPRIVATE(min_icevalue)
70     
[697]71      character(len=50) :: string = 'thermal'
[669]72     
73      character (len=100) :: zedatafile
[1130]74
[1212]75#ifdef MESOSCALE
76
77      do ig=1,ngrid
78
79         !write(*,*) "all qsurf to zero. dirty."
80         do iq=1,nqmx
81         qsurf(ig,iq)=0.  !! on jette les inputs GCM
82                          !! on regle juste watercaptag
83                          !! il faudrait garder les inputs GCM
84                          !! si elles sont consequentes
85         enddo
[1541]86         if ( ( latitude(ig)*180./pi .gt. 70. ) .and.
[1212]87     .        ( albedodat(ig) .ge. 0.26   ) )  then
88                 write(*,*)"outlier ",ig
89                 watercaptag(ig)  = .true.
90                 dryness(ig)      = 1.
[2508]91                 albedodat(ig)    = albedo_h2o_cap  !! pour output
[1212]92         else
93                 watercaptag(ig)  = .false.
94                 dryness(ig)      = 1.
95         endif
96
97      enddo
98#endif
99! problem with nested precompiling flags
100
101#ifndef MESOSCALE
[1130]102      ! to handle parallel cases
103#if CPP_PARA
104      logical watercaptag_glo(klon_glo)
105      real dryness_glo(klon_glo)
106      real lati_glo(klon_glo)
107      real long_glo(klon_glo)
108#else
109      logical watercaptag_glo(ngrid)
110      real dryness_glo(ngrid)
111      real lati_glo(ngrid)
112      real long_glo(ngrid)
113#endif
[1212]114#endif
[1130]115
[1212]116#ifndef MESOSCALE
117
[38]118c
119c=======================================================================
[1130]120! Initialize watercaptag (default is false)
121      watercaptag_glo(:)=.false.
[38]122
[283]123c     water ice outliers
[38]124c     ------------------------------------------
125
[2741]126      IF (water) THEN
[285]127     
[283]128c Perennial H20 north cap defined by watercaptag=true (allows surface to be
129c hollowed by sublimation in vdifc).
[38]130
[283]131c We might not want albedodat to be modified because it is used to write
132c restart files. Instead, albedo is directly modified when needed (i.e.
133c if we have watercaptag and no co2 ice), below and in albedocaps.F90
134
135c       "Dryness coefficient" controlling the evaporation and
136c        sublimation from the ground water ice (close to 1)
137c        HERE, the goal is to correct for the fact
138c        that the simulated permanent water ice polar caps
139c        is larger than the actual cap and the atmospheric
140c        opacity not always realistic.
141
142         alternate = 0
[520]143         
[1047]144         if (ngrid .ne. 1) then
[633]145           watercaptag(:) = .false.
146           longwatercaptag(:) = .false.
147         endif
148         
[1130]149         write(*,*) "surfini: Ice dryness ?"
[520]150         icedryness=1. ! default value
[2311]151         call getin_p("icedryness",icedryness)
[1130]152         write(*,*) "surfini: icedryness = ",icedryness
[669]153         dryness (:) = icedryness
[633]154         
[1130]155      ! To be able to run in parallel, we work on the full grid
156      ! and dispatch results afterwards
[669]157
[1130]158      ! start by geting latitudes and logitudes on full grid
159      ! (in serial mode, this is just a copy)
[1541]160      call gather(latitude,lati_glo)
161      call gather(longitude,long_glo)
[669]162
[1130]163      if (is_master) then
164
165        IF (ngrid .eq. 1) THEN ! special case for 1d --> do nothing
[669]166     
[1047]167         print*, 'ngrid = 1, do no put ice caps in surfini.F'
[669]168
[1130]169        ELSE IF (icelocationmode .eq. 1) THEN
[669]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           
[1918]182         zedatafile = trim(datadir)
[669]183 
184       
[1130]185         ierr=nf90_open(trim(zedatafile)//'/surface.nc',
186     &   NF90_NOWRITE,nid)
[669]187     
[1130]188         IF (ierr.NE.nf90_noerr) THEN
[669]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:'
[1381]196       write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
[2311]197       call abort_physic("surfini","missing surface.nc file",1)
[1130]198         ENDIF
[669]199     
200     
[1130]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))
[2311]206          call abort_physic("surfini","missing "//trim(string),1)
[1130]207         endif
[697]208
[1130]209         ierr=nf90_get_var(nid, nvarid, zdata)
[697]210
[1130]211         if (ierr.ne.nf90_noerr) then
212          write(*,*) 'surfini: error failed loading ',trim(string)
213          write(*,*)trim(nf90_strerror(ierr))
[2311]214          call abort_physic("surfini","failed loading "//trim(string),1)
[1130]215         endif
[669]216 
217                     
[1130]218         ierr=nf90_close(nid)
[669]219 
[285]220
[1130]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
[669]228
[1130]229         ! loop over the GCM grid - except for poles (ig=1 and ngrid)
230         do ig=2,klon_glo-1
[520]231     
[1130]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
[1528]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
[669]241              ! count all points in that GCM grid point
242              nb_ice(ig,1) = nb_ice(ig,1) + 1
[712]243              if (zdata(i,j) > min_icevalue)
[669]244                 ! count all detected points in that GCM grid point
245     &           nb_ice(ig,2) = nb_ice(ig,2) + 1
[712]246             endif
[1130]247           enddo
248          enddo 
[712]249
[669]250        ! projection of nb_ice on GCM lat and lon axes
[1528]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 ...
[669]255
[1130]256         enddo ! of do ig=2,klon_glo-1
[669]257     
258
259     
[1130]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,:)
[1528]264         latice(nbp_lat-1,:) = nb_ice(ngrid,:)
265         lonice(nbp_lon,:) = nb_ice(ngrid,:)
[669]266     
267     
[712]268!      print*, 'latice TOT', latice(:,1)
269!      print*, 'latice FOUND', latice(:,2)
270!      print*, 'lonice TOT', lonice(:,1)
271!      print*, 'lonice FOUND', lonice(:,2)
[669]272     
[712]273!      print*, 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
274!      print*, 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
[669]275     
[712]276!      print*,''
277!      print*,'sum lat', sum(latice(:,1)), sum(lonice(:,1))
278!      print*,'sum lon', sum(latice(:,2)), sum(lonice(:,2))
[669]279     
280   
[1130]281         ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
[1528]282         do i=1,(nbp_lat-1)/2
[1130]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
[669]288     
[1130]289          ! put ice caps while there is not enough ice,
290          ! as long as the threshold is above 20%
[1528]291          do while ((count.le.ratiolat*nbp_lon).and.(step.ge.0.2))
[1130]292           count = 0.
293           ! loop over GCM longitudes
[1528]294           do j=1,nbp_lon
[669]295            ! if the detected ice ratio in the GCM grid point
296            ! is more than 'step', then add ice
[1528]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.
[669]300                  count = count + 1
301            endif
[1528]302           enddo ! of do j=1,nbp_lon
303           !print*, 'step',step,count,ratiolat*nbp_lon
[1130]304           step = step - 0.01
305          enddo ! of do while
[1528]306          !print*, 'step',step,count,ratiolat*nbp_lon
[669]307
[1130]308         enddo ! of do i=1,jjm/2
[669]309           
310
[1130]311        ELSE IF (icelocationmode .eq. 2) THEN
[669]312     
[1130]313         print*,'Surfini: predefined ice caps'
[669]314     
[1528]315         if ((nbp_lon.eq.32).and.((nbp_lat-1).eq.24)) then ! 32x24
[633]316           
[669]317          print*,'water ice caps distribution for 32x24 resolution'
[633]318          longwatercaptag(1:9)    = .true. ! central cap - core
319          longwatercaptag(26:33)  = .true. ! central cap
[669]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  ----------------------------
[633]327
[1528]328         else if ((nbp_lon.eq.64).and.((nbp_lat-1).eq.48)) then ! 64x48
[633]329
[669]330          print*,'water ice caps distribution for 64x48 resolution'
[633]331          longwatercaptag(1:65)   = .true. ! central cap - core
332          longwatercaptag(75:85)  = .true. ! central cap
333          longwatercaptag(93:114) = .true. ! central cap
[669]334!---------------------   OUTLIERS  ----------------------------
335          if (.true.) then
[633]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
[740]357          longwatercaptag(256)    = .true. ! outlier, lat = 75
[669]358          endif
359!--------------------------------------------------------------       
[633]360
[669]361
[520]362           
[1130]363         else if (klon_glo .ne. 1) then
[669]364       
[1130]365          print*,'No predefined ice location for this resolution :',
[1528]366     &           nbp_lon,nbp_lat-1
[1130]367          print*,'Please change icelocationmode in surfini.F'
368          print*,'Or add some new definitions ...'
[2311]369          call abort_physic("surfini",
370     &         "no pre-definitions for this resolution",1)
[520]371         
[1130]372         endif
[669]373
[1130]374         do ig=1,klon_glo
375          if (longwatercaptag(ig)) watercaptag_glo(ig) = .true.
376         enddo
[283]377
[633]378
[1130]379        ELSE IF (icelocationmode .eq. 3) THEN
[669]380     
[1130]381         print*,'Surfini: ice caps defined by lat and lon values'
[520]382
[1130]383         do ig=1,klon_glo
[669]384         
385c-------- Towards olympia planitia water caps -----------
386c--------------------------------------------------------
[520]387
[1130]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. ) )
[520]392     .         .or.
[283]393
[1130]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. ) )
[520]398     .         .or.
[1130]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.) ) )
[520]403     .         then
404             
405               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
[669]406              !    watercaptag(ig)=.true.
[520]407                  alternate = 1
[669]408               else
[520]409                  alternate = 0
[669]410               endif !end if alternate = 0
411               
[1130]412          endif
[633]413
[669]414c----------- Opposite olympia planitia water cap --------
415c--------------------------------------------------------
[520]416
[1130]417          if ( ( ( lati_glo(ig)*180./pi     .ge.  80 ) .and.
418     .         ( lati_glo(ig)*180./pi     .le.  84 ) )
[669]419     .         .and.
[1130]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
[520]428
429
[669]430c -------------------- Central cap ----------------------
[283]431c--------------------------------------------------------
432
[1130]433          if (abs(lati_glo(ig)*180./pi).gt.80)
434     .          watercaptag_glo(ig)=.true.
[669]435           
436c--------------------------------------------------------
437c--------------------------------------------------------
[1130]438         end do ! of (klon_glo)
[669]439
[2502]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)
[669]447
[2502]448!         print*,'watercaptag_glo(:), ',watercaptag_glo(:)
449
[1130]450        ELSE
[669]451     
452         print*, 'In surfini.F, icelocationmode is ', icelocationmode
[2502]453         print*, 'It should be 1, 2, 3 or 4 (default is 4)'
[2311]454         call abort_physic("surfini","wrong icelocationmode",1)
[669]455
[1130]456        ENDIF ! of if (icelocation)
[669]457       
458       
459        ! print caps locations - useful for plots too
[1130]460        print*,'surfini: latitude | longitude | ig'
461        do ig=1,klon_glo
462          dryness_glo(ig) = icedryness
[669]463
[1130]464          if (watercaptag_glo(ig)) then
465             print*,'surfini: ice water cap', lati_glo(ig)*180./pi,
466     &              long_glo(ig)*180./pi, ig
[2502]467!             write(1,*),ig, lati_glo(ig)*180./pi,
468!     &              cell_area(ig)
469!             write(2,*), lati_glo(ig)*180./pi,
470!     &              long_glo(ig)*180./pi,cell_area(ig)
471!             write(3,*), ig, lati_glo(ig)*180./pi,
472!     &              long_glo(ig)*180./pi,cell_area(ig)
[669]473          endif
474        enddo
475       
[1130]476       endif !of if (is_master)
477       
[2182]478       if (ngrid.gt.1) then
479        ! Now scatter fields watercaptag and dryness from master to all
480        ! (is just a plain copy in serial mode)
481        call scatter(dryness_glo,dryness)
482        call scatter(watercaptag_glo,watercaptag)
483       endif
[1944]484       ELSE
485         watercaptag(:) = .false.
[2741]486       ENDIF ! water
[1944]487! end of #else of #ifndef MESOSCALE
[1212]488#endif       
[2502]489!      END SUBROUTINE surfini(ngrid,piceco2,qsurf)
490      END !SUBROUTINE surfini(ngrid,piceco2,qsurf)
[520]491
[2502]492      SUBROUTINE locate_watercaptag(klon_glo,lati_glo,
493     &            long_glo,watercaptag_glo)
494
495      USE comcstfi_h, ONLY: pi
496
497      integer, intent(in) :: klon_glo
498      real, intent(in) :: lati_glo(klon_glo)
499      real, intent(in) :: long_glo(klon_glo)
500      logical, intent(out) :: watercaptag_glo(klon_glo)
501      integer :: ig,i
502!      real, dimension(klon_glo,120) :: wcap
503      real, dimension(120,2) :: latedge
504      real, dimension(120,2) :: lonedge
505
506! In icelocationmode=2 there are 120 manually predefined grid points where
507! watercaptag is true (for the 64x48 resolution). The grid cells corners
508! coordinates in latitude and longitude are written below. With this
509! routine, we check if the grid cell center is in between any of those
510! points. If so, watercaptag = true.
511
512
513
514
515      latedge(:,1)=(/
516     & 88.125, 84.375, 84.375, 84.375, 84.375, 84.375,84.375, 84.375,
517     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
518     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
519     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
520     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
521     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
522     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
523     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
524     & 84.375, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
525     & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
526     & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
527     & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
528     & 80.625, 80.625, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875,
529     & 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 73.125,
530     & 73.125, 73.125, 73.125, 73.125, 73.125, 73.125, 73.125, 73.125/)
531
532
533      latedge(:,2)=(/
534     & 90. , 88.125, 88.125, 88.125, 88.125, 88.125,88.125, 88.125,   
535     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
536     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
537     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
538     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
539     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
540     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
541     & 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125, 88.125,
542     & 88.125, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
543     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
544     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
545     & 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375, 84.375,
546     & 84.375, 84.375, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625,
547     & 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 80.625, 76.875,
548     & 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875, 76.875/)
549
550
551      lonedge(:,1)=(/
552     &-180.    , -180.    , -177.1875, -171.5625,-165.9375, -160.3125,
553     &-154.6875, -149.0625, -143.4375, -137.8125, -132.1875,-126.5625,
554     &-120.9375, -115.3125, -109.6875, -104.0625,  -98.4375, -92.8125,
555     & -87.1875,  -81.5625,  -75.9375,  -70.3125,  -64.6875, -59.0625,
556     & -53.4375,  -47.8125,  -42.1875,  -36.5625,  -30.9375, -25.3125,
557     & -19.6875,  -14.0625,   -8.4375,   -2.8125,    2.8125,   8.4375,
558     &  14.0625,   19.6875,   25.3125,   30.9375,   36.5625,  42.1875,
559     &  47.8125,   53.4375,   59.0625,   64.6875,   70.3125,  75.9375,
560     &  81.5625,   87.1875,   92.8125,   98.4375,  104.0625, 109.6875,
561     & 115.3125,  120.9375,  126.5625,  132.1875,  137.8125, 143.4375,
562     & 149.0625,  154.6875,  160.3125,  165.9375,  171.5625,-132.1875,
563     &-126.5625, -120.9375, -115.3125, -109.6875, -104.0625, -98.4375,
564     & -92.8125,  -87.1875,  -81.5625,  -75.9375,  -30.9375, -25.3125,
565     & -19.6875,  -14.0625,   -8.4375,   -2.8125,    2.8125,   8.4375,
566     &  14.0625,   19.6875,   25.3125,   30.9375,   36.5625,  42.1875,
567     &  47.8125,   53.4375,   59.0625,   64.6875,   70.3125,  75.9375,
568     &  81.5625,   87.1875, -149.0625, -137.8125, -126.5625,-115.3125,
569     &  -8.4375,    2.8125,   14.0625,  115.3125,  126.5625, 137.8125,
570     & 149.0625,  160.3125,  171.5625, -180.    , -132.1875,-109.6875,
571     &  98.4375,  109.6875,  132.1875,  143.4375,  154.6875,165.9375/)
572
573      lonedge(:,2)=(/
574     & 180.    , -180.    , -171.5625, -165.9375,-160.3125, -154.6875,
575     &-149.0625,-143.4375, -137.8125, -132.1875, -126.5625, -120.9375,
576     &-115.3125,-109.6875, -104.0625,  -98.4375,  -92.8125,  -87.1875,
577     & -81.5625, -75.9375,  -70.3125,  -64.6875,  -59.0625,  -53.4375,
578     & -47.8125, -42.1875,  -36.5625,  -30.9375,  -25.3125,  -19.6875,
579     & -14.0625,  -8.4375,   -2.8125,    2.8125,    8.4375,   14.0625,
580     &  19.6875,  25.3125,   30.9375,   36.5625,   42.1875,   47.8125,
581     &  53.4375,  59.0625,   64.6875,   70.3125,   75.9375,   81.5625,
582     &  87.1875,  92.8125,   98.4375,  104.0625,  109.6875,  115.3125,
583     & 120.9375, 126.5625,  132.1875,  137.8125,  143.4375,  149.0625,
584     & 154.6875, 160.3125,  165.9375,  171.5625,  177.1875, -126.5625,
585     &-120.9375,-115.3125, -109.6875, -104.0625,  -98.4375,  -92.8125,
586     & -87.1875, -81.5625,  -75.9375,  -70.3125,  -25.3125,  -19.6875,
587     & -14.0625,  -8.4375,   -2.8125,    2.8125,    8.4375,   14.0625,
588     &  19.6875,  25.3125,   30.9375,   36.5625,   42.1875,   47.8125,
589     &  53.4375,  59.0625,   64.6875,   70.3125,   75.9375,   81.5625,
590     &  87.1875,  92.8125, -143.4375, -132.1875, -120.9375, -109.6875,
591     &  -2.8125,   8.4375,   19.6875,  120.9375,  132.1875,  143.4375,
592     & 154.6875, 165.9375,  177.1875, -177.1875, -126.5625, -104.0625,
593     & 104.0625, 115.3125,  137.8125,  149.0625,  160.3125,171.5625/)
594
595
596      watercaptag_glo(:) = .false.
597      DO ig=1, klon_glo
598        DO i=1, 120
599           if ((long_glo(ig)*180./pi.ge.lonedge(i,1))
600     &         .and.(long_glo(ig)*180./pi.le.lonedge(i,2))
601     &         .and.(lati_glo(ig)*180./pi.ge.latedge(i,1))
602     &         .and.(lati_glo(ig)*180./pi.le.latedge(i,2))) then
603             watercaptag_glo(ig) = .true.
604           endif
605        ENDDO !i=1, 120
606      ENDDO ! ig=1, klon_glo
607
608      END SUBROUTINE locate_watercaptag
Note: See TracBrowser for help on using the repository browser.