- 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/domain.F90
r6624 r6667 20 20 !! dom_nam : read and contral domain namelists 21 21 !! dom_ctl : control print for the ocean domain 22 !! dom_stiff : diagnose maximum grid stiffness/hydrostatic consistency (s-coordinate)23 22 !! cfg_wri : create the "domain_cfg.nc" file containing all required configuration information 24 23 !!---------------------------------------------------------------------- … … 67 66 !! and scale factors, and the coriolis factor 68 67 !! - dom_zgr: define the vertical coordinate and the bathymetry 69 !! - dom_wri: create the meshmask file if n msh=168 !! - dom_wri: create the meshmask file if nn_msh=1 70 69 !! - 1D configuration, move Coriolis, u and v at T-point 71 70 !!---------------------------------------------------------------------- 72 INTEGER :: jk ! dummy loop indices 73 INTEGER :: iconf = 0 ! local integers 74 CHARACTER (len=64) :: cform = "(A12, 3(A13, I7) )" 75 REAL(wp), POINTER, DIMENSION(:,:) :: z1_hu_0, z1_hv_0 71 INTEGER :: ji, jj, jk ! dummy loop indices 72 INTEGER :: iconf = 0 ! local integers 73 CHARACTER (len=64) :: cform = "(A12, 3(A13, I7))" 74 INTEGER, DIMENSION(jpi,jpj) :: ik_top, ik_bot ! top and bottom ocean level 75 REAL(wp), POINTER, DIMENSION(:,:) :: zht, z1_hu_0, z1_hv_0 76 76 !!---------------------------------------------------------------------- 77 77 ! 78 78 IF( nn_timing == 1 ) CALL timing_start('dom_init') 79 79 ! 80 IF(lwp) THEN 80 IF(lwp) THEN ! Ocean domain Parameters (control print) 81 81 WRITE(numout,*) 82 82 WRITE(numout,*) 'dom_init : domain initialization' 83 83 WRITE(numout,*) '~~~~~~~~' 84 ENDIF 85 ! 86 ! Ocean domain Parameters (control print) 87 ! ----------------------- 88 IF(lwp) THEN 84 ! 89 85 WRITE(numout,*) ' Domain info' 90 86 WRITE(numout,*) ' dimension of model' … … 100 96 WRITE(numout,*) ' lateral boundary of the Global domain : jperio = ', jperio 101 97 ENDIF 102 103 ! !== Reference coordinate system ==! 104 ! 105 CALL dom_nam ! read namelist ( namrun, namdom ) 106 CALL dom_clo ! Closed seas and lake 107 CALL dom_hgr ! Horizontal mesh 108 CALL dom_zgr ! Vertical mesh and bathymetry 109 CALL dom_msk ! Masks 110 ! 111 IF( ln_sco ) CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency 98 ! 99 ! !== Reference coordinate system ==! 100 ! 101 CALL dom_nam ! read namelist ( namrun, namdom ) 102 CALL dom_clo ! Closed seas and lake 103 CALL dom_hgr ! Horizontal mesh 104 CALL dom_zgr( ik_top, ik_bot ) ! Vertical mesh and bathymetry 105 CALL dom_msk( ik_top, ik_bot ) ! Masks 106 ! 107 IF( ln_sco ) CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency 108 ! 109 DO jj = 1, jpj ! depth of the iceshelves 110 DO ji = 1, jpj 111 risfdep(ji,jj) = gdepw_0(ji,jj,mikt(ji,jj)) 112 END DO 113 END DO 112 114 ! 113 115 ht_0(:,:) = e3t_0(:,:,1) * tmask(:,:,1) ! Reference ocean thickness … … 120 122 END DO 121 123 ! 122 ! 124 ! !== time varying part of coordinate system ==! 123 125 ! 124 126 IF( ln_linssh ) THEN ! Fix in time : set to the reference one for all … … 158 160 IF( lk_c1d ) CALL cor_c1d ! 1D configuration: Coriolis set at T-point 159 161 ! 160 IF( n msh > 0 .AND. .NOT. ln_iscpl ) CALL dom_wri ! Create a domain file161 IF( n msh > 0 .AND. ln_iscpl .AND. .NOT. ln_rstart ) CALL dom_wri ! Create a domain file162 IF( nn_msh > 0 .AND. .NOT. ln_iscpl ) CALL dom_wri ! Create a domain file 163 IF( nn_msh > 0 .AND. ln_iscpl .AND. .NOT. ln_rstart ) CALL dom_wri ! Create a domain file 162 164 IF( .NOT.ln_rstart ) CALL dom_ctl ! Domain control 163 165 ! … … 165 167 IF(lwp) THEN 166 168 WRITE(numout,*) 167 WRITE(numout,*) 'dom_init : end of domain initialization n msh=', nmsh169 WRITE(numout,*) 'dom_init : end of domain initialization nn_msh=', nn_msh 168 170 WRITE(numout,*) 169 171 ENDIF 170 172 ! 171 IF( nmsh == -1) CALL cfg_wri ! create the configuration file173 IF( ln_write_cfg ) CALL cfg_wri ! create the configuration file 172 174 ! 173 175 IF( nn_timing == 1 ) CALL timing_stop('dom_init') … … 192 194 & nn_stock, nn_write , ln_mskland , ln_clobber , nn_chunksz, nn_euler , & 193 195 & ln_cfmeta, ln_iscpl 194 NAMELIST/namdom/ nn_bathy, rn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh, rn_hmin, rn_isfhmin, & 195 & rn_atfp , rn_rdt , nn_closea , ln_crs , & 196 & ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, & 197 & ppa2, ppkth2, ppacr2 196 NAMELIST/namdom/ ln_linssh, nn_closea, nn_msh, rn_isfhmin, rn_rdt, rn_atfp, ln_crs 198 197 #if defined key_netcdf4 199 198 NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip … … 261 260 neuler = 0 262 261 ENDIF 263 264 262 ! ! control of output frequency 265 263 IF ( nstock == 0 .OR. nstock > nitend ) THEN … … 305 303 WRITE(numout,*) 306 304 WRITE(numout,*) ' Namelist namdom : space & time domain' 307 WRITE(numout,*) ' flag read/compute bathymetry nn_bathy = ', nn_bathy 308 WRITE(numout,*) ' Depth (if =0 bathy=jpkm1) rn_bathy = ', rn_bathy 309 WRITE(numout,*) ' min depth of the ocean (>0) or rn_hmin = ', rn_hmin 310 WRITE(numout,*) ' min number of ocean level (<0) ' 311 WRITE(numout,*) ' treshold to open the isf cavity rn_isfhmin = ', rn_isfhmin, ' (m)' 312 WRITE(numout,*) ' minimum thickness of partial rn_e3zps_min = ', rn_e3zps_min, ' (m)' 313 WRITE(numout,*) ' step level rn_e3zps_rat = ', rn_e3zps_rat 314 WRITE(numout,*) ' create mesh/mask file(s) nn_msh = ', nn_msh 305 WRITE(numout,*) ' linear free surface (=T) ln_linssh = ', ln_linssh 306 WRITE(numout,*) ' suppression of closed seas (=0) nn_closea = ', nn_closea 307 WRITE(numout,*) ' create mesh/mask file(s) nn_msh = ', nn_msh 315 308 WRITE(numout,*) ' = 0 no file created ' 316 309 WRITE(numout,*) ' = 1 mesh_mask ' 317 310 WRITE(numout,*) ' = 2 mesh and mask ' 318 311 WRITE(numout,*) ' = 3 mesh_hgr, msh_zgr and mask' 319 WRITE(numout,*) ' ocean time step rn_rdt = ', rn_rdt 320 WRITE(numout,*) ' asselin time filter parameter rn_atfp = ', rn_atfp 321 WRITE(numout,*) ' suppression of closed seas (=0) nn_closea = ', nn_closea 322 WRITE(numout,*) ' online coarsening of dynamical fields ln_crs = ', ln_crs 323 324 WRITE(numout,*) ' ORCA r4, r2 and r05 coefficients ppsur = ', ppsur 325 WRITE(numout,*) ' ppa0 = ', ppa0 326 WRITE(numout,*) ' ppa1 = ', ppa1 327 WRITE(numout,*) ' ppkth = ', ppkth 328 WRITE(numout,*) ' ppacr = ', ppacr 329 WRITE(numout,*) ' Minimum vertical spacing ppdzmin = ', ppdzmin 330 WRITE(numout,*) ' Maximum depth pphmax = ', pphmax 331 WRITE(numout,*) ' Use double tanf function for vertical coordinates ldbletanh = ', ldbletanh 332 WRITE(numout,*) ' Double tanh function parameters ppa2 = ', ppa2 333 WRITE(numout,*) ' ppkth2 = ', ppkth2 334 WRITE(numout,*) ' ppacr2 = ', ppacr2 312 WRITE(numout,*) ' treshold to open the isf cavity rn_isfhmin = ', rn_isfhmin, ' (m)' 313 WRITE(numout,*) ' ocean time step rn_rdt = ', rn_rdt 314 WRITE(numout,*) ' asselin time filter parameter rn_atfp = ', rn_atfp 315 WRITE(numout,*) ' online coarsening of dynamical fields ln_crs = ', ln_crs 335 316 ENDIF 336 317 337 318 call flush( numout ) 338 319 ! 339 ntopo = nn_bathy ! conversion DOCTOR names into model names (this should disappear soon) 340 e3zps_min = rn_e3zps_min 341 e3zps_rat = rn_e3zps_rat 342 nmsh = nn_msh 320 ! ! ! conversion DOCTOR names into model names (this should disappear soon) 343 321 atfp = rn_atfp 344 322 rdt = rn_rdt … … 373 351 snc4set%luse = .FALSE. ! No NetCDF 4 case 374 352 #endif 375 376 call flush( numout )377 378 353 ! 379 354 END SUBROUTINE dom_nam … … 403 378 ze1max = MAXVAL( e1t(:,:), mask = tmask_i(:,:) == 1._wp ) 404 379 ze2max = MAXVAL( e2t(:,:), mask = tmask_i(:,:) == 1._wp ) 405 380 ! 406 381 iloc = MINLOC( e1t(:,:), mask = tmask_i(:,:) == 1._wp ) 407 382 iimi1 = iloc(1) + nimpp - 1 … … 430 405 431 406 432 SUBROUTINE dom_stiff433 !!----------------------------------------------------------------------434 !! *** ROUTINE dom_stiff ***435 !!436 !! ** Purpose : Diagnose maximum grid stiffness/hydrostatic consistency437 !!438 !! ** Method : Compute Haney (1991) hydrostatic condition ratio439 !! Save the maximum in the vertical direction440 !! (this number is only relevant in s-coordinates)441 !!442 !! Haney, R. L., 1991: On the pressure gradient force443 !! over steep topography in sigma coordinate ocean models.444 !! J. Phys. Oceanogr., 21, 610-619.445 !!----------------------------------------------------------------------446 INTEGER :: ji, jj, jk447 REAL(wp) :: zrxmax448 REAL(wp), DIMENSION(4) :: zr1449 !!----------------------------------------------------------------------450 rx1(:,:) = 0._wp451 zrxmax = 0._wp452 zr1(:) = 0._wp453 !454 DO ji = 2, jpim1455 DO jj = 2, jpjm1456 DO jk = 1, jpkm1457 zr1(1) = ABS( ( gdepw_0(ji ,jj,jk )-gdepw_0(ji-1,jj,jk ) &458 & +gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) ) &459 & / ( gdepw_0(ji ,jj,jk )+gdepw_0(ji-1,jj,jk ) &460 & -gdepw_0(ji ,jj,jk+1)-gdepw_0(ji-1,jj,jk+1) + rsmall ) ) * umask(ji-1,jj,jk)461 zr1(2) = ABS( ( gdepw_0(ji+1,jj,jk )-gdepw_0(ji ,jj,jk ) &462 & +gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) ) &463 & / ( gdepw_0(ji+1,jj,jk )+gdepw_0(ji ,jj,jk ) &464 & -gdepw_0(ji+1,jj,jk+1)-gdepw_0(ji ,jj,jk+1) + rsmall ) ) * umask(ji ,jj,jk)465 zr1(3) = ABS( ( gdepw_0(ji,jj+1,jk )-gdepw_0(ji,jj ,jk ) &466 & +gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) ) &467 & / ( gdepw_0(ji,jj+1,jk )+gdepw_0(ji,jj ,jk ) &468 & -gdepw_0(ji,jj+1,jk+1)-gdepw_0(ji,jj ,jk+1) + rsmall ) ) * vmask(ji,jj ,jk)469 zr1(4) = ABS( ( gdepw_0(ji,jj ,jk )-gdepw_0(ji,jj-1,jk ) &470 & +gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) ) &471 & / ( gdepw_0(ji,jj ,jk )+gdepw_0(ji,jj-1,jk ) &472 & -gdepw_0(ji,jj ,jk+1)-gdepw_0(ji,jj-1,jk+1) + rsmall ) ) * vmask(ji,jj-1,jk)473 zrxmax = MAXVAL( zr1(1:4) )474 rx1(ji,jj) = MAX( rx1(ji,jj) , zrxmax )475 END DO476 END DO477 END DO478 CALL lbc_lnk( rx1, 'T', 1. )479 !480 zrxmax = MAXVAL( rx1 )481 !482 IF( lk_mpp ) CALL mpp_max( zrxmax ) ! max over the global domain483 !484 IF(lwp) THEN485 WRITE(numout,*)486 WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax487 WRITE(numout,*) '~~~~~~~~~'488 ENDIF489 !490 END SUBROUTINE dom_stiff491 492 493 407 SUBROUTINE cfg_wri 494 408 !!---------------------------------------------------------------------- … … 503 417 !! domhgr, domzgr, and dommsk. Note: the file contain depends on 504 418 !! the vertical coord. used (z-coord, partial steps, s-coord) 505 !! MOD(n msh, 3) = 1 : 'mesh_mask.nc' file419 !! MOD(nn_msh, 3) = 1 : 'mesh_mask.nc' file 506 420 !! = 2 : 'mesh.nc' and mask.nc' files 507 421 !! = 0 : 'mesh_hgr.nc', 'mesh_zgr.nc' and … … 510 424 !! vertical coordinate. 511 425 !! 512 !! if n msh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw]513 !! if 3 < n msh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays426 !! if nn_msh <= 3: write full 3D arrays for e3[tuvw] and gdep[tuvw] 427 !! if 3 < nn_msh <= 6: write full 3D arrays for e3[tuvw] and 2D arrays 514 428 !! corresponding to the depth of the bottom t- and w-points 515 !! if 6 < n msh <= 9: write 2D arrays corresponding to the depth and the429 !! if 6 < nn_msh <= 9: write 2D arrays corresponding to the depth and the 516 430 !! thickness (e3[tw]_ps) of the bottom points 517 431 !! … … 523 437 INTEGER :: inum ! temprary units for 'domain_cfg.nc' file 524 438 CHARACTER(len=21) :: clnam ! filename (mesh and mask informations) 439 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! workspace 525 440 !!---------------------------------------------------------------------- 526 441 ! … … 537 452 538 453 ! !== global domain size ==! 454 ! 539 455 CALL iom_rstput( 0, 0, inum, 'jpiglo', REAL( jpiglo, wp), ktype = jp_i4 ) 540 456 CALL iom_rstput( 0, 0, inum, 'jpjglo', REAL( jpjglo, wp), ktype = jp_i4 ) 541 457 CALL iom_rstput( 0, 0, inum, 'jpkglo', REAL( jpk , wp), ktype = jp_i4 ) 542 458 ! 543 459 ! !== domain characteristics ==! 460 ! 544 461 ! ! lateral boundary of the global domain 545 462 CALL iom_rstput( 0, 0, inum, 'jperio', REAL( jperio, wp), ktype = jp_i4 ) 463 ! 546 464 ! ! type of vertical coordinate 547 465 IF( ln_zco ) THEN ; izco = 1 ; ELSE ; izco = 0 ; ENDIF … … 551 469 CALL iom_rstput( 0, 0, inum, 'ln_zps' , REAL( izps, wp), ktype = jp_i4 ) 552 470 CALL iom_rstput( 0, 0, inum, 'ln_sco' , REAL( isco, wp), ktype = jp_i4 ) 471 ! 553 472 ! ! ocean cavities under iceshelves 554 473 IF( ln_isfcav ) THEN ; icav = 1 ; ELSE ; icav = 0 ; ENDIF 555 474 CALL iom_rstput( 0, 0, inum, 'ln_isfcav', REAL( icav, wp), ktype = jp_i4 ) 556 475 ! 557 476 ! !== horizontal mesh ! 558 477 ! … … 579 498 CALL iom_rstput( 0, 0, inum, 'ff_f' , ff_f , ktype = jp_r8 ) ! coriolis factor 580 499 CALL iom_rstput( 0, 0, inum, 'ff_t' , ff_t , ktype = jp_r8 ) 581 582 500 ! 583 501 ! !== vertical mesh - 3D mask ==! 584 502 ! … … 594 512 CALL iom_rstput( 0, 0, inum, 'e3u_0' , e3u_0 , ktype = jp_r8 ) 595 513 CALL iom_rstput( 0, 0, inum, 'e3v_0' , e3v_0 , ktype = jp_r8 ) 514 CALL iom_rstput( 0, 0, inum, 'e3f_0' , e3f_0 , ktype = jp_r8 ) 596 515 CALL iom_rstput( 0, 0, inum, 'e3w_0' , e3w_0 , ktype = jp_r8 ) 597 ! 598 CALL iom_rstput( 0, 0, inum, 'tmask' , tmask , ktype = jp_i1 ) ! masks (in bytes) 599 CALL iom_rstput( 0, 0, inum, 'umask' , umask , ktype = jp_i1 ) 600 CALL iom_rstput( 0, 0, inum, 'vmask' , vmask , ktype = jp_i1 ) 601 CALL iom_rstput( 0, 0, inum, 'fmask' , fmask , ktype = jp_i1 ) 602 603 !!gm Probably not required fields : 604 CALL iom_rstput( 0, 0, inum, 'bathy' , bathy , ktype = jp_r8 ) ! depth of the ocean at T-points 605 CALL iom_rstput( 0, 0, inum, 'mbathy' , REAL( mbathy , wp ) , ktype = jp_i4 ) ! nb of ocean T-points 606 CALL iom_rstput( 0, 0, inum, 'mbkt' , REAL( mbkt , wp ) , ktype = jp_i4 ) ! nb of ocean T-points 607 608 ! 609 CALL iom_rstput( 0, 0, inum, 'bathy' , risfdep , ktype = jp_r8 ) ! depth of the iceshelves at T-points 610 CALL iom_rstput( 0, 0, inum, 'misfdep', REAL( misfdep, wp ) , ktype = jp_i4 ) ! nb of ocean T-points (ISF) 611 612 !!gm end 613 614 !!gm ? 615 ! CALL iom_rstput( 0, 0, inum, 'hbatt', hbatt ) 616 ! CALL iom_rstput( 0, 0, inum, 'hbatu', hbatu ) 617 ! CALL iom_rstput( 0, 0, inum, 'hbatv', hbatv ) 618 ! CALL iom_rstput( 0, 0, inum, 'hbatf', hbatf ) 619 !!gm ? 620 !!gm ? 621 ! CALL iom_rstput( 0, 0, inum, 'rx1' , rx1 ) ! Max. grid stiffness ratio 622 !!gm ? 623 ! 624 ! DO jk = 1,jpk 625 ! DO jj = 1, jpjm1 626 ! DO ji = 1, fs_jpim1 ! vector opt. 627 ! zdepu(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj ,jk) ) 628 ! zdepv(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji ,jj+1,jk) ) 629 ! END DO 630 ! END DO 631 ! END DO 632 ! CALL lbc_lnk( zdepu, 'U', 1. ) ; CALL lbc_lnk( zdepv, 'V', 1. ) 633 ! CALL iom_rstput( 0, 0, inum, 'gdepu' , zdepu , ktype = jp_r8 ) 634 ! CALL iom_rstput( 0, 0, inum, 'gdepv' , zdepv , ktype = jp_r8 ) 635 516 CALL iom_rstput( 0, 0, inum, 'e3uw_0' , e3uw_0 , ktype = jp_r8 ) 517 CALL iom_rstput( 0, 0, inum, 'e3vw_0' , e3vw_0 , ktype = jp_r8 ) 518 ! 519 ! !== ocean top and bottom level ==! 520 ! 521 CALL iom_rstput( 0, 0, inum, 'bottom level' , REAL( mbkt, wp )*ssmask , ktype = jp_i4 ) ! nb of ocean T-points 522 CALL iom_rstput( 0, 0, inum, 'top level' , REAL( mikt, wp )*ssmask , ktype = jp_i4 ) ! nb of ocean T-points (ISF) 523 ! 524 IF( ln_sco ) THEN ! s-coordinate: store grid stiffness ratio (Not required anyway) 525 CALL dom_stiff( z2d ) 526 CALL iom_rstput( 0, 0, inum, 'stiffness', z2d ) ! ! Max. grid stiffness ratio 527 ENDIF 528 ! 636 529 ! ! ============================ 637 530 ! ! close the files
Note: See TracChangeset
for help on using the changeset viewer.
