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

Last change on this file since 2609 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
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
[669]47       
48       
49      !in case icelocationmode == 1
50      INTEGER i,j
51      INTEGER     imd,jmd
52      PARAMETER   (imd=360,jmd=180)
[712]53      REAL        zdata(imd,jmd)
[669]54      REAL        zelat,zelon
55
[1212]56#ifndef MESOSCALE
[1207]57      INTEGER nb_ice(klon_glo,2)   ! number of counts | detected ice for GCM grid
[1212]58#endif
[1528]59      INTEGER latice(nbp_lat-1,2),lonice (nbp_lon,2) ! number of counts | detected ice along lat & lon axis
[669]60
61      REAL step,count,ratiolat
62
63      INTEGER   ierr,nid,nvarid
64     
65      REAL,SAVE :: min_icevalue = 500.
[697]66      character(len=50) :: string = 'thermal'
[669]67     
68      character (len=100) :: zedatafile
[1130]69
[1212]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
[1541]81         if ( ( latitude(ig)*180./pi .gt. 70. ) .and.
[1212]82     .        ( albedodat(ig) .ge. 0.26   ) )  then
83                 write(*,*)"outlier ",ig
84                 watercaptag(ig)  = .true.
85                 dryness(ig)      = 1.
[2508]86                 albedodat(ig)    = albedo_h2o_cap  !! pour output
[1212]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
[1130]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
[1212]109#endif
[1130]110
[1212]111#ifndef MESOSCALE
112
[38]113c
114c=======================================================================
[1130]115! Initialize watercaptag (default is false)
116      watercaptag_glo(:)=.false.
[38]117
[283]118c     water ice outliers
[38]119c     ------------------------------------------
120
[283]121      IF ((water) .and. (caps)) THEN
[285]122     
[283]123c Perennial H20 north cap defined by watercaptag=true (allows surface to be
124c hollowed by sublimation in vdifc).
[38]125
[283]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
[520]138         
[1047]139         if (ngrid .ne. 1) then
[633]140           watercaptag(:) = .false.
141           longwatercaptag(:) = .false.
142         endif
143         
[1130]144         write(*,*) "surfini: Ice dryness ?"
[520]145         icedryness=1. ! default value
[2311]146         call getin_p("icedryness",icedryness)
[1130]147         write(*,*) "surfini: icedryness = ",icedryness
[669]148         dryness (:) = icedryness
[633]149         
[1130]150      ! To be able to run in parallel, we work on the full grid
151      ! and dispatch results afterwards
[669]152
[1130]153      ! start by geting latitudes and logitudes on full grid
154      ! (in serial mode, this is just a copy)
[1541]155      call gather(latitude,lati_glo)
156      call gather(longitude,long_glo)
[669]157
[1130]158      if (is_master) then
159
160        IF (ngrid .eq. 1) THEN ! special case for 1d --> do nothing
[669]161     
[1047]162         print*, 'ngrid = 1, do no put ice caps in surfini.F'
[669]163
[1130]164        ELSE IF (icelocationmode .eq. 1) THEN
[669]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           
[1918]177         zedatafile = trim(datadir)
[669]178 
179       
[1130]180         ierr=nf90_open(trim(zedatafile)//'/surface.nc',
181     &   NF90_NOWRITE,nid)
[669]182     
[1130]183         IF (ierr.NE.nf90_noerr) THEN
[669]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:'
[1381]191       write(*,*)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
[2311]192       call abort_physic("surfini","missing surface.nc file",1)
[1130]193         ENDIF
[669]194     
195     
[1130]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))
[2311]201          call abort_physic("surfini","missing "//trim(string),1)
[1130]202         endif
[697]203
[1130]204         ierr=nf90_get_var(nid, nvarid, zdata)
[697]205
[1130]206         if (ierr.ne.nf90_noerr) then
207          write(*,*) 'surfini: error failed loading ',trim(string)
208          write(*,*)trim(nf90_strerror(ierr))
[2311]209          call abort_physic("surfini","failed loading "//trim(string),1)
[1130]210         endif
[669]211 
212                     
[1130]213         ierr=nf90_close(nid)
[669]214 
[285]215
[1130]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
[669]223
[1130]224         ! loop over the GCM grid - except for poles (ig=1 and ngrid)
225         do ig=2,klon_glo-1
[520]226     
[1130]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
[1528]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
[669]236              ! count all points in that GCM grid point
237              nb_ice(ig,1) = nb_ice(ig,1) + 1
[712]238              if (zdata(i,j) > min_icevalue)
[669]239                 ! count all detected points in that GCM grid point
240     &           nb_ice(ig,2) = nb_ice(ig,2) + 1
[712]241             endif
[1130]242           enddo
243          enddo 
[712]244
[669]245        ! projection of nb_ice on GCM lat and lon axes
[1528]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 ...
[669]250
[1130]251         enddo ! of do ig=2,klon_glo-1
[669]252     
253
254     
[1130]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,:)
[1528]259         latice(nbp_lat-1,:) = nb_ice(ngrid,:)
260         lonice(nbp_lon,:) = nb_ice(ngrid,:)
[669]261     
262     
[712]263!      print*, 'latice TOT', latice(:,1)
264!      print*, 'latice FOUND', latice(:,2)
265!      print*, 'lonice TOT', lonice(:,1)
266!      print*, 'lonice FOUND', lonice(:,2)
[669]267     
[712]268!      print*, 'lat ratio', int(real(latice(:,2))/real(latice(:,1))*iim)
269!      print*, 'lon ratio', int(real(lonice(:,2))/real(lonice(:,1))*jjm)
[669]270     
[712]271!      print*,''
272!      print*,'sum lat', sum(latice(:,1)), sum(lonice(:,1))
273!      print*,'sum lon', sum(latice(:,2)), sum(lonice(:,2))
[669]274     
275   
[1130]276         ! loop over GCM latitudes. CONSIDER ONLY NORTHERN HEMISPHERE
[1528]277         do i=1,(nbp_lat-1)/2
[1130]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
[669]283     
[1130]284          ! put ice caps while there is not enough ice,
285          ! as long as the threshold is above 20%
[1528]286          do while ((count.le.ratiolat*nbp_lon).and.(step.ge.0.2))
[1130]287           count = 0.
288           ! loop over GCM longitudes
[1528]289           do j=1,nbp_lon
[669]290            ! if the detected ice ratio in the GCM grid point
291            ! is more than 'step', then add ice
[1528]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.
[669]295                  count = count + 1
296            endif
[1528]297           enddo ! of do j=1,nbp_lon
298           !print*, 'step',step,count,ratiolat*nbp_lon
[1130]299           step = step - 0.01
300          enddo ! of do while
[1528]301          !print*, 'step',step,count,ratiolat*nbp_lon
[669]302
[1130]303         enddo ! of do i=1,jjm/2
[669]304           
305
[1130]306        ELSE IF (icelocationmode .eq. 2) THEN
[669]307     
[1130]308         print*,'Surfini: predefined ice caps'
[669]309     
[1528]310         if ((nbp_lon.eq.32).and.((nbp_lat-1).eq.24)) then ! 32x24
[633]311           
[669]312          print*,'water ice caps distribution for 32x24 resolution'
[633]313          longwatercaptag(1:9)    = .true. ! central cap - core
314          longwatercaptag(26:33)  = .true. ! central cap
[669]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  ----------------------------
[633]322
[1528]323         else if ((nbp_lon.eq.64).and.((nbp_lat-1).eq.48)) then ! 64x48
[633]324
[669]325          print*,'water ice caps distribution for 64x48 resolution'
[633]326          longwatercaptag(1:65)   = .true. ! central cap - core
327          longwatercaptag(75:85)  = .true. ! central cap
328          longwatercaptag(93:114) = .true. ! central cap
[669]329!---------------------   OUTLIERS  ----------------------------
330          if (.true.) then
[633]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
[740]352          longwatercaptag(256)    = .true. ! outlier, lat = 75
[669]353          endif
354!--------------------------------------------------------------       
[633]355
[669]356
[520]357           
[1130]358         else if (klon_glo .ne. 1) then
[669]359       
[1130]360          print*,'No predefined ice location for this resolution :',
[1528]361     &           nbp_lon,nbp_lat-1
[1130]362          print*,'Please change icelocationmode in surfini.F'
363          print*,'Or add some new definitions ...'
[2311]364          call abort_physic("surfini",
365     &         "no pre-definitions for this resolution",1)
[520]366         
[1130]367         endif
[669]368
[1130]369         do ig=1,klon_glo
370          if (longwatercaptag(ig)) watercaptag_glo(ig) = .true.
371         enddo
[283]372
[633]373
[1130]374        ELSE IF (icelocationmode .eq. 3) THEN
[669]375     
[1130]376         print*,'Surfini: ice caps defined by lat and lon values'
[520]377
[1130]378         do ig=1,klon_glo
[669]379         
380c-------- Towards olympia planitia water caps -----------
381c--------------------------------------------------------
[520]382
[1130]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. ) )
[520]387     .         .or.
[283]388
[1130]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. ) )
[520]393     .         .or.
[1130]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.) ) )
[520]398     .         then
399             
400               if ((alternate .eq. 0)) then  ! 1/2 en 64x48 sinon trop large en lat
[669]401              !    watercaptag(ig)=.true.
[520]402                  alternate = 1
[669]403               else
[520]404                  alternate = 0
[669]405               endif !end if alternate = 0
406               
[1130]407          endif
[633]408
[669]409c----------- Opposite olympia planitia water cap --------
410c--------------------------------------------------------
[520]411
[1130]412          if ( ( ( lati_glo(ig)*180./pi     .ge.  80 ) .and.
413     .         ( lati_glo(ig)*180./pi     .le.  84 ) )
[669]414     .         .and.
[1130]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
[520]423
424
[669]425c -------------------- Central cap ----------------------
[283]426c--------------------------------------------------------
427
[1130]428          if (abs(lati_glo(ig)*180./pi).gt.80)
429     .          watercaptag_glo(ig)=.true.
[669]430           
431c--------------------------------------------------------
432c--------------------------------------------------------
[1130]433         end do ! of (klon_glo)
[669]434
[2502]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)
[669]442
[2502]443!         print*,'watercaptag_glo(:), ',watercaptag_glo(:)
444
[1130]445        ELSE
[669]446     
447         print*, 'In surfini.F, icelocationmode is ', icelocationmode
[2502]448         print*, 'It should be 1, 2, 3 or 4 (default is 4)'
[2311]449         call abort_physic("surfini","wrong icelocationmode",1)
[669]450
[1130]451        ENDIF ! of if (icelocation)
[669]452       
453       
454        ! print caps locations - useful for plots too
[1130]455        print*,'surfini: latitude | longitude | ig'
456        do ig=1,klon_glo
457          dryness_glo(ig) = icedryness
[669]458
[1130]459          if (watercaptag_glo(ig)) then
460             print*,'surfini: ice water cap', lati_glo(ig)*180./pi,
461     &              long_glo(ig)*180./pi, ig
[2502]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)
[669]468          endif
469        enddo
470       
[1130]471       endif !of if (is_master)
472       
[2182]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
[1944]479       ELSE
480         watercaptag(:) = .false.
[520]481       ENDIF ! (caps & water)
[1944]482! end of #else of #ifndef MESOSCALE
[1212]483#endif       
[2502]484!      END SUBROUTINE surfini(ngrid,piceco2,qsurf)
485      END !SUBROUTINE surfini(ngrid,piceco2,qsurf)
[520]486
[2502]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.