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

Last change on this file since 2903 was 2892, checked in by romain.vande, 2 years ago

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