New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
coast_dist.F90 in branches/UKMO/dev_r5518_DMP_TOOLS/NEMOGCM/TOOLS/DMP_TOOLS/src – NEMO

source: branches/UKMO/dev_r5518_DMP_TOOLS/NEMOGCM/TOOLS/DMP_TOOLS/src/coast_dist.F90 @ 10199

Last change on this file since 10199 was 10199, checked in by jenniewaters, 6 years ago

Alllow a distance to coast file to be read in. Also modify code to prevent multiple calculations of the surface distance to coast.
File size: 9.8 KB
RevLine 
[4738]1MODULE coastdist
2
3   USE utils
4   USE netcdf
5
6   IMPLICIT NONE
7   PUBLIC
8
9   CONTAINS
10
[10199]11   SUBROUTINE coast_dist_weight( presto, ln_read , klev ) 
[4738]12      !!----------------------------------------------------------------------
13      !!                 *** ROUTINE coast_dist_weight ***
14      !!
15      !! ** Purpose: Weight restoration coefficient by distance to coast
16      !!
17      !! ** Method: 1) Calculate distance to coast
18      !!            2) Reduce resto with 1000km of coast
19      !!
20      IMPLICIT NONE
[4739]21      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) :: presto
[10199]22      LOGICAL, INTENT( in )                         :: ln_read
23      INTEGER, INTENT( in )                         :: klev
[4739]24      REAL(wp), DIMENSION(jpi,jpj) :: zdct
25      REAL(wp) :: zinfl = 1000.e3_wp  ! Distance of influence of coast line (could be
[4738]26                                  ! a namelist setting)
27      INTEGER :: jj, ji           ! dummy loop indices
[10199]28      INTEGER :: ncin, tdist_id
[4738]29
[10199]30      IF ( ln_read ) THEN
31
32         CALL check_nf90( nf90_open('dist_coast_uvtf.nc', NF90_NOWRITE, ncin), 'Error opening dist to coast file' )
33         CALL check_nf90( nf90_inq_varid( ncin, 'tdist', tdist_id ), 'Cannot get variable ID for tmask')
34         CALL check_nf90( nf90_get_var( ncin, tdist_id, zdct, (/ 1,1,klev /), (/ jpi, jpj,1 /) ) )
35         CALL check_nf90( nf90_close(ncin) ) 
36
37      ELSE
38
39         IF (klev == 1 ) THEN ! Since only a the surface distance to coast is used, only calculate once.
40            CALL cofdis( zdct )
41         ENDIF
42
43      ENDIF
44
[4738]45      DO jj = 1, jpj
46         DO ji = 1, jpi
47            zdct(ji,jj) = MIN( zinfl, zdct(ji,jj) )
[4739]48            presto(ji,jj) = presto(ji, jj) * 0.5_wp * (  1._wp - COS( rpi*zdct(ji,jj)/zinfl) )
[4738]49         END DO
50      END DO
51
52   END SUBROUTINE coast_dist_weight
53
54
55   SUBROUTINE cofdis( pdct )
56      !!----------------------------------------------------------------------
57      !!                 ***  ROUTINE cofdis  ***
58      !!
59      !! ** Purpose :   Compute the distance between ocean T-points and the
60      !!      ocean model coastlines.
61      !!
62      !! ** Method  :   For each model level, the distance-to-coast is
63      !!      computed as follows :
64      !!       - The coastline is defined as the serie of U-,V-,F-points
65      !!      that are at the ocean-land bound.
66      !!       - For each ocean T-point, the distance-to-coast is then
67      !!      computed as the smallest distance (on the sphere) between the
68      !!      T-point and all the coastline points.
69      !!       - For land T-points, the distance-to-coast is set to zero.
70      !!
71      !! ** Action  : - pdct, distance to the coastline (argument)
72      !!              - NetCDF file 'dist.coast.nc'
73      !!----------------------------------------------------------------------
74      !!
[4739]75      REAL(wp), DIMENSION(jpi,jpj), INTENT( out ) ::   pdct   ! distance to the coastline
[4738]76      !!
77      INTEGER ::   ji, jj, jl   ! dummy loop indices
78      INTEGER ::   iju, ijt, icoast, itime, ierr, icot   ! local integers
79      CHARACTER (len=32) ::   clname                     ! local name
[4739]80      REAL(wp) ::   zdate0                               ! local scalar
81      REAL(wp), POINTER, DIMENSION(:,:) ::  zxt, zyt, zzt, zmask
82      REAL(wp), POINTER, DIMENSION(:  ) ::  zxc, zyc, zzc, zdis    ! temporary workspace
[4738]83      LOGICAL , ALLOCATABLE, DIMENSION(:,:) ::  llcotu, llcotv, llcotf   ! 2D logical workspace
84
85      !!----------------------------------------------------------------------
86      !
87      ALLOCATE( zxt(jpi,jpj) , zyt(jpi,jpj) , zzt(jpi,jpj) , zmask(jpi,jpj)    )
88      ALLOCATE(zxc(3*jpi*jpj), zyc(3*jpi*jpj), zzc(3*jpi*jpj), zdis(3*jpi*jpj)     )
89      ALLOCATE( llcotu(jpi,jpj), llcotv(jpi,jpj), llcotf(jpi,jpj)  )
90      ALLOCATE( gphiu(jpi,jpj), gphiv(jpi,jpj), gphif(jpi,jpj)  )
91      ALLOCATE( glamu(jpi,jpj), glamv(jpi,jpj), glamf(jpi,jpj), glamt(jpi,jpj)  )
92      ALLOCATE( umask(jpi,jpj), vmask(jpi,jpj), fmask(jpi,jpj)  )
93      !
94
95      CALL check_nf90( nf90_get_var( ncin, gphit_id, gphit, (/ 1,1 /), (/ jpi, jpj /) ) )
96      CALL check_nf90( nf90_get_var( ncin, gphiu_id, gphiu, (/ 1,1 /), (/ jpi, jpj /) ) )
97      CALL check_nf90( nf90_get_var( ncin, gphiv_id, gphiv, (/ 1,1 /), (/ jpi, jpj /) ) )
98      CALL check_nf90( nf90_get_var( ncin, gphif_id, gphif, (/ 1,1 /), (/ jpi, jpj /) ) )
99      CALL check_nf90( nf90_get_var( ncin, glamt_id, glamt, (/ 1,1 /), (/ jpi, jpj /) ) )
100      CALL check_nf90( nf90_get_var( ncin, glamu_id, glamu, (/ 1,1 /), (/ jpi, jpj /) ) )
101      CALL check_nf90( nf90_get_var( ncin, glamv_id, glamv, (/ 1,1 /), (/ jpi, jpj /) ) )
102      CALL check_nf90( nf90_get_var( ncin, glamf_id, glamf, (/ 1,1 /), (/ jpi, jpj /) ) )
103      CALL check_nf90( nf90_get_var( ncin, tmask_id, tmask, (/ 1,1 /), (/ jpi, jpj /) ) )
104      CALL check_nf90( nf90_get_var( ncin, umask_id, umask, (/ 1,1 /), (/ jpi, jpj /) ) )
105      CALL check_nf90( nf90_get_var( ncin, vmask_id, vmask, (/ 1,1 /), (/ jpi, jpj /) ) )
106      CALL check_nf90( nf90_get_var( ncin, fmask_id, fmask, (/ 1,1 /), (/ jpi, jpj /) ) )
107
[4739]108      pdct(:,:) = 0._wp
[4738]109      zxt(:,:) = COS( rad * gphit(:,:) ) * COS( rad * glamt(:,:) )
110      zyt(:,:) = COS( rad * gphit(:,:) ) * SIN( rad * glamt(:,:) )
111      zzt(:,:) = SIN( rad * gphit(:,:) )
112
113
114         ! Define the coastline points (U, V and F)
115         DO jj = 2, jpj-1
116            DO ji = 2, jpi-1
117               zmask(ji,jj) =  ( tmask(ji,jj+1) + tmask(ji+1,jj+1) &
118                   &           + tmask(ji,jj  ) + tmask(ji+1,jj  ) )
[4739]119               llcotu(ji,jj) = ( tmask(ji,jj ) + tmask(ji+1,jj  ) == 1._wp ) 
120               llcotv(ji,jj) = ( tmask(ji,jj  ) + tmask(ji  ,jj+1) == 1._wp ) 
121               llcotf(ji,jj) = ( zmask(ji,jj) > 0._wp ) .AND. ( zmask(ji,jj) < 4._wp )
[4738]122            END DO
123         END DO
124
125         ! Lateral boundaries conditions
126         llcotu(:, 1 ) = umask(:,  2  ) == 1
127         llcotu(:,jpj) = umask(:,jpj-1) == 1
128         llcotv(:, 1 ) = vmask(:,  2  ) == 1
129         llcotv(:,jpj) = vmask(:,jpj-1) == 1
130         llcotf(:, 1 ) = fmask(:,  2  ) == 1
131         llcotf(:,jpj) = fmask(:,jpj-1) == 1
132
133         IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
134            llcotu( 1 ,:) = llcotu(jpi-1,:)
135            llcotu(jpi,:) = llcotu(  2  ,:)
136            llcotv( 1 ,:) = llcotv(jpi-1,:)
137            llcotv(jpi,:) = llcotv(  2  ,:)
138            llcotf( 1 ,:) = llcotf(jpi-1,:)
139            llcotf(jpi,:) = llcotf(  2  ,:)
140         ELSE
141            llcotu( 1 ,:) = umask(  2  ,:) == 1
142            llcotu(jpi,:) = umask(jpi-1,:) == 1
143            llcotv( 1 ,:) = vmask(  2  ,:) == 1
144            llcotv(jpi,:) = vmask(jpi-1,:) == 1
145            llcotf( 1 ,:) = fmask(  2  ,:) == 1
146            llcotf(jpi,:) = fmask(jpi-1,:) == 1
147         ENDIF
148         IF( jperio == 3 .OR. jperio == 4 ) THEN
149            DO ji = 1, jpi-1
150               iju = jpi - ji + 1
151               llcotu(ji,jpj  ) = llcotu(iju,jpj-2)
152               llcotf(ji,jpj-1) = llcotf(iju,jpj-2)
153               llcotf(ji,jpj  ) = llcotf(iju,jpj-3)
154            END DO
155            DO ji = jpi/2, jpi-1
156               iju = jpi - ji + 1
157               llcotu(ji,jpj-1) = llcotu(iju,jpj-1)
158            END DO
159            DO ji = 2, jpi
160               ijt = jpi - ji + 2
161               llcotv(ji,jpj-1) = llcotv(ijt,jpj-2)
162               llcotv(ji,jpj  ) = llcotv(ijt,jpj-3)
163            END DO
164         ENDIF
165         IF( jperio == 5 .OR. jperio == 6 ) THEN
166            DO ji = 1, jpi-1
167               iju = jpi - ji
168               llcotu(ji,jpj  ) = llcotu(iju,jpj-1)
169               llcotf(ji,jpj  ) = llcotf(iju,jpj-2)
170            END DO
171            DO ji = jpi/2, jpi-1
172               iju = jpi - ji
173               llcotf(ji,jpj-1) = llcotf(iju,jpj-1)
174            END DO
175            DO ji = 1, jpi
176               ijt = jpi - ji + 1
177               llcotv(ji,jpj  ) = llcotv(ijt,jpj-1)
178            END DO
179            DO ji = jpi/2+1, jpi
180               ijt = jpi - ji + 1
181               llcotv(ji,jpj-1) = llcotv(ijt,jpj-1)
182            END DO
183         ENDIF
184
185         ! Compute cartesian coordinates of coastline points
186         ! and the number of coastline points
187         icoast = 0
188         DO jj = 1, jpj
189            DO ji = 1, jpi
190               IF( llcotf(ji,jj) ) THEN
191                  icoast = icoast + 1
192                  zxc(icoast) = COS( rad*gphif(ji,jj) ) * COS( rad*glamf(ji,jj) )
193                  zyc(icoast) = COS( rad*gphif(ji,jj) ) * SIN( rad*glamf(ji,jj) )
194                  zzc(icoast) = SIN( rad*gphif(ji,jj) )
195               ENDIF
196               IF( llcotu(ji,jj) ) THEN
197                  icoast = icoast+1
198                  zxc(icoast) = COS( rad*gphiu(ji,jj) ) * COS( rad*glamu(ji,jj) )
199                  zyc(icoast) = COS( rad*gphiu(ji,jj) ) * SIN( rad*glamu(ji,jj) )
200                  zzc(icoast) = SIN( rad*gphiu(ji,jj) )
201               ENDIF
202               IF( llcotv(ji,jj) ) THEN
203                  icoast = icoast+1
204                  zxc(icoast) = COS( rad*gphiv(ji,jj) ) * COS( rad*glamv(ji,jj) )
205                  zyc(icoast) = COS( rad*gphiv(ji,jj) ) * SIN( rad*glamv(ji,jj) )
206                  zzc(icoast) = SIN( rad*gphiv(ji,jj) )
207               ENDIF
208            END DO
209         END DO
210
211         ! Distance for the T-points
212         DO jj = 1, jpj
213            DO ji = 1, jpi
[4739]214               IF( tmask(ji,jj) == 0._wp ) THEN
215                  pdct(ji,jj) = 0._wp
[4738]216               ELSE
217                  DO jl = 1, icoast
218                     zdis(jl) = ( zxt(ji,jj) - zxc(jl) )**2   &
219                        &     + ( zyt(ji,jj) - zyc(jl) )**2   &
220                        &     + ( zzt(ji,jj) - zzc(jl) )**2
221                  END DO
222                  pdct(ji,jj) = ra * SQRT( MINVAL( zdis(1:icoast) ) )
223               ENDIF
224            END DO
225         END DO
226
227      DEALLOCATE( zxt , zyt , zzt , zmask    )
228      DEALLOCATE(zxc, zyc, zzc, zdis    )
229      DEALLOCATE( llcotu, llcotv, llcotf  )
230      DEALLOCATE( gphiu, gphiv, gphif  )
231      DEALLOCATE( glamu, glamv, glamf, glamt  )
232      DEALLOCATE( umask, vmask, fmask  )
233
234   END SUBROUTINE cofdis
235
236END MODULE coastdist
Note: See TracBrowser for help on using the repository browser.