source: trunk/LMDZ.MARS/libf/phymars/surfini_mod.F90 @ 3325

Last change on this file since 3325 was 3142, checked in by jbclement, 12 months ago

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