Changeset 5257
- Timestamp:
- 2015-05-08T16:18:50+02:00 (8 years ago)
- Location:
- branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src
- Files:
-
- 2 edited
-
scoord_gen.F90 (modified) (42 diffs)
-
utils.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/scoord_gen.F90
r5255 r5257 44 44 !!---------------------------------------------------------------------- 45 45 ! 46 USE utils 46 47 IMPLICIT NONE 47 48 48 USE utils49 49 50 50 !!---------------------------------------------------------------------- … … 54 54 READ( numnam, namzgr_sco ) 55 55 CLOSE( numnam ) 56 jpk = rn_jpk57 58 WRITE( numout,*)59 WRITE( numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'60 WRITE( numout,*) '~~~~~~~~~~~'61 WRITE( numout,*) ' Namelist namzgr_sco'62 WRITE( numout,*) ' stretching coeffs '63 WRITE( numout,*) ' Number of levels rn_jpk = ',rn_jpk64 WRITE( numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max65 WRITE( numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min66 WRITE( numout,*) ' Critical depth rn_hc = ',rn_hc67 WRITE( numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax68 WRITE( numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh9469 WRITE( numout,*) ' Song and Haidvogel 1994 stretching coefficients'70 WRITE( numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta71 WRITE( numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb72 WRITE( numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb73 WRITE( numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf1274 WRITE( numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit75 WRITE( numout,*) ' Siddorn and Furner 2012 stretching coefficients'76 WRITE( numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha77 WRITE( numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold78 WRITE( numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs79 WRITE( numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a80 WRITE( numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b81 WRITE( numout,*) ' Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'56 jpk = INT(rn_jpk) 57 58 WRITE(*,*) 59 WRITE(*,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate' 60 WRITE(*,*) '~~~~~~~~~~~' 61 WRITE(*,*) ' Namelist namzgr_sco' 62 WRITE(*,*) ' stretching coeffs ' 63 WRITE(*,*) ' Number of levels rn_jpk = ',rn_jpk 64 WRITE(*,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max 65 WRITE(*,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min 66 WRITE(*,*) ' Critical depth rn_hc = ',rn_hc 67 WRITE(*,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax 68 WRITE(*,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94 69 WRITE(*,*) ' Song and Haidvogel 1994 stretching coefficients' 70 WRITE(*,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta 71 WRITE(*,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb 72 WRITE(*,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb 73 WRITE(*,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12 74 WRITE(*,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit 75 WRITE(*,*) ' Siddorn and Furner 2012 stretching coefficients' 76 WRITE(*,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha 77 WRITE(*,*) ' e-fold length scale for transition region rn_efold = ',rn_efold 78 WRITE(*,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs 79 WRITE(*,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a 80 WRITE(*,*) ' Offset for Zb rn_zb_b = ',rn_zb_b 81 WRITE(*,*) ' Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b' 82 82 83 83 … … 102 102 DO jj = 1, jpj 103 103 DO ji = 1, jpi 104 IF( bathy(ji,jj) > 0. _wp) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )104 IF( bathy(ji,jj) > 0. ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 105 105 END DO 106 106 END DO … … 114 114 DO jj = 1, jpj 115 115 DO ji = 1, jpi 116 IF( bathy(ji,jj) == 0. _wp) THEN116 IF( bathy(ji,jj) == 0. ) THEN 117 117 iip1 = MIN( ji+1, jpi ) 118 118 ijp1 = MIN( jj+1, jpj ) … … 120 120 ijm1 = MAX( jj-1, 1 ) 121 121 IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + & 122 & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0. _wp) THEN122 & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0. ) THEN 123 123 zenv(ji,jj) = rn_sbot_min 124 124 ENDIF … … 128 128 ! 129 129 ! smooth the bathymetry (if required) 130 scosrf(:,:) = 0. _wp! ocean surface depth (here zero: no under ice-shelf sea)130 scosrf(:,:) = 0. ! ocean surface depth (here zero: no under ice-shelf sea) 131 131 scobot(:,:) = bathy(:,:) ! ocean bottom depth 132 132 ! 133 133 jl = 0 134 zrmax = 1. _wp134 zrmax = 1. 135 135 ! 136 136 ! 137 137 ! set scaling factor used in reducing vertical gradients 138 zrfact = ( 1. _wp - rn_rmax ) / ( 1._wp+ rn_rmax )138 zrfact = ( 1. - rn_rmax ) / ( 1. + rn_rmax ) 139 139 ! 140 140 ! initialise temporary evelope depth arrays … … 145 145 ! 146 146 ! initialise temporary r-value arrays 147 zri(:,:) = 1. _wp148 zrj(:,:) = 1. _wp147 zri(:,:) = 1. 148 zrj(:,:) = 1. 149 149 ! ! ================ ! 150 DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8 _wp) ! Iterative loop !150 DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8 ) ! Iterative loop ! 151 151 ! ! ================ ! 152 152 jl = jl + 1 153 zrmax = 0. _wp153 zrmax = 0. 154 154 ! we set zrmax from previous r-values (zri and zrj) first 155 155 ! if set after current r-value calculation (as previously) 156 156 ! we could exit DO WHILE prematurely before checking r-value 157 157 ! of current zenv 158 DO jj = 1, nlcj 159 DO ji = 1, nlci 158 ! DO jj = 1, nlcj 159 ! DO ji = 1, nlci 160 DO jj = 1, jpi !jpi or jpim1? 161 DO ji = 1, jpj 160 162 zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) ) 161 163 END DO 162 164 END DO 163 zri(:,:) = 0._wp 164 zrj(:,:) = 0._wp 165 DO jj = 1, nlcj 166 DO ji = 1, nlci 167 iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) 168 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) 169 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN 165 zri(:,:) = 0. 166 zrj(:,:) = 0. 167 ! DO jj = 1, nlci 168 ! DO ji = 1, nlcj 169 ! iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) 170 ! ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) 171 ! IF( (zenv(ji,jj) > 0.) .AND. (zenv(iip1,jj) > 0.)) THEN 172 ! zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 173 ! END IF 174 ! IF( (zenv(ji,jj) > 0.) .AND. (zenv(ji,ijp1) > 0.)) THEN 175 ! zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 176 ! END IF 177 ! IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact 178 ! IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact 179 ! IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact 180 ! IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact 181 ! END DO 182 ! END DO 183 DO jj = 1, jpi-1 184 DO ji = 1, jpj-1 185 iip1 = ji+1 186 ijp1 = jj+1 187 IF( (zenv(ji,jj) > 0.) .AND. (zenv(iip1,jj) > 0.)) THEN 170 188 zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 171 189 END IF 172 IF( (zenv(ji,jj) > 0. _wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN190 IF( (zenv(ji,jj) > 0.) .AND. (zenv(ji,ijp1) > 0.)) THEN 173 191 zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 174 192 END IF … … 182 200 WRITE(*,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax 183 201 ! 184 DO jj = 1, nlcj185 DO ji = 1, nlci202 DO jj = 1, jpi 203 DO ji = 1, jpj 186 204 zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) ) 187 205 END DO 188 206 END DO 189 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero190 207 ! ! ================ ! 191 208 END DO ! End loop ! … … 198 215 ! 199 216 ! Envelope bathymetry saved in hbatt 217 ! TODO - get this section to work 200 218 hbatt(:,:) = zenv(:,:) 201 IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp) THEN202 CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )203 DO jj = 1, jpj204 DO ji = 1, jpi205 ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp)206 hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp- ztaper )207 END DO208 END DO209 ENDIF219 ! IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0. ) THEN 220 ! CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) 221 ! DO jj = 1, jpj 222 ! DO ji = 1, jpi 223 ! ztaper = EXP( -(gphit(ji,jj)/8.)**2. ) 224 ! hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1. - ztaper ) 225 ! END DO 226 ! END DO 227 ! ENDIF 210 228 ! 211 229 ! ! ============================== 212 230 ! ! hbatu, hbatv, hbatf fields 213 231 ! ! ============================== 214 WRITE( numout,*)215 WRITE( numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min232 WRITE(*,*) 233 WRITE(*,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 216 234 hbatu(:,:) = rn_sbot_min 217 235 hbatv(:,:) = rn_sbot_min 218 236 hbatf(:,:) = rn_sbot_min 219 DO jj = 1, jpj m1220 DO ji = 1, jpi m1 ! NO vector opt.221 hbatu(ji,jj) = 0.50 _wp* ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )222 hbatv(ji,jj) = 0.50 _wp* ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )223 hbatf(ji,jj) = 0.25 _wp* ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &237 DO jj = 1, jpj-1 238 DO ji = 1, jpi-1 ! NO vector opt. 239 hbatu(ji,jj) = 0.50 * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) ) 240 hbatv(ji,jj) = 0.50 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) ) 241 hbatf(ji,jj) = 0.25 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) & 224 242 & + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) ) 225 243 END DO … … 229 247 DO jj = 1, jpj 230 248 DO ji = 1, jpi 231 IF( hbatu(ji,jj) == 0. _wp) THEN232 IF( zhbat(ji,jj) == 0. _wp) hbatu(ji,jj) = rn_sbot_min233 IF( zhbat(ji,jj) /= 0. _wp) hbatu(ji,jj) = zhbat(ji,jj)249 IF( hbatu(ji,jj) == 0. ) THEN 250 IF( zhbat(ji,jj) == 0. ) hbatu(ji,jj) = rn_sbot_min 251 IF( zhbat(ji,jj) /= 0. ) hbatu(ji,jj) = zhbat(ji,jj) 234 252 ENDIF 235 253 END DO … … 238 256 DO jj = 1, jpj 239 257 DO ji = 1, jpi 240 IF( hbatv(ji,jj) == 0. _wp) THEN241 IF( zhbat(ji,jj) == 0. _wp) hbatv(ji,jj) = rn_sbot_min242 IF( zhbat(ji,jj) /= 0. _wp) hbatv(ji,jj) = zhbat(ji,jj)258 IF( hbatv(ji,jj) == 0. ) THEN 259 IF( zhbat(ji,jj) == 0. ) hbatv(ji,jj) = rn_sbot_min 260 IF( zhbat(ji,jj) /= 0. ) hbatv(ji,jj) = zhbat(ji,jj) 243 261 ENDIF 244 262 END DO … … 247 265 DO jj = 1, jpj 248 266 DO ji = 1, jpi 249 IF( hbatf(ji,jj) == 0. _wp) THEN250 IF( zhbat(ji,jj) == 0. _wp) hbatf(ji,jj) = rn_sbot_min251 IF( zhbat(ji,jj) /= 0. _wp) hbatf(ji,jj) = zhbat(ji,jj)267 IF( hbatf(ji,jj) == 0. ) THEN 268 IF( zhbat(ji,jj) == 0. ) hbatf(ji,jj) = rn_sbot_min 269 IF( zhbat(ji,jj) /= 0. ) hbatf(ji,jj) = zhbat(ji,jj) 252 270 ENDIF 253 271 END DO … … 260 278 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 261 279 262 WRITE( numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &280 WRITE(*,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), & 263 281 & ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) ) 264 WRITE( numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &282 WRITE(*,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), & 265 283 & ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) ) 266 WRITE( numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &284 WRITE(*,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), & 267 285 & ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) ) 268 WRITE( numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &286 WRITE(*,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), & 269 287 & ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) ) 270 288 !! helsinki … … 326 344 DO ji = 1, jpi 327 345 DO jk = 1, jpk-1 328 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )329 END DO 330 IF( scobot(ji,jj) == 0. _wp) mbathy(ji,jj) = 0331 END DO 332 END DO 333 WRITE( numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )346 IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 347 END DO 348 IF( scobot(ji,jj) == 0. ) mbathy(ji,jj) = 0 349 END DO 350 END DO 351 WRITE(*,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 334 352 335 353 ! min max values over the domain 336 WRITE( numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )337 WRITE( numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &354 WRITE(*,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 355 WRITE(*,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 338 356 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gdep3w_0(:,:,:) ) 339 WRITE( numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), &357 WRITE(*,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), & 340 358 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), & 341 359 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), & 342 360 & ' w ', MINVAL( e3w_0 (:,:,:) ) 343 361 344 WRITE( numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &362 WRITE(*,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 345 363 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gdep3w_0(:,:,:) ) 346 WRITE( numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), &364 WRITE(*,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), & 347 365 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), & 348 366 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), & … … 350 368 351 369 ! selected vertical profiles 352 WRITE( numout,*)353 WRITE( numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)354 WRITE( numout,*) ' ~~~~~~ --------------------'355 WRITE( numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")356 WRITE( numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), &370 WRITE(*,*) 371 WRITE(*,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1) 372 WRITE(*,*) ' ~~~~~~ --------------------' 373 WRITE(*,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')") 374 WRITE(*,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), & 357 375 & e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk ) 358 DO jj = mj0(20), mj1(20)359 DO ji = mi0(20), mi1(20)360 WRITE( numout,*)361 WRITE( numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)362 WRITE( numout,*) ' ~~~~~~ --------------------'363 WRITE( numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")364 WRITE( numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &376 DO jj = 20, 20 377 DO ji = 20, 20 378 WRITE(*,*) 379 WRITE(*,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 380 WRITE(*,*) ' ~~~~~~ --------------------' 381 WRITE(*,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')") 382 WRITE(*,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), & 365 383 & e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 366 384 END DO 367 385 END DO 368 DO jj = mj0(74), mj1(74)369 DO ji = mi0(100), mi1(100)370 WRITE( numout,*)371 WRITE( numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)372 WRITE( numout,*) ' ~~~~~~ --------------------'373 WRITE( numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")374 WRITE( numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &386 DO jj = 74,74 387 DO ji = 100, 100 388 WRITE(*,*) 389 WRITE(*,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 390 WRITE(*,*) ' ~~~~~~ --------------------' 391 WRITE(*,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')") 392 WRITE(*,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), & 375 393 & e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 376 394 END DO … … 383 401 DO jj = 1, jpj 384 402 385 IF( hbatt(ji,jj) > 0. _wp) THEN403 IF( hbatt(ji,jj) > 0.) THEN 386 404 DO jk = 1, mbathy(ji,jj) 387 405 ! check coordinate is monotonically increasing 388 IF ( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp) THEN389 WRITE( numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk390 WRITE( numout,*) 'e3w',fse3w(ji,jj,:)391 WRITE( numout,*) 'e3t',fse3t(ji,jj,:)392 CALLSTOP 1406 IF (e3w_0(ji,jj,jk) <= 0. .OR. e3t_0(ji,jj,jk) <= 0. ) THEN 407 WRITE(*,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 408 WRITE(*,*) 'e3w',e3w_0(ji,jj,:) 409 WRITE(*,*) 'e3t',e3t_0(ji,jj,:) 410 STOP 1 393 411 ENDIF 394 412 ! and check it has never gone negative 395 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp) THEN396 WRITE( numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk397 WRITE( numout,*) 'gdepw',fsdepw(ji,jj,:)398 WRITE( numout,*) 'gdept',fsdept(ji,jj,:)399 CALLSTOP 1413 IF( gdepw_0(ji,jj,jk) < 0. .OR. gdept_0(ji,jj,jk) < 0. ) THEN 414 WRITE(*,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 415 WRITE(*,*) 'gdepw',gdepw_0(ji,jj,:) 416 WRITE(*,*) 'gdept',gdept_0(ji,jj,:) 417 STOP 1 400 418 ENDIF 401 419 ! and check it never exceeds the total depth 402 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN403 WRITE( numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk404 WRITE( numout,*) 'gdepw',fsdepw(ji,jj,:)405 CALLSTOP 1420 IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN 421 WRITE(*,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 422 WRITE(*,*) 'gdepw',gdepw_0(ji,jj,:) 423 STOP 1 406 424 ENDIF 407 425 END DO … … 409 427 DO jk = 1, mbathy(ji,jj)-1 410 428 ! and check it never exceeds the total depth 411 IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN412 WRITE( numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk413 WRITE( numout,*) 'gdept',fsdept(ji,jj,:)414 CALLSTOP 1429 IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN 430 WRITE(*,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 431 WRITE(*,*) 'gdept',gdept_0(ji,jj,:) 432 STOP 1 415 433 ENDIF 416 434 END DO … … 421 439 END DO 422 440 ! 423 ! 441 ! Write output file 442 CALL write_coord_file() 443 ! 444 ! 445 END PROGRAM SCOORD_GEN 424 446 425 447 !!====================================================================== … … 437 459 !!---------------------------------------------------------------------- 438 460 ! 439 INTEGER :: ji, jj, jk ! dummy loop argument461 USE utils 440 462 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 441 463 ! … … 449 471 ALLOCATE( z_esigwv3(jpi,jpj,jpk) ) 450 472 451 z_gsigw3 = 0. _wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp452 z_esigt3 = 0. _wp ; z_esigw3 = 0._wp453 z_esigtu3 = 0. _wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp454 z_esigwu3 = 0. _wp ; z_esigwv3 = 0._wp473 z_gsigw3 = 0. ; z_gsigt3 = 0. ; z_gsi3w3 = 0. 474 z_esigt3 = 0. ; z_esigw3 = 0. 475 z_esigtu3 = 0. ; z_esigtv3 = 0. ; z_esigtf3 = 0. 476 z_esigwu3 = 0. ; z_esigwv3 = 0. 455 477 456 478 DO ji = 1, jpi … … 459 481 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 460 482 DO jk = 1, jpk 461 z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5 _wp, rn_bb )483 z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5, rn_bb ) 462 484 z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) 463 485 END DO … … 465 487 DO jk = 1, jpk 466 488 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp) 467 z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5 _wp) / REAL(jpk-1,wp)489 z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5 ) / REAL(jpk-1,wp) 468 490 END DO 469 491 ENDIF … … 473 495 z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 474 496 END DO 475 z_esigw3(ji,jj,1 ) = 2. _wp* ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) )476 z_esigt3(ji,jj,jpk) = 2. _wp* ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )497 z_esigw3(ji,jj,1 ) = 2. * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) ) 498 z_esigt3(ji,jj,jpk) = 2. * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) ) 477 499 ! 478 500 ! Coefficients for vertical depth as the sum of e3w scale factors 479 z_gsi3w3(ji,jj,1) = 0.5 _wp* z_esigw3(ji,jj,1)501 z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1) 480 502 DO jk = 2, jpk 481 503 z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) … … 483 505 ! 484 506 DO jk = 1, jpk 485 zcoeft = ( REAL(jk,wp) - 0.5 _wp) / REAL(jpk-1,wp)486 zcoefw = ( REAL(jk,wp) - 1.0 _wp) / REAL(jpk-1,wp)487 gdept_0 488 gdepw_0 507 zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp) 508 zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp) 509 gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 510 gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 489 511 gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 490 512 END DO … … 543 565 !!---------------------------------------------------------------------- 544 566 ! 545 INTEGER :: ji, jj, jk ! dummy loop argument567 USE utils 546 568 REAL(wp) :: zsmth ! smoothing around critical depth 547 569 REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space … … 556 578 ALLOCATE( z_esigwv3(jpi,jpj,jpk) ) 557 579 558 z_gsigw3 = 0. _wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp559 z_esigt3 = 0. _wp ; z_esigw3 = 0._wp560 z_esigtu3 = 0. _wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp561 z_esigwu3 = 0. _wp ; z_esigwv3 = 0._wp580 z_gsigw3 = 0. ; z_gsigt3 = 0. ; z_gsi3w3 = 0. 581 z_esigt3 = 0. ; z_esigw3 = 0. 582 z_esigtu3 = 0. ; z_esigtv3 = 0. ; z_esigtf3 = 0. 583 z_esigwu3 = 0. ; z_esigwv3 = 0. 562 584 563 585 DO ji = 1, jpi … … 568 590 zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,. 569 591 ! could be changed by users but care must be taken to do so carefully 570 zzb = 1.0 _wp-(zzb/hbatt(ji,jj))592 zzb = 1.0-(zzb/hbatt(ji,jj)) 571 593 572 594 zzs = rn_zs / hbatt(ji,jj) 573 595 574 IF (rn_efold /= 0.0 _wp) THEN596 IF (rn_efold /= 0.0) THEN 575 597 zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) 576 598 ELSE 577 zsmth = 1.0 _wp599 zsmth = 1.0 578 600 ENDIF 579 601 580 602 DO jk = 1, jpk 581 603 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) 582 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5 _wp)/REAL(jpk-1,wp)604 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 583 605 ENDDO 584 606 z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth ) 585 607 z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth ) 586 608 587 ELSE IF 609 ELSE IF(ln_sigcrit) THEN ! shallow water, uniform sigma 588 610 589 611 DO jk = 1, jpk … … 596 618 DO jk = 1, jpk 597 619 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 598 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5 _wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))620 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 599 621 END DO 600 622 … … 605 627 z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 606 628 END DO 607 z_esigw3(ji,jj,1 ) = 2.0 _wp* (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ))608 z_esigt3(ji,jj,jpk) = 2.0 _wp* (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))629 z_esigw3(ji,jj,1 ) = 2.0 * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 )) 630 z_esigt3(ji,jj,jpk) = 2.0 * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk)) 609 631 610 632 ! Coefficients for vertical depth as the sum of e3w scale factors … … 615 637 616 638 DO jk = 1, jpk 617 gdept_0 618 gdepw_0 639 gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 640 gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 619 641 gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 620 642 END DO … … 671 693 !!---------------------------------------------------------------------- 672 694 673 INTEGER :: ji, jj, jk ! dummy loop argument695 USE utils 674 696 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 675 697 … … 680 702 ALLOCATE( z_esigt(jpk), z_esigw(jpk) ) 681 703 682 z_gsigw = 0. _wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp683 z_esigt = 0. _wp ; z_esigw = 0._wp704 z_gsigw = 0. ; z_gsigt = 0. ; z_gsi3w = 0. 705 z_esigt = 0. ; z_esigw = 0. 684 706 685 707 DO jk = 1, jpk 686 z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5 _wp)708 z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5 ) 687 709 z_gsigt(jk) = -fssig( REAL(jk,wp) ) 688 710 END DO … … 696 718 z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk) 697 719 END DO 698 z_esigw( 1 ) = 2. _wp* ( z_gsigt(1 ) - z_gsigw(1 ) )699 z_esigt(jpk) = 2. _wp* ( z_gsigt(jpk) - z_gsigw(jpk) )720 z_esigw( 1 ) = 2. * ( z_gsigt(1 ) - z_gsigw(1 ) ) 721 z_esigt(jpk) = 2. * ( z_gsigt(jpk) - z_gsigw(jpk) ) 700 722 ! 701 723 ! Coefficients for vertical depth as the sum of e3w scale factors 702 z_gsi3w(1) = 0.5 _wp* z_esigw(1)724 z_gsi3w(1) = 0.5 * z_esigw(1) 703 725 DO jk = 2, jpk 704 726 z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk) … … 706 728 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 707 729 DO jk = 1, jpk 708 zcoeft = ( REAL(jk,wp) - 0.5 _wp) / REAL(jpk-1,wp)709 zcoefw = ( REAL(jk,wp) - 1.0 _wp) / REAL(jpk-1,wp)710 gdept_0 711 gdepw_0 730 zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp) 731 zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp) 732 gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 733 gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 712 734 gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 713 735 END DO … … 721 743 e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpk-1,wp) ) 722 744 ! 723 e3w_0 745 e3w_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) 724 746 e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpk-1,wp) ) 725 747 e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpk-1,wp) ) … … 744 766 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 745 767 !!---------------------------------------------------------------------- 768 USE utils, ONLY : wp 746 769 REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate 747 770 REAL(wp) :: pf ! sigma value 748 771 !!---------------------------------------------------------------------- 749 772 ! 750 pf = ( TANH( rn_theta * ( -(pk-0.5 _wp) / REAL(jpk-1) + rn_thetb ) ) &773 pf = ( TANH( rn_theta * ( -(pk-0.5) / REAL(jpk-1) + rn_thetb ) ) & 751 774 & - TANH( rn_thetb * rn_theta ) ) & 752 775 & * ( COSH( rn_theta ) & 753 & + COSH( rn_theta * ( 2. _wp * rn_thetb - 1._wp) ) ) &754 & / ( 2. _wp* SINH( rn_theta ) )776 & + COSH( rn_theta * ( 2. * rn_thetb - 1. ) ) ) & 777 & / ( 2. * SINH( rn_theta ) ) 755 778 ! 756 779 END FUNCTION fssig … … 768 791 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 769 792 !!---------------------------------------------------------------------- 793 USE utils, ONLY : wp 770 794 REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate 771 795 REAL(wp), INTENT(in) :: pbb ! Stretching coefficient … … 774 798 ! 775 799 IF ( rn_theta == 0 ) then ! uniform sigma 776 pf1 = - ( pk1 - 0.5 _wp) / REAL( jpk-1 )800 pf1 = - ( pk1 - 0.5 ) / REAL( jpk-1 ) 777 801 ELSE ! stretched sigma 778 pf1 = ( 1. _wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpk-1)) ) ) / SINH( rn_theta ) &779 & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5 _wp)/REAL(jpk-1)) + 0.5_wp) ) - TANH( 0.5_wp* rn_theta ) ) &780 & / ( 2. _wp * TANH( 0.5_wp* rn_theta ) ) )802 pf1 = ( 1. - pbb ) * ( SINH( rn_theta*(-(pk1-0.5)/REAL(jpk-1)) ) ) / SINH( rn_theta ) & 803 & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5)/REAL(jpk-1)) + 0.5) ) - TANH( 0.5 * rn_theta ) ) & 804 & / ( 2. * TANH( 0.5 * rn_theta ) ) ) 781 805 ENDIF 782 806 ! … … 801 825 !! Reference : Siddorn and Furner, in prep 802 826 !!---------------------------------------------------------------------- 827 USE utils, ONLY : jpk,jk,wp 803 828 REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate 804 829 REAL(wp) :: p_gamma(jpk) ! stretched coordinate … … 809 834 REAL(wp) :: zn1,zn2 ! local variables 810 835 REAL(wp) :: za,zb,zx ! local variables 811 integer :: jk812 836 !!---------------------------------------------------------------------- 813 837 ! … … 816 840 zn2 = 1. - zn1 817 841 818 za1 = (rn_alpha+2.0 _wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp)819 za2 = (rn_alpha+2.0 _wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)820 za3 = (zn2**3.0 _wp - za2)/( zn1**3.0_wp- za1)842 za1 = (rn_alpha+2.0)*zn1**(rn_alpha+1.0)-(rn_alpha+1.0)*zn1**(rn_alpha+2.0) 843 za2 = (rn_alpha+2.0)*zn2**(rn_alpha+1.0)-(rn_alpha+1.0)*zn2**(rn_alpha+2.0) 844 za3 = (zn2**3.0 - za2)/( zn1**3.0 - za1) 821 845 822 846 za = pzb - za3*(pzs-za1)-za2 823 za = za/( zn2-0.5 _wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )824 zb = (pzs - za1 - za*( zn1-0.5 _wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp- za1)825 zx = 1.0 _wp-za/2.0_wp-zb847 za = za/( zn2-0.5*(za2+zn2**2.0) - za3*(zn1-0.5*(za1+zn1**2.0) ) ) 848 zb = (pzs - za1 - za*( zn1-0.5*(za1+zn1**2.0 ) ) ) / (zn1**3.0 - za1) 849 zx = 1.0-za/2.0-zb 826 850 827 851 DO jk = 1, jpk 828 p_gamma(jk) = za*(pk1(jk)*(1.0 _wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp+ &829 & zx*( (rn_alpha+2.0 _wp)*pk1(jk)**(rn_alpha+1.0_wp)- &830 & (rn_alpha+1.0 _wp)*pk1(jk)**(rn_alpha+2.0_wp) )831 p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0 _wp-psmth)852 p_gamma(jk) = za*(pk1(jk)*(1.0-pk1(jk)/2.0))+zb*pk1(jk)**3.0 + & 853 & zx*( (rn_alpha+2.0)*pk1(jk)**(rn_alpha+1.0)- & 854 & (rn_alpha+1.0)*pk1(jk)**(rn_alpha+2.0) ) 855 p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0-psmth) 832 856 ENDDO 833 857 … … 835 859 END FUNCTION fgamma 836 860 837 END PROGRAM SCOORD_GEN -
branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/utils.F90
r5255 r5257 1 1 MODULE utils 2 2 3 IMPLICIT NONE4 3 USE netcdf 5 4 5 IMPLICIT NONE 6 6 PUBLIC ! allows the acces to par_oce when dom_oce is used 7 7 ! ! exception to coding rules... to be suppressed ??? 8 8 9 PUBLIC dom_oce_alloc 10 PUBLIC read_bathy 9 ! PUBLIC dom_oce_alloc 11 10 11 INTEGER, PARAMETER :: dp=8 , sp=4, wp=dp 12 12 13 13 !! All coordinates … … 45 45 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers 46 46 INTEGER :: ios ! Local integer output status for namelist read and allocation 47 INTEGER :: numnam ! File handle for namelist 47 48 REAL(wp) :: zrmax, ztaper ! temporary scalars 48 49 REAL(wp) :: zrfact 49 50 ! 50 REAL(wp), DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 51 REAL(wp), DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 51 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 52 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zenv, ztmp, zmsk, zri, zrj, zhbat 53 54 !Namelist variables 55 REAL(wp) :: rn_jpk, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax, rn_theta 56 REAL(wp) :: rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 57 LOGICAL :: ln_s_sh94, ln_s_sf12, ln_sigcrit 52 58 53 59 NAMELIST/namzgr_sco/rn_jpk, ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 54 60 rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 55 61 62 CONTAINS 56 63 57 64 INTEGER FUNCTION dom_oce_alloc() 58 65 !!---------------------------------------------------------------------- 59 INTEGER, DIMENSION( 12) :: ierr66 INTEGER, DIMENSION(4) :: ierr 60 67 !!---------------------------------------------------------------------- 61 68 ierr(:) = 0 62 69 ! 63 70 ALLOCATE( zenv(jpi,jpj), ztmp(jpi,jpj), zmsk(jpi,jpj), zri(jpi,jpj), zrj(jpi,jpj), & 64 & zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj) )71 & zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj), STAT=ierr(1) ) 65 72 ! 66 ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0 67 & gdept_0 68 & gdepw_0 (jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , STAT=ierr(4) )73 ALLOCATE( gdep3w_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) , & 74 & gdept_0(jpi,jpj,jpk) , e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) , & 75 & gdepw_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , STAT=ierr(2) ) 69 76 ! 70 77 ! … … 74 81 & scosrf(jpi,jpj) , scobot(jpi,jpj) , & 75 82 & hifv (jpi,jpj) , hiff (jpi,jpj) , & 76 & hift (jpi,jpj) , hifu (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr( 8) )83 & hift (jpi,jpj) , hifu (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(3) ) 77 84 78 ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(9) )79 85 ALLOCATE( mbathy(jpi,jpj) , STAT=ierr(4) ) 86 ! 80 87 dom_oce_alloc = MAXVAL(ierr) 81 88 ! 82 89 END FUNCTION dom_oce_alloc 83 90 84 91 85 92 SUBROUTINE read_bathy() 86 93 !! Read bathymetry from input netcdf file 87 INTEGER :: var_id 94 INTEGER :: var_id, ncin 88 95 89 CALL check_nf90( nf90_open('bathy.nc', NF90_NOWRITE, ncin), 'Error opening mesh_maskfile' )96 CALL check_nf90( nf90_open('bathy.nc', NF90_NOWRITE, ncin), 'Error opening bathy.nc file' ) 90 97 91 98 ! Find the size of the input bathymetry … … 96 103 97 104 ! Read the bathymetry variable from file 98 CALL check_nf90( nf90_inq_varid( ncin, 'bathymetry', tmask_id ), 'Cannot get variable ID for bathymetry')105 CALL check_nf90( nf90_inq_varid( ncin, 'bathymetry', var_id ), 'Cannot get variable ID for bathymetry') 99 106 CALL check_nf90( nf90_get_var( ncin, var_id, bathy, (/ 1,1 /), (/ jpi, jpj /) ) ) 107 108 CALL check_nf90( nf90_close(ncin), 'Error closing bathy.nc file' ) 100 109 101 110 END SUBROUTINE read_bathy … … 113 122 CALL check_nf90( nf90_inquire_dimension(ncid,id_var,len=len)) 114 123 115 END SUBROUTINE dimlen 124 END SUBROUTINE dimlen 125 126 127 SUBROUTINE write_coord_file() 128 ! Write out variables to the a netcdf coordinates file 129 130 INTEGER :: id_x, id_y, id_z 131 INTEGER :: ncout 132 INTEGER, DIMENSION(11) :: var_ids !Array to contain all variable IDs 116 133 134 !Create the file 135 CALL check_nf90( nf90_create('coord_zgr.nc', NF90_CLOBBER, ncout), 'Could not create output file') 136 ! 137 !Define dimensions 138 CALL check_nf90( nf90_def_dim(ncout, 'x', jpi, id_x) ) 139 CALL check_nf90( nf90_def_dim(ncout, 'y', jpj, id_y) ) 140 CALL check_nf90( nf90_def_dim(ncout, 'z', jpk, id_z) ) 141 ! 142 !Define variables 143 CALL check_nf90( nf90_def_var(ncout, 'gdept_0', nf90_double, (/id_x, id_y,id_x/), var_ids(1)) ) 144 CALL check_nf90( nf90_def_var(ncout, 'gdepw_0', nf90_double, (/id_x, id_y,id_x/), var_ids(2)) ) 145 CALL check_nf90( nf90_def_var(ncout, 'gdep3w_0', nf90_double, (/id_x, id_y,id_x/), var_ids(3)) ) 146 CALL check_nf90( nf90_def_var(ncout, 'e3f_0', nf90_double, (/id_x, id_y,id_x/), var_ids(4)) ) 147 CALL check_nf90( nf90_def_var(ncout, 'e3t_0', nf90_double, (/id_x, id_y,id_x/), var_ids(5)) ) 148 CALL check_nf90( nf90_def_var(ncout, 'e3u_0', nf90_double, (/id_x, id_y,id_x/), var_ids(6)) ) 149 CALL check_nf90( nf90_def_var(ncout, 'e3v_0', nf90_double, (/id_x, id_y,id_x/), var_ids(7)) ) 150 CALL check_nf90( nf90_def_var(ncout, 'e3w_0', nf90_double, (/id_x, id_y,id_x/), var_ids(8)) ) 151 CALL check_nf90( nf90_def_var(ncout, 'e3uw_0', nf90_double, (/id_x, id_y,id_x/), var_ids(9)) ) 152 CALL check_nf90( nf90_def_var(ncout, 'e3vw_0', nf90_double, (/id_x, id_y,id_x/), var_ids(10)) ) 153 CALL check_nf90( nf90_def_var(ncout, 'mbathy', nf90_double, (/id_x, id_y,id_x/), var_ids(11)) ) 154 155 ! End define mode 156 CALL check_nf90( nf90_enddef(ncout) ) 117 157 118 119 SUBROUTINE write_coord_file() 158 !Write variables to file 159 CALL check_nf90( nf90_put_var(ncout, var_ids(1), gdept_0) ) 160 CALL check_nf90( nf90_put_var(ncout, var_ids(2), gdepw_0) ) 161 CALL check_nf90( nf90_put_var(ncout, var_ids(3), gdep3w_0) ) 162 CALL check_nf90( nf90_put_var(ncout, var_ids(4), e3f_0) ) 163 CALL check_nf90( nf90_put_var(ncout, var_ids(5), e3t_0) ) 164 CALL check_nf90( nf90_put_var(ncout, var_ids(6), e3u_0) ) 165 CALL check_nf90( nf90_put_var(ncout, var_ids(7), e3v_0) ) 166 CALL check_nf90( nf90_put_var(ncout, var_ids(8), e3w_0) ) 167 CALL check_nf90( nf90_put_var(ncout, var_ids(9), e3uw_0) ) 168 CALL check_nf90( nf90_put_var(ncout, var_ids(10), e3vw_0) ) 169 CALL check_nf90( nf90_put_var(ncout, var_ids(11), mbathy) ) 170 171 CALL check_nf90( nf90_close(ncout) ) 120 172 121 173 END SUBROUTINE write_coord_file … … 127 179 128 180 IF (istat /= nf90_noerr) THEN 129 WRITE( numerr,*) 'ERROR! : '//TRIM(nf90_strerror(istat))130 IF ( PRESENT(message) ) THEN ; WRITE( numerr,*) message ; ENDIF181 WRITE(*,*) 'ERROR! : '//TRIM(nf90_strerror(istat)) 182 IF ( PRESENT(message) ) THEN ; WRITE(*,*) message ; ENDIF 131 183 STOP 132 184 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
