- Timestamp:
- 2016-06-06T07:57:00+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
r6624 r6667 13 13 !! dom_wri : create and write mesh and mask file(s) 14 14 !! dom_uniq : identify unique point of a grid (TUVF) 15 !! dom_stiff : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate) 15 16 !!---------------------------------------------------------------------- 16 17 USE dom_oce ! ocean space and time domain … … 27 28 PUBLIC dom_wri ! routine called by inidom.F90 28 29 PUBLIC dom_wri_coordinate ! routine called by domhgr.F90 29 30 PUBLIC dom_stiff ! routine called by inidom.F90 31 30 32 !! * Substitutions 31 33 # include "vectopt_loop_substitute.h90" … … 114 116 !! domhgr, domzgr, and dommsk. Note: the file contain depends on 115 117 !! the vertical coord. used (z-coord, partial steps, s-coord) 116 !! MOD(n msh, 3) = 1 : 'mesh_mask.nc' file118 !! MOD(nn_msh, 3) = 1 : 'mesh_mask.nc' file 117 119 !! = 2 : 'mesh.nc' and mask.nc' files 118 120 !! = 0 : 'mesh_hgr.nc', 'mesh_zgr.nc' and … … 121 123 !! vertical coordinate. 122 124 !! 123 !! if n msh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]124 !! if 3 < n msh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays125 !! if nn_msh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw] 126 !! if 3 < nn_msh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays 125 127 !! corresponding to the depth of the bottom t- and w-points 126 !! if 6 < n msh <= 9: write 2D arrays corresponding to the depth and the128 !! if 6 < nn_msh <= 9: write 2D arrays corresponding to the depth and the 127 129 !! thickness (e3[tw]_ps) of the bottom points 128 130 !! … … 149 151 IF( nn_timing == 1 ) CALL timing_start('dom_wri') 150 152 ! 151 CALL wrk_alloc( jpi, jpj, zprt, zprw)152 CALL wrk_alloc( jpi, jpj, jpk,zdepu, zdepv )153 CALL wrk_alloc( jpi,jpj, zprt , zprw ) 154 CALL wrk_alloc( jpi,jpj,jpk, zdepu, zdepv ) 153 155 ! 154 156 IF(lwp) WRITE(numout,*) … … 162 164 clnam4 = 'mesh_zgr' ! filename (vertical mesh informations) 163 165 164 SELECT CASE ( MOD(n msh, 3) )166 SELECT CASE ( MOD(nn_msh, 3) ) 165 167 ! ! ============================ 166 168 CASE ( 1 ) ! create 'mesh_mask.nc' file … … 201 203 IF( ln_zps ) THEN ; izps = 1 ; ELSE ; izps = 0 ; ENDIF 202 204 IF( ln_sco ) THEN ; isco = 1 ; ELSE ; isco = 0 ; ENDIF 203 CALL iom_rstput( 0, 0, inum , 'ln_zco' , REAL( izco, wp), ktype = jp_i4 )204 CALL iom_rstput( 0, 0, inum , 'ln_zps' , REAL( izps, wp), ktype = jp_i4 )205 CALL iom_rstput( 0, 0, inum , 'ln_sco' , REAL( isco, wp), ktype = jp_i4 )205 CALL iom_rstput( 0, 0, inum2, 'ln_zco' , REAL( izco, wp), ktype = jp_i4 ) 206 CALL iom_rstput( 0, 0, inum2, 'ln_zps' , REAL( izps, wp), ktype = jp_i4 ) 207 CALL iom_rstput( 0, 0, inum2, 'ln_sco' , REAL( isco, wp), ktype = jp_i4 ) 206 208 ! ! ocean cavities under iceshelves 207 209 IF( ln_isfcav ) THEN ; icav = 1 ; ELSE ; icav = 0 ; ENDIF 208 CALL iom_rstput( 0, 0, inum , 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 )210 CALL iom_rstput( 0, 0, inum2, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 ) 209 211 210 212 ! ! masks (inum2) … … 280 282 281 283 IF( ln_sco ) THEN ! s-coordinate 282 CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt )283 CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu )284 CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv )285 CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf )286 !287 CALL iom_rstput( 0, 0, inum4, 'gsigt', gsigt ) ! ! scaling coef.288 CALL iom_rstput( 0, 0, inum4, 'gsigw', gsigw )289 CALL iom_rstput( 0, 0, inum4, 'gsi3w', gsi3w )290 CALL iom_rstput( 0, 0, inum4, 'esigt', esigt )291 CALL iom_rstput( 0, 0, inum4, 'esigw', esigw )292 !293 284 CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 ) ! ! scale factors 294 285 CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 ) 295 286 CALL iom_rstput( 0, 0, inum4, 'e3v_0', e3v_0 ) 296 287 CALL iom_rstput( 0, 0, inum4, 'e3w_0', e3w_0 ) 297 CALL iom_rstput( 0, 0, inum4, 'rx1', rx1 ) ! ! Max. grid stiffness ratio 288 ! 289 CALL dom_stiff( zprt ) 290 CALL iom_rstput( 0, 0, inum4, 'stiffness', zprt ) ! ! Max. grid stiffness ratio 298 291 ! 299 292 CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d ) ! ! stretched system … … 305 298 IF( ln_zps ) THEN ! z-coordinate - partial steps 306 299 ! 307 IF( n msh <= 6 ) THEN ! ! 3D vertical scale factors300 IF( nn_msh <= 6 ) THEN ! ! 3D vertical scale factors 308 301 CALL iom_rstput( 0, 0, inum4, 'e3t_0', e3t_0 ) 309 302 CALL iom_rstput( 0, 0, inum4, 'e3u_0', e3u_0 ) … … 321 314 END IF 322 315 ! 323 IF( n msh <= 3 ) THEN ! ! 3D depth316 IF( nn_msh <= 3 ) THEN ! ! 3D depth 324 317 CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r8 ) 325 318 DO jk = 1,jpk … … 362 355 ! ! close the files 363 356 ! ! ============================ 364 SELECT CASE ( MOD(n msh, 3) )357 SELECT CASE ( MOD(nn_msh, 3) ) 365 358 CASE ( 1 ) 366 359 CALL iom_close( inum0 ) … … 425 418 END SUBROUTINE dom_uniq 426 419 420 421 SUBROUTINE dom_stiff( px1 ) 422 !!---------------------------------------------------------------------- 423 !! *** ROUTINE dom_stiff *** 424 !! 425 !! ** Purpose : Diagnose maximum grid stiffness/hydrostatic consistency 426 !! 427 !! ** Method : Compute Haney (1991) hydrostatic condition ratio 428 !! Save the maximum in the vertical direction 429 !! (this number is only relevant in s-coordinates) 430 !! 431 !! Haney, 1991, J. Phys. Oceanogr., 21, 610-619. 432 !!---------------------------------------------------------------------- 433 REAL(wp), DIMENSION(:,:), INTENT(out), OPTIONAL :: px1 ! stiffness 434 ! 435 INTEGER :: ji, jj, jk 436 REAL(wp) :: zrxmax 437 REAL(wp), DIMENSION(4) :: zr1 438 REAL(wp), DIMENSION(jpi,jpj) :: zx1 439 !!---------------------------------------------------------------------- 440 zx1(:,:) = 0._wp 441 zrxmax = 0._wp 442 zr1(:) = 0._wp 443 ! 444 DO ji = 2, jpim1 445 DO jj = 2, jpjm1 446 DO jk = 1, jpkm1 447 !!gm remark: dk(gdepw) = e3t ===>>> possible simplification of the following calculation.... 448 !! especially since it is gde3w which is used to compute the pressure gradient 449 !! furthermore, I think gdept_0 should be used below instead of w point in the numerator 450 !! so that the ratio is computed at the same point (i.e. uw and vw) .... 451 zr1(1) = ABS( ( gdepw_0(ji ,jj,jk )-gdepw_0(ji-1,jj,jk ) & 452 & +gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) ) & 453 & / ( gdepw_0(ji ,jj,jk )+gdepw_0(ji-1,jj,jk ) & 454 & -gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall ) ) * umask(ji-1,jj,jk) 455 zr1(2) = ABS( ( gdepw_0(ji+1,jj,jk )-gdepw_0(ji ,jj,jk ) & 456 & +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) ) & 457 & / ( gdepw_0(ji+1,jj,jk )+gdepw_0(ji ,jj,jk ) & 458 & -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) + rsmall ) ) * umask(ji ,jj,jk) 459 zr1(3) = ABS( ( gdepw_0(ji,jj+1,jk )-gdepw_0(ji,jj ,jk ) & 460 & +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) ) & 461 & / ( gdepw_0(ji,jj+1,jk )+gdepw_0(ji,jj ,jk ) & 462 & -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) + rsmall ) ) * vmask(ji,jj ,jk) 463 zr1(4) = ABS( ( gdepw_0(ji,jj ,jk )-gdepw_0(ji,jj-1,jk ) & 464 & +gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) ) & 465 & / ( gdepw_0(ji,jj ,jk )+gdepw_0(ji,jj-1,jk ) & 466 & -gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall ) ) * vmask(ji,jj-1,jk) 467 zrxmax = MAXVAL( zr1(1:4) ) 468 zx1(ji,jj) = MAX( zx1(ji,jj) , zrxmax ) 469 END DO 470 END DO 471 END DO 472 CALL lbc_lnk( zx1, 'T', 1. ) 473 ! 474 IF( PRESENT( px1 ) ) px1 = zx1 475 ! 476 zrxmax = MAXVAL( zx1 ) 477 ! 478 IF( lk_mpp ) CALL mpp_max( zrxmax ) ! max over the global domain 479 ! 480 IF(lwp) THEN 481 WRITE(numout,*) 482 WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax 483 WRITE(numout,*) '~~~~~~~~~' 484 ENDIF 485 ! 486 END SUBROUTINE dom_stiff 487 427 488 !!====================================================================== 428 489 END MODULE domwri
Note: See TracChangeset
for help on using the changeset viewer.
