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