Changeset 5269 for branches/2015/dev_r5187_UKMO13_simplification
- Timestamp:
- 2015-05-13T15:28:30+02:00 (9 years ago)
- Location:
- branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src
- Files:
-
- 2 edited
-
scoord_gen.F90 (modified) (9 diffs)
-
utils.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/scoord_gen.F90
r5257 r5269 286 286 WRITE(*,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), & 287 287 & ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) ) 288 !! helsinki 288 289 ! Create the output file 290 CALL make_coord_file() 289 291 290 292 ! ! ======================= … … 293 295 ! 294 296 ! non-dimensional "sigma" for model level depth at w- and t-levels 297 295 298 296 299 !======================================================================== … … 298 301 ! Siddorn and Furner 2012 (ln_sf12=T) 299 302 ! or tanh function (both false) 303 ! To reduce memory loop over jpk and write each level to file 300 304 !======================================================================== 301 305 IF ( ln_s_sh94 ) THEN … … 305 309 ELSE 306 310 CALL s_tanh() 307 ENDIF308 309 !310 where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0311 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0312 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0313 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0314 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0315 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0316 where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0317 318 #if defined key_agrif319 ! Ensure meaningful vertical scale factors in ghost lines/columns320 ! TODO: How can we deal with this in offline tool?321 IF( .NOT. Agrif_Root() ) THEN322 !323 IF((nbondi == -1).OR.(nbondi == 2)) THEN324 e3u_0(1,:,:) = e3u_0(2,:,:)325 ENDIF326 !327 IF((nbondi == 1).OR.(nbondi == 2)) THEN328 e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)329 ENDIF330 !331 IF((nbondj == -1).OR.(nbondj == 2)) THEN332 e3v_0(:,1,:) = e3v_0(:,2,:)333 ENDIF334 !335 IF((nbondj == 1).OR.(nbondj == 2)) THEN336 e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)337 ENDIF338 !339 311 ENDIF 340 #endif 341 !! 342 ! HYBRID : 343 DO jj = 1, jpj 344 DO ji = 1, jpi 345 DO jk = 1, jpk-1 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(:,:) ) 352 353 ! min max values over the domain 354 WRITE(*,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 355 WRITE(*,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 356 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gdep3w_0(:,:,:) ) 357 WRITE(*,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), & 358 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), & 359 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), & 360 & ' w ', MINVAL( e3w_0 (:,:,:) ) 361 362 WRITE(*,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 363 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gdep3w_0(:,:,:) ) 364 WRITE(*,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), & 365 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), & 366 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), & 367 & ' w ', MAXVAL( e3w_0 (:,:,:) ) 368 369 ! selected vertical profiles 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), & 375 & e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk ) 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), & 383 & e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 384 END DO 385 END DO 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), & 393 & e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 394 END DO 395 END DO 396 397 !================================================================================ 398 ! check the coordinate makes sense 399 !================================================================================ 400 DO ji = 1, jpi 401 DO jj = 1, jpj 402 403 IF( hbatt(ji,jj) > 0.) THEN 404 DO jk = 1, mbathy(ji,jj) 405 ! check coordinate is monotonically increasing 406 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 411 ENDIF 412 ! and check it has never gone negative 413 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 418 ENDIF 419 ! and check it never exceeds the total depth 420 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 424 ENDIF 425 END DO 426 427 DO jk = 1, mbathy(ji,jj)-1 428 ! and check it never exceeds the total depth 429 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 433 ENDIF 434 END DO 435 436 ENDIF 437 438 END DO 439 END DO 440 ! 441 ! Write output file 442 CALL write_coord_file() 443 ! 312 313 CALL check_nf90( nf90_put_var(ncout, var_ids(11), mbathy) ) 314 CALL check_nf90( nf90_close(ncout) ) 315 316 444 317 ! 445 318 END PROGRAM SCOORD_GEN … … 462 335 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 463 336 ! 464 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 465 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3 466 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 467 468 ALLOCATE( z_gsigw3(jpi,jpj,jpk), z_gsigt3(jpi,jpj,jpk), z_gsi3w3(jpi,jpj,jpk) ) 469 ALLOCATE( z_esigt3(jpi,jpj,jpk), z_esigw3(jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk)) 470 ALLOCATE( z_esigtv3(jpi,jpj,jpk), z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk)) 471 ALLOCATE( z_esigwv3(jpi,jpj,jpk) ) 337 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 338 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1 339 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigt3, z_esigw3, z_esigtu3 340 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 341 342 ALLOCATE( z_gsigw3(jpi,jpj), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) ) 343 ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj)) 344 ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj)) 345 ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) ) 346 ALLOCATE( z_gsi3w3m1(jpi,jpj) ) 347 348 z_gsigw3 = 0. ; z_gsigt3 = 0. ; z_gsi3w3 = 0. 349 z_esigt3 = 0. ; z_esigw3 = 0. 350 z_esigt3p1 = 0. ; z_esigw3p1 = 0. 351 z_esigtu3 = 0. ; z_esigtv3 = 0. ; z_esigtf3 = 0. 352 z_esigwu3 = 0. ; z_esigwv3 = 0. 353 354 DO jk = 1,jpk 355 DO ji = 1, jpi 356 DO jj = 1, jpj 357 358 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 359 z_gsigw3(ji,jj) = -fssig1( REAL(jk,wp)-0.5, rn_bb ) 360 z_gsigw3p1(ji,jj) = -fssig1( REAL(jk+1,wp)-0.5, rn_bb ) 361 z_gsigt3(ji,jj) = -fssig1( REAL(jk,wp) , rn_bb ) 362 ELSE ! shallow water, uniform sigma 363 z_gsigw3(ji,jj) = REAL(jk-1,wp) / REAL(jpk-1,wp) 364 z_gsigw3p1(ji,jj) = REAL(jk,wp) / REAL(jpk-1,wp) 365 z_gsigt3(ji,jj) = ( REAL(jk-1,wp) + 0.5 ) / REAL(jpk-1,wp) 366 ENDIF 367 ! 368 !gsi3w3m1 & gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk 369 IF( jk .EQ. 1) THEN 370 z_esigw3(ji,jj ) = 2. * ( z_gsigt3(ji,jj ) - z_gsigw3(ji,jj ) ) 371 z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj) 372 ELSE 373 z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj) 374 z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj) 375 ENDIF 376 IF( jk .EQ. jpk) THEN 377 z_esigt3(ji,jj) = 2. * ( z_gsigt3(ji,jj) - z_gsigw3(ji,jj) ) 378 ELSE 379 z_esigt3(ji,jj ) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj) 380 ENDIF 381 ! 382 383 zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp) 384 zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp) 385 gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj)+rn_hc*zcoeft ) 386 gdepw_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj)+rn_hc*zcoefw ) 387 gdep3w_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj)+rn_hc*zcoeft ) 388 ! 389 END DO ! for all jj's 390 END DO ! for all ji's 391 392 DO ji = 1, jpim1 393 DO jj = 1, jpjm1 394 z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) ) & 395 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 396 z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) ) & 397 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 398 z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) & 399 & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) ) & 400 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 401 z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) ) & 402 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 403 z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) ) & 404 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 405 ! 406 e3t_0(ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj) + rn_hc/REAL(jpk-1,wp) ) 407 e3u_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) 408 e3v_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) 409 e3f_0(ji,jj) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) 410 ! 411 e3w_0 (ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj) + rn_hc/REAL(jpk-1,wp) ) 412 e3uw_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) 413 e3vw_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj) + rn_hc/REAL(jpk-1,wp) ) 414 END DO 415 END DO 416 417 z_gsigt3m1 = z_gsigt3 418 z_gsiw3m1 = z_gsiw3 419 420 where (e3t_0 (:,:).eq.0.0) e3t_0(:,:) = 1.0 421 where (e3u_0 (:,:).eq.0.0) e3u_0(:,:) = 1.0 422 where (e3v_0 (:,:).eq.0.0) e3v_0(:,:) = 1.0 423 where (e3f_0 (:,:).eq.0.0) e3f_0(:,:) = 1.0 424 where (e3w_0 (:,:).eq.0.0) e3w_0(:,:) = 1.0 425 where (e3uw_0 (:,:).eq.0.0) e3uw_0(:,:) = 1.0 426 where (e3vw_0 (:,:).eq.0.0) e3vw_0(:,:) = 1.0 427 428 CALL write_netcdf_vars(jk) 429 DO jj = 1, jpj 430 DO ji = 1, jpi 431 IF( scobot(ji,jj) >= gdept_0(ji,jj) ) mbathy(ji,jj) = MAX( 2, jk ) 432 IF( scobot(ji,jj) == 0. ) mbathy(ji,jj) = 0 433 END DO 434 END DO 435 END DO !End of loop over jk 436 437 DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1) 438 DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 439 440 END SUBROUTINE s_sh94 441 442 SUBROUTINE s_sf12 443 444 !!---------------------------------------------------------------------- 445 !! *** ROUTINE s_sf12 *** 446 !! 447 !! ** Purpose : stretch the s-coordinate system 448 !! 449 !! ** Method : s-coordinate stretch using the Siddorn and Furner 2012? 450 !! mixed S/sigma/Z coordinate 451 !! 452 !! This method allows the maintenance of fixed surface and or 453 !! bottom cell resolutions (cf. geopotential coordinates) 454 !! within an analytically derived stretched S-coordinate framework. 455 !! 456 !! 457 !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 458 !!---------------------------------------------------------------------- 459 ! 460 USE utils 461 REAL(wp) :: zsmth ! smoothing around critical depth 462 REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space 463 ! 464 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 465 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1 466 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigt3, z_esigw3, z_esigtu3 467 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 468 469 ALLOCATE( z_gsigw3(jpi,jpj), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) ) 470 ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj)) 471 ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj)) 472 ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) ) 473 ALLOCATE( z_gsi3w3m1(jpi,jpj) ) 472 474 473 475 z_gsigw3 = 0. ; z_gsigt3 = 0. ; z_gsi3w3 = 0. … … 475 477 z_esigtu3 = 0. ; z_esigtv3 = 0. ; z_esigtf3 = 0. 476 478 z_esigwu3 = 0. ; z_esigwv3 = 0. 477 478 DO ji = 1, jpi 479 DO jj = 1, jpj 480 481 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 482 DO jk = 1, jpk 483 z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5, rn_bb ) 484 z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) 485 END DO 486 ELSE ! shallow water, uniform sigma 487 DO jk = 1, jpk 488 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp) 489 z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5 ) / REAL(jpk-1,wp) 490 END DO 491 ENDIF 492 ! 493 DO jk = 1, jpk-1 494 z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk) 495 z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 496 END DO 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) ) 499 ! 500 ! Coefficients for vertical depth as the sum of e3w scale factors 501 z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1) 502 DO jk = 2, jpk 503 z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 504 END DO 505 ! 506 DO jk = 1, jpk 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 ) 511 gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 512 END DO 513 ! 514 END DO ! for all jj's 515 END DO ! for all ji's 516 517 DO ji = 1, jpim1 518 DO jj = 1, jpjm1 519 DO jk = 1, jpk 520 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 521 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 522 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 523 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 524 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) & 525 & + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 526 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 527 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 528 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 529 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 530 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 479 480 481 DO jk = 1,jpk 482 DO ji = 1, jpi 483 DO jj = 1, jpj 484 IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma 485 486 zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,. 487 ! could be changed by users but care must be taken to do so carefully 488 zzb = 1.0-(zzb/hbatt(ji,jj)) 489 490 zzs = rn_zs / hbatt(ji,jj) 491 492 IF (rn_efold /= 0.0) THEN 493 zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) 494 ELSE 495 zsmth = 1.0 496 ENDIF 497 498 z_gsigw3(ji,jj) = REAL(jk-1,wp) /REAL(jpk-1,wp) 499 z_gsigw3p1(ji,jj) = REAL(jk,wp) /REAL(jpk-1,wp) 500 z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 501 z_gsigw3(ji,jj) = fgamma( z_gsigw3(ji,jj), zzb, zzs, zsmth ) 502 z_gsigw3p1(ji,jj) = fgamma( z_gsigw3p1(ji,jj), zzb, zzs, zsmth ) 503 z_gsigt3(ji,jj) = fgamma( z_gsigt3(ji,jj), zzb, zzs, zsmth ) 504 505 ELSE IF(ln_sigcrit) THEN ! shallow water, uniform sigma 506 507 z_gsigw3(ji,jj) = REAL(jk-1,wp) /REAL(jpk-1,wp) 508 z_gsigw3p1(ji,jj) = REAL(jk,wp) /REAL(jpk-1,wp) 509 z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 510 511 ELSE ! shallow water, z coordinates 512 513 z_gsigw3(ji,jj) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 514 z_gsigw3p1(ji,jj) = REAL(jk,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 515 z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 516 517 ENDIF 518 519 !gsi3w3m1 & z_gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk 520 IF( jk .EQ. 1) THEN 521 z_esigw3(ji,jj ) = 2.0 * (z_gsigt3(ji,jj ) - z_gsigw3(ji,jj )) 522 z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj) 523 ELSE 524 z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj) 525 z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj) 526 ENDIF 527 IF( jk .EQ. jpk) THEN 528 z_esigt3(ji,jj) = 2.0 * (z_gsigt3(ji,jj) - z_gsigw3(ji,jj)) 529 ELSE 530 z_esigt3(ji,jj) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj) 531 ENDIF 532 533 gdept_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj) 534 gdepw_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj) 535 gdep3w_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj) 536 537 ENDDO ! for all jj's 538 ENDDO ! for all ji's 539 540 DO ji=1,jpi-1 541 DO jj=1,jpj-1 542 543 z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) ) / & 544 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 545 z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) ) / & 546 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 547 z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) + & 548 hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) ) / & 549 ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 550 z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) ) / & 551 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 552 z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) ) / & 553 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 554 555 e3t_0(ji,jj)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj) 556 e3u_0(ji,jj)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj) 557 e3v_0(ji,jj)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj) 558 e3f_0(ji,jj)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj) 531 559 ! 532 e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) 533 e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) 534 e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) 535 e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) 536 ! 537 e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) 538 e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) 539 e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) ) 540 END DO 541 END DO 542 END DO 543 544 DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3) 545 DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 546 547 END SUBROUTINE s_sh94 548 549 SUBROUTINE s_sf12 550 551 !!---------------------------------------------------------------------- 552 !! *** ROUTINE s_sf12 *** 553 !! 554 !! ** Purpose : stretch the s-coordinate system 555 !! 556 !! ** Method : s-coordinate stretch using the Siddorn and Furner 2012? 557 !! mixed S/sigma/Z coordinate 558 !! 559 !! This method allows the maintenance of fixed surface and or 560 !! bottom cell resolutions (cf. geopotential coordinates) 561 !! within an analytically derived stretched S-coordinate framework. 562 !! 563 !! 564 !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 565 !!---------------------------------------------------------------------- 566 ! 567 USE utils 568 REAL(wp) :: zsmth ! smoothing around critical depth 569 REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space 570 ! 571 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 572 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3 573 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 574 ! 575 ALLOCATE( z_gsigw3(jpi,jpj,jpk), z_gsigt3(jpi,jpj,jpk), z_gsi3w3(jpi,jpj,jpk) ) 576 ALLOCATE( z_esigt3(jpi,jpj,jpk), z_esigw3(jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk)) 577 ALLOCATE( z_esigtv3(jpi,jpj,jpk), z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk)) 578 ALLOCATE( z_esigwv3(jpi,jpj,jpk) ) 579 580 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. 584 585 DO ji = 1, jpi 586 DO jj = 1, jpj 587 588 IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma 589 590 zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,. 591 ! could be changed by users but care must be taken to do so carefully 592 zzb = 1.0-(zzb/hbatt(ji,jj)) 593 594 zzs = rn_zs / hbatt(ji,jj) 595 596 IF (rn_efold /= 0.0) THEN 597 zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) 598 ELSE 599 zsmth = 1.0 600 ENDIF 601 602 DO jk = 1, jpk 603 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) 604 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 605 ENDDO 606 z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth ) 607 z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth ) 560 e3w_0(ji,jj)=hbatt(ji,jj)*z_esigw3(ji,jj) 561 e3uw_0(ji,jj)=hbatu(ji,jj)*z_esigwu3(ji,jj) 562 e3vw_0(ji,jj)=hbatv(ji,jj)*z_esigwv3(ji,jj) 563 ENDDO 564 ENDDO 565 ! Keep some arrays for next level 566 z_gsigt3m1 = z_gsigt3 567 z_gsiw3m1 = z_gsiw3 608 568 609 ELSE IF(ln_sigcrit) THEN ! shallow water, uniform sigma 610 611 DO jk = 1, jpk 612 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) 613 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 614 END DO 615 616 ELSE ! shallow water, z coordinates 617 618 DO jk = 1, jpk 619 z_gsigw3(ji,jj,jk) = REAL(jk-1,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)) 621 END DO 622 623 ENDIF 624 625 DO jk = 1, jpk-1 626 z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk) 627 z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 628 END DO 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)) 631 632 ! Coefficients for vertical depth as the sum of e3w scale factors 633 z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1) 634 DO jk = 2, jpk 635 z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 636 END DO 637 638 DO jk = 1, jpk 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) 641 gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 642 END DO 643 644 ENDDO ! for all jj's 645 ENDDO ! for all ji's 646 647 DO ji=1,jpi-1 648 DO jj=1,jpj-1 649 650 DO jk = 1, jpk 651 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / & 652 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 653 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / & 654 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 655 z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) + & 656 hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / & 657 ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 658 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / & 659 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 660 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / & 661 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 662 663 e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) 664 e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk) 665 e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk) 666 e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk) 667 ! 668 e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk) 669 e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk) 670 e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk) 671 END DO 672 673 ENDDO 674 ENDDO 675 ! 676 ! ! ============= 677 678 DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3) 569 where (e3t_0 (:,:).eq.0.0) e3t_0(:,:) = 1.0 570 where (e3u_0 (:,:).eq.0.0) e3u_0(:,:) = 1.0 571 where (e3v_0 (:,:).eq.0.0) e3v_0(:,:) = 1.0 572 where (e3f_0 (:,:).eq.0.0) e3f_0(:,:) = 1.0 573 where (e3w_0 (:,:).eq.0.0) e3w_0(:,:) = 1.0 574 where (e3uw_0 (:,:).eq.0.0) e3uw_0(:,:) = 1.0 575 where (e3vw_0 (:,:).eq.0.0) e3vw_0(:,:) = 1.0 576 577 WRITE(*,*) 'Writing level ',jk,' to file' 578 CALL write_netcdf_vars(jk) 579 WRITE(*,*) 'Written level ',jk,' to file' 580 ENDDO ! End of loop over jk 581 582 DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1) 679 583 DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 680 584 … … 728 632 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 729 633 DO jk = 1, jpk 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 )734 gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )735 634 END DO 736 635 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 737 DO jj = 1, jpj 738 DO ji = 1, jpi 739 DO jk = 1, jpk 740 e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) 741 e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpk-1,wp) ) 742 e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpk-1,wp) ) 743 e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpk-1,wp) ) 636 DO jk = 1, jpk 637 DO jj = 1, jpj 638 DO ji = 1, jpi 639 zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp) 640 zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp) 641 gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsigt(jk) + hift(ji,jj)*zcoeft ) 642 gdepw_0(:,:) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsigw(jk) + hift(ji,jj)*zcoefw ) 643 gdep3w_0(:,:) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsi3w(jk) + hift(ji,jj)*zcoeft ) 644 e3t_0(ji,jj) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) 645 e3u_0(ji,jj) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpk-1,wp) ) 646 e3v_0(ji,jj) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpk-1,wp) ) 647 e3f_0(ji,jj) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpk-1,wp) ) 744 648 ! 745 e3w_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) 746 e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpk-1,wp) ) 747 e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpk-1,wp) ) 649 e3w_0(ji,jj) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) ) 650 e3uw_0(ji,jj) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpk-1,wp) ) 651 e3vw_0(ji,jj) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpk-1,wp) ) 652 IF( scobot(ji,jj) >= gdept_0(ji,jj) ) mbathy(ji,jj) = MAX( 2, jk ) 653 IF( scobot(ji,jj) == 0. ) mbathy(ji,jj) = 0 748 654 END DO 749 655 END DO 750 END DO 656 where (e3t_0 (:,:).eq.0.0) e3t_0(:,:) = 1.0 657 where (e3u_0 (:,:).eq.0.0) e3u_0(:,:) = 1.0 658 where (e3v_0 (:,:).eq.0.0) e3v_0(:,:) = 1.0 659 where (e3f_0 (:,:).eq.0.0) e3f_0(:,:) = 1.0 660 where (e3w_0 (:,:).eq.0.0) e3w_0(:,:) = 1.0 661 where (e3uw_0 (:,:).eq.0.0) e3uw_0(:,:) = 1.0 662 where (e3vw_0 (:,:).eq.0.0) e3vw_0(:,:) = 1.0 663 664 CALL write_netcdf_vars(jk) 665 ENDDO ! End of loop over jk 666 751 667 752 668 DEALLOCATE( z_gsigw, z_gsigt, z_gsi3w ) … … 826 742 !!---------------------------------------------------------------------- 827 743 USE utils, ONLY : jpk,jk,wp 828 REAL(wp), INTENT(in ) :: pk1 (jpk)! continuous "k" coordinate829 REAL(wp) :: p_gamma (jpk)! stretched coordinate744 REAL(wp), INTENT(in ) :: pk1 ! continuous "k" coordinate 745 REAL(wp) :: p_gamma ! stretched coordinate 830 746 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth 831 747 REAL(wp), INTENT(in ) :: pzs ! surface box depth 832 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter833 REAL(wp) :: za1,za2,za3 834 REAL(wp) :: zn1,zn2 835 REAL(wp) :: za,zb,zx 748 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter 749 REAL(wp) :: za1,za2,za3 ! local variables 750 REAL(wp) :: zn1,zn2 ! local variables 751 REAL(wp) :: za,zb,zx ! local variables 836 752 !!---------------------------------------------------------------------- 837 753 ! … … 849 765 zx = 1.0-za/2.0-zb 850 766 851 DO jk = 1, jpk 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) 856 ENDDO 767 p_gamma = za*(pk1*(1.0-pk1/2.0))+zb*pk1**3.0 + & 768 & zx*( (rn_alpha+2.0)*pk1**(rn_alpha+1.0)- & 769 & (rn_alpha+1.0)*pk1**(rn_alpha+2.0) ) 770 p_gamma = p_gamma*psmth+pk1*(1.0-psmth) 857 771 858 772 ! -
branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/utils.F90
r5257 r5269 13 13 !! All coordinates 14 14 !! --------------- 15 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,: ,:) :: gdep3w_0 !: depth of t-points (sum of e3w) (m)16 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,: ,:) :: gdept_0, gdepw_0 !: analytical (time invariant) depth at t-w points (m)17 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,: ,:) :: e3v_0 , e3f_0 !: analytical (time invariant) vertical scale factors at v-f18 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,: ,:) :: e3t_0 , e3u_0 !: t-u points (m)19 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,: ,:) :: e3vw_0 !: analytical (time invariant) vertical scale factors at vw20 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,: ,:) :: e3w_0 , e3uw_0 !: w-uw points (m)15 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gdep3w_0 !: depth of t-points (sum of e3w) (m) 16 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: gdept_0, gdepw_0 !: analytical (time invariant) depth at t-w points (m) 17 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3v_0 , e3f_0 !: analytical (time invariant) vertical scale factors at v-f 18 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3t_0 , e3u_0 !: t-u points (m) 19 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3vw_0 !: analytical (time invariant) vertical scale factors at vw 20 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: e3w_0 , e3uw_0 !: w-uw points (m) 21 21 22 22 !! s-coordinate and hybrid z-s-coordinate … … 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 namelist47 INTEGER,PARAMETER :: numnam=8 ! File handle for namelist 48 48 REAL(wp) :: zrmax, ztaper ! temporary scalars 49 49 REAL(wp) :: zrfact … … 60 60 rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 61 61 62 ! IDs for output netcdf file 63 INTEGER :: id_x, id_y, id_z 64 INTEGER :: ncout 65 INTEGER, DIMENSION(11) :: var_ids !Array to contain all variable IDs 66 62 67 CONTAINS 63 68 … … 71 76 & zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj), STAT=ierr(1) ) 72 77 ! 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) )78 ALLOCATE( gdep3w_0(jpi,jpj) , e3v_0(jpi,jpj) , e3f_0(jpi,jpj) , & 79 & gdept_0(jpi,jpj) , e3t_0(jpi,jpj) , e3u_0 (jpi,jpj) , & 80 & gdepw_0(jpi,jpj) , e3w_0(jpi,jpj) , e3vw_0(jpi,jpj) , e3uw_0(jpi,jpj) , STAT=ierr(2) ) 76 81 ! 77 82 ! … … 97 102 98 103 ! Find the size of the input bathymetry 99 CALL dimlen(ncin, ' x', jpi)100 CALL dimlen(ncin, ' y', jpj)104 CALL dimlen(ncin, 'lon', jpi) 105 CALL dimlen(ncin, 'lat', jpj) 101 106 102 107 ALLOCATE( bathy(jpi, jpj) ) 103 108 104 109 ! Read the bathymetry variable from file 105 CALL check_nf90( nf90_inq_varid( ncin, ' bathymetry', var_id ), 'Cannot get variable ID for bathymetry')110 CALL check_nf90( nf90_inq_varid( ncin, 'Bathymetry', var_id ), 'Cannot get variable ID for bathymetry') 106 111 CALL check_nf90( nf90_get_var( ncin, var_id, bathy, (/ 1,1 /), (/ jpi, jpj /) ) ) 107 112 … … 125 130 126 131 127 SUBROUTINE write_coord_file() 128 ! Write out variables to the a netcdf coordinates file 132 SUBROUTINE make_coord_file() 133 ! Create new coordinates file and define dimensions and variables ready for 134 ! writing 129 135 130 INTEGER :: id_x, id_y, id_z131 INTEGER :: ncout132 INTEGER, DIMENSION(11) :: var_ids !Array to contain all variable IDs133 136 134 137 !Create the file … … 141 144 ! 142 145 !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)) )146 CALL check_nf90( nf90_def_var(ncout, 'gdept_0', nf90_double, (/id_x, id_y,id_z/), var_ids(1)) ) 147 CALL check_nf90( nf90_def_var(ncout, 'gdepw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(2)) ) 148 CALL check_nf90( nf90_def_var(ncout, 'gdep3w_0', nf90_double, (/id_x, id_y,id_z/), var_ids(3)) ) 149 CALL check_nf90( nf90_def_var(ncout, 'e3f_0', nf90_double, (/id_x, id_y,id_z/), var_ids(4)) ) 150 CALL check_nf90( nf90_def_var(ncout, 'e3t_0', nf90_double, (/id_x, id_y,id_z/), var_ids(5)) ) 151 CALL check_nf90( nf90_def_var(ncout, 'e3u_0', nf90_double, (/id_x, id_y,id_z/), var_ids(6)) ) 152 CALL check_nf90( nf90_def_var(ncout, 'e3v_0', nf90_double, (/id_x, id_y,id_z/), var_ids(7)) ) 153 CALL check_nf90( nf90_def_var(ncout, 'e3w_0', nf90_double, (/id_x, id_y,id_z/), var_ids(8)) ) 154 CALL check_nf90( nf90_def_var(ncout, 'e3uw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(9)) ) 155 CALL check_nf90( nf90_def_var(ncout, 'e3vw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(10)) ) 156 CALL check_nf90( nf90_def_var(ncout, 'mbathy', nf90_double, (/id_x, id_y/), var_ids(11)) ) 154 157 155 158 ! End define mode 156 159 CALL check_nf90( nf90_enddef(ncout) ) 160 161 WRITE(*,*) 'Opened coord_zgr.nc file and defined variables' 157 162 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) ) 163 END SUBROUTINE make_coord_file 172 164 173 END SUBROUTINE write_coord_file 165 SUBROUTINE write_netcdf_vars(kk) 166 ! Write variables to the netcdf file at level kk 167 INTEGER, INTENT(in) :: kk 168 169 CALL check_nf90( nf90_put_var(ncout, var_ids(1), gdept_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 170 CALL check_nf90( nf90_put_var(ncout, var_ids(2), gdepw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 171 CALL check_nf90( nf90_put_var(ncout, var_ids(3), gdep3w_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 172 CALL check_nf90( nf90_put_var(ncout, var_ids(4), e3f_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 173 CALL check_nf90( nf90_put_var(ncout, var_ids(5), e3t_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 174 CALL check_nf90( nf90_put_var(ncout, var_ids(6), e3u_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 175 CALL check_nf90( nf90_put_var(ncout, var_ids(7), e3v_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 176 CALL check_nf90( nf90_put_var(ncout, var_ids(8), e3w_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 177 CALL check_nf90( nf90_put_var(ncout, var_ids(9), e3uw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 178 CALL check_nf90( nf90_put_var(ncout, var_ids(10), e3vw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) ) 179 180 END SUBROUTINE write_netcdf_vars 174 181 175 182 SUBROUTINE check_nf90( istat, message )
Note: See TracChangeset
for help on using the changeset viewer.
