- Timestamp:
- 2017-03-17T08:46:30+01:00 (7 years ago)
- Location:
- branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OFF_SRC
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OFF_SRC/domrea.F90
r7217 r7806 17 17 USE lib_mpp ! distributed memory computing library 18 18 19 USE iom 19 20 USE domstp ! domain: set the time-step 20 21 … … 73 74 74 75 CALL dom_nam ! read namelist ( namrun, namdom, namcla ) 76 CALL dom_msk ! Masks 77 CALL dom_hgr ! Horizontal grid 75 78 CALL dom_zgr ! Vertical mesh and bathymetry option 76 CALL dom_grd ! Create a domain file 77 78 ! 79 ! - ML - Used in dom_vvl_sf_nxt and lateral diffusion routines 80 ! but could be usefull in many other routines 79 ! 81 80 e12t (:,:) = e1t(:,:) * e2t(:,:) 82 81 e1e2t (:,:) = e1t(:,:) * e2t(:,:) … … 91 90 re1v_e2v(:,:) = e1v(:,:) / e2v(:,:) 92 91 ! 93 hu(:,:) = 0._wp ! Ocean depth at U- and V-points94 hv(:,:) = 0._wp95 DO jk = 1, jpk96 hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)97 hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)98 END DO99 ! ! Inverse of the local depth100 hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1)101 hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)102 103 92 CALL dom_stp ! Time step 104 CALL dom_msk ! Masks105 93 CALL dom_ctl ! Domain control 106 94 … … 178 166 nstocklist = nn_stocklist 179 167 nwrite = nn_write 180 181 182 168 ! ! control of output frequency 183 169 IF ( nstock == 0 .OR. nstock > nitend ) THEN … … 222 208 904 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdom in configuration namelist', lwp ) 223 209 IF(lwm) WRITE ( numond, namdom ) 210 224 211 225 212 IF(lwp) THEN … … 321 308 END SUBROUTINE dom_nam 322 309 310 SUBROUTINE dom_msk 311 !!--------------------------------------------------------------------- 312 !! *** ROUTINE dom_msk *** 313 !! ** Purpose : Read the NetCDF file(s) which contain(s) all the 314 !! ocean mask informations and defines the interior domain T-mask. 315 !! 316 !! ** Method : Read in a file all the arrays generated in routines 317 !! dommsk: 'mask.nc' file 318 !! The interior ocean/land mask is computed from tmask 319 !! setting to zero the duplicated row and lines due to 320 !! MPP exchange halos, est-west cyclic and north fold 321 !! boundary conditions. 322 !! 323 !! ** Action : tmask_i : interiorland/ocean mask at t-point 324 !! tpol : ??? 325 !!---------------------------------------------------------------------- 326 ! 327 INTEGER :: inum ! local integers 328 INTEGER :: ji, jj, jk ! dummy loop indices 329 INTEGER :: iif, iil, ijf, ijl ! local integers 330 REAL(wp), POINTER, DIMENSION(:,:) :: zmbk 331 ! 332 !!--------------------------------------------------------------------- 333 334 335 336 IF(lwp) WRITE(numout,*) 337 IF(lwp) WRITE(numout,*) 'dom_rea : read NetCDF mesh and mask information file(s)' 338 IF(lwp) WRITE(numout,*) '~~~~~~~' 339 340 CALL wrk_alloc( jpi, jpj, zmbk ) 341 zmbk(:,:) = 0._wp 342 343 IF(lwp) WRITE(numout,*) ' one file in "mesh_mask.nc" ' 344 CALL iom_open( 'mask', inum ) 345 346 ! ! masks (inum2) 347 CALL iom_get( inum, jpdom_data, 'tmask', tmask ) 348 CALL iom_get( inum, jpdom_data, 'umask', umask ) 349 CALL iom_get( inum, jpdom_data, 'vmask', vmask ) 350 CALL iom_get( inum, jpdom_data, 'fmask', fmask ) 351 352 CALL lbc_lnk( tmask, 'T', 1._wp ) ! Lateral boundary conditions 353 CALL lbc_lnk( umask, 'U', 1._wp ) 354 CALL lbc_lnk( vmask, 'V', 1._wp ) 355 CALL lbc_lnk( fmask, 'F', 1._wp ) 356 357 #if defined key_c1d 358 ! set umask and vmask equal tmask in 1D configuration 359 IF(lwp) WRITE(numout,*) 360 IF(lwp) WRITE(numout,*) '********** 1D configuration : set umask and vmask equal tmask ********' 361 IF(lwp) WRITE(numout,*) '********** ********' 362 363 umask(:,:,:) = tmask(:,:,:) 364 vmask(:,:,:) = tmask(:,:,:) 365 #endif 366 367 #if defined key_degrad 368 CALL iom_get( inum, jpdom_data, 'facvolt', facvol ) 369 #endif 370 371 CALL iom_get( inum, jpdom_data, 'mbathy', zmbk ) ! number of ocean t-points 372 mbathy (:,:) = INT( zmbk(:,:) ) 373 misfdep(:,:) = 1 ! ice shelf case not yet done 374 375 CALL zgr_bot_level ! mbk. arrays (deepest ocean t-, u- & v-points 376 377 ! ! ============================ 378 ! ! close the files 379 ! ! ============================ 380 381 ! 382 ! Interior domain mask (used for global sum) 383 ! -------------------- 384 ssmask(:,:) = tmask(:,:,1) 385 tmask_i(:,:) = tmask(:,:,1) 386 iif = jpreci ! thickness of exchange halos in i-axis 387 iil = nlci - jpreci + 1 388 ijf = jprecj ! thickness of exchange halos in j-axis 389 ijl = nlcj - jprecj + 1 390 ! 391 tmask_i( 1 :iif, : ) = 0._wp ! first columns 392 tmask_i(iil:jpi, : ) = 0._wp ! last columns (including mpp extra columns) 393 tmask_i( : , 1 :ijf) = 0._wp ! first rows 394 tmask_i( : ,ijl:jpj) = 0._wp ! last rows (including mpp extra rows) 395 ! 396 ! ! north fold mask 397 tpol(1:jpiglo) = 1._wp 398 ! 399 IF( jperio == 3 .OR. jperio == 4 ) tpol(jpiglo/2+1:jpiglo) = 0._wp ! T-point pivot 400 IF( jperio == 5 .OR. jperio == 6 ) tpol( 1 :jpiglo) = 0._wp ! F-point pivot 401 IF( jperio == 3 .OR. jperio == 4 ) THEN ! T-point pivot: only half of the nlcj-1 row 402 IF( mjg(ijl-1) == jpjglo-1 ) THEN 403 DO ji = iif+1, iil-1 404 tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji)) 405 END DO 406 ENDIF 407 ENDIF 408 ! 409 ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at 410 ! least 1 wet u point 411 DO jj = 1, jpjm1 412 DO ji = 1, fs_jpim1 ! vector loop 413 umask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:))) 414 vmask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:))) 415 END DO 416 DO ji = 1, jpim1 ! NO vector opt. 417 fmask_i(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) & 418 & * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 419 END DO 420 END DO 421 CALL lbc_lnk( umask_i, 'U', 1._wp ) ! Lateral boundary conditions 422 CALL lbc_lnk( vmask_i, 'V', 1._wp ) 423 CALL lbc_lnk( fmask_i, 'F', 1._wp ) 424 425 ! 3. Ocean/land mask at wu-, wv- and w points 426 !---------------------------------------------- 427 wmask (:,:,1) = tmask(:,:,1) ! ???????? 428 wumask(:,:,1) = umask(:,:,1) ! ???????? 429 wvmask(:,:,1) = vmask(:,:,1) ! ???????? 430 DO jk = 2, jpk 431 wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1) 432 wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1) 433 wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1) 434 END DO 435 ! 436 CALL wrk_dealloc( jpi, jpj, zmbk ) 437 ! 438 CALL iom_close( inum ) 439 ! 440 END SUBROUTINE dom_msk 441 442 SUBROUTINE zgr_bot_level 443 !!---------------------------------------------------------------------- 444 !! *** ROUTINE zgr_bot_level *** 445 !! 446 !! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays) 447 !! 448 !! ** Method : computes from mbathy with a minimum value of 1 over land 449 !! 450 !! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest 451 !! ocean level at t-, u- & v-points 452 !! (min value = 1 over land) 453 !!---------------------------------------------------------------------- 454 ! 455 INTEGER :: ji, jj ! dummy loop indices 456 REAL(wp), POINTER, DIMENSION(:,:) :: zmbk 457 !!---------------------------------------------------------------------- 458 459 ! 460 IF(lwp) WRITE(numout,*) 461 IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels ' 462 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~' 463 ! 464 CALL wrk_alloc( jpi, jpj, zmbk ) 465 ! 466 mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land) 467 mikt(:,:) = 1 ; miku(:,:) = 1; mikv(:,:) = 1; ! top k-index of T-level (=1 over open ocean; >1 beneath ice shelf) 468 ! ! bottom k-index of W-level = mbkt+1 469 DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level 470 DO ji = 1, jpim1 471 mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) ) 472 mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) ) 473 END DO 474 END DO 475 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 476 zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 477 zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 478 ! 479 CALL wrk_dealloc( jpi, jpj, zmbk ) 480 ! 481 END SUBROUTINE zgr_bot_level 482 483 SUBROUTINE dom_hgr 484 !!---------------------------------------------------------------------- 485 !! *** ROUTINE dom_hgr *** 486 !! 487 !! ** Purpose : Read the NetCDF file(s) which contain(s) all the 488 !! ocean horizontal mesh informations 489 !! 490 !! ** Method : Read in a file all the arrays generated in routines 491 !! domhgr: 'mesh_hgr.nc' file 492 !!---------------------------------------------------------------------- 493 !! 494 INTEGER :: ji, jj ! dummy loop indices 495 INTEGER :: inum ! local integers 496 !!---------------------------------------------------------------------- 497 498 IF(lwp) WRITE(numout,*) 499 IF(lwp) WRITE(numout,*) 'dom_grd_hgr : read NetCDF mesh and mask information file(s)' 500 IF(lwp) WRITE(numout,*) '~~~~~~~' 501 502 IF(lwp) WRITE(numout,*) ' one file in "mesh_mask.nc" ' 503 CALL iom_open( 'mesh_hgr', inum ) 504 505 ! ! horizontal mesh (inum3) 506 CALL iom_get( inum, jpdom_data, 'glamt', glamt ) 507 CALL iom_get( inum, jpdom_data, 'glamu', glamu ) 508 CALL iom_get( inum, jpdom_data, 'glamv', glamv ) 509 CALL iom_get( inum, jpdom_data, 'glamf', glamf ) 510 511 CALL iom_get( inum, jpdom_data, 'gphit', gphit ) 512 CALL iom_get( inum, jpdom_data, 'gphiu', gphiu ) 513 CALL iom_get( inum, jpdom_data, 'gphiv', gphiv ) 514 CALL iom_get( inum, jpdom_data, 'gphif', gphif ) 515 516 CALL iom_get( inum, jpdom_data, 'e1t', e1t ) 517 CALL iom_get( inum, jpdom_data, 'e1u', e1u ) 518 CALL iom_get( inum, jpdom_data, 'e1v', e1v ) 519 520 CALL iom_get( inum, jpdom_data, 'e2t', e2t ) 521 CALL iom_get( inum, jpdom_data, 'e2u', e2u ) 522 CALL iom_get( inum, jpdom_data, 'e2v', e2v ) 523 524 CALL iom_get( inum, jpdom_data, 'ff', ff ) 525 526 527 ! Control printing : Grid informations (if not restart) 528 ! ---------------- 529 530 IF(lwp .AND. .NOT.ln_rstart ) THEN 531 WRITE(numout,*) 532 WRITE(numout,*) ' longitude and e1 scale factors' 533 WRITE(numout,*) ' ------------------------------' 534 WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1), & 535 glamv(ji,1), glamf(ji,1), & 536 e1t(ji,1), e1u(ji,1), & 537 e1v(ji,1), ji = 1, jpi,10) 538 539 WRITE(numout,*) 540 WRITE(numout,*) ' latitude and e2 scale factors' 541 WRITE(numout,*) ' -----------------------------' 542 WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj), & 543 & gphiv(1,jj), gphif(1,jj), & 544 & e2t (1,jj), e2u (1,jj), & 545 & e2v (1,jj), jj = 1, jpj, 10 ) 546 ENDIF 547 548 ! ! ============================ 549 ! ! close the files 550 ! ! ============================ 551 CALL iom_close( inum ) 552 ! 553 9300 FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x, & 554 f19.10, 1x, f19.10, 1x, f19.10 ) 555 END SUBROUTINE dom_hgr 556 557 323 558 SUBROUTINE dom_zgr 324 559 !!---------------------------------------------------------------------- 325 560 !! *** ROUTINE dom_zgr *** 326 561 !! 327 !! ** Purpose : set the depth of model levels and the resulting 328 !! vertical scale factors. 562 !! ** Purpose : Read the NetCDF file(s) which contain(s) all the 563 !! ocean horizontal mesh informations and/or set the depth of model levels 564 !! and the resulting vertical scale factors. 329 565 !! 330 566 !! ** Method : - reference 1D vertical coordinate (gdep._1d, e3._1d) … … 338 574 !! ** Action : define gdep., e3., mbathy and bathy 339 575 !!---------------------------------------------------------------------- 340 INTEGER :: ioptio = 0 ! temporary integer 341 INTEGER :: ios 576 INTEGER :: ioptio = 0 ! temporary integer 577 INTEGER :: inum, ios 578 INTEGER :: ji, jj, jk, ik 579 REAL(wp) :: zrefdep 342 580 !! 343 581 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 582 REAL(wp), POINTER, DIMENSION(:,:) :: zprt, zprw 344 583 !!---------------------------------------------------------------------- 345 584 … … 372 611 IF ( ioptio == 33 ) CALL ctl_stop( ' isf cavity with off line module not yet done ' ) 373 612 374 END SUBROUTINE dom_zgr 375 376 SUBROUTINE dom_ctl 377 !!---------------------------------------------------------------------- 378 !! *** ROUTINE dom_ctl *** 379 !! 380 !! ** Purpose : Domain control. 381 !! 382 !! ** Method : compute and print extrema of masked scale factors 383 !! 384 !! History : 385 !! 8.5 ! 02-08 (G. Madec) Original code 386 !!---------------------------------------------------------------------- 387 !! * Local declarations 388 INTEGER :: iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2 389 INTEGER, DIMENSION(2) :: iloc ! 390 REAL(wp) :: ze1min, ze1max, ze2min, ze2max 391 !!---------------------------------------------------------------------- 392 393 ! Extrema of the scale factors 394 395 IF(lwp)WRITE(numout,*) 396 IF(lwp)WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors' 397 IF(lwp)WRITE(numout,*) '~~~~~~~' 398 399 IF (lk_mpp) THEN 400 CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 ) 401 CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 ) 402 CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 ) 403 CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 ) 613 IF(lwp) WRITE(numout,*) ' one file in "mesh_mask.nc" ' 614 CALL iom_open( 'mesh_zgr', inum ) 615 616 CALL iom_get( inum, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth 617 CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', gdepw_1d ) 618 IF( ln_zco .OR. ln_zps ) THEN 619 CALL iom_get( inum, jpdom_unknown, 'e3t_1d' , e3t_1d ) ! reference scale factors 620 CALL iom_get( inum, jpdom_unknown, 'e3w_1d' , e3w_1d ) 621 ENDIF 622 623 !!gm BUG in s-coordinate this does not work! 624 ! deepest/shallowest W level Above/Below ~10m 625 zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_1d) ) ! ref. depth with tolerance (10% of minimum layer thickness) 626 nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 627 nla10 = nlb10 - 1 ! deepest W level Above ~10m 628 !!gm end bug 629 630 IF(lwp) THEN 631 WRITE(numout,*) 632 WRITE(numout,*) ' Reference z-coordinate depth and scale factors:' 633 WRITE(numout, "(9x,' level gdept gdepw e3t e3w ')" ) 634 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk ) 635 ENDIF 636 637 DO jk = 1, jpk 638 IF( e3w_1d (jk) <= 0._wp .OR. e3t_1d (jk) <= 0._wp ) CALL ctl_stop( ' e3w_1d or e3t_1d =< 0 ' ) 639 IF( gdepw_1d(jk) < 0._wp .OR. gdept_1d(jk) < 0._wp ) CALL ctl_stop( ' gdepw_1d or gdept_1d < 0 ' ) 640 END DO 641 642 IF( lk_vvl ) THEN 643 CALL iom_get( inum, jpdom_data, 'e3t_0', e3t_0(:,:,:) ) 644 CALL iom_get( inum, jpdom_data, 'e3u_0', e3u_0(:,:,:) ) 645 CALL iom_get( inum, jpdom_data, 'e3v_0', e3v_0(:,:,:) ) 646 CALL iom_get( inum, jpdom_data, 'e3w_0', e3w_0(:,:,:) ) 647 CALL iom_get( inum, jpdom_data, 'gdept_0', gdept_0(:,:,:) ) 648 CALL iom_get( inum, jpdom_data, 'gdepw_0', gdepw_0(:,:,:) ) 649 ht_0(:,:) = 0.0_wp ! Reference ocean depth at T-points 650 DO jk = 1, jpk 651 ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) 652 END DO 404 653 ELSE 405 ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )406 ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )407 ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )408 ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )409 410 iloc = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )411 iimi1 = iloc(1) + nimpp - 1412 ijmi1 = iloc(2) + njmpp - 1413 iloc = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )414 iimi2 = iloc(1) + nimpp - 1415 ijmi2 = iloc(2) + njmpp - 1416 iloc = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 )417 iima1 = iloc(1) + nimpp - 1418 ijma1 = iloc(2) + njmpp - 1419 iloc = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 )420 iima2 = iloc(1) + nimpp - 1421 ijma2 = iloc(2) + njmpp - 1422 ENDIF423 424 IF(lwp) THEN425 WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1426 WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1427 WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2428 WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2429 ENDIF430 431 END SUBROUTINE dom_ctl432 433 SUBROUTINE dom_grd434 !!----------------------------------------------------------------------435 !! *** ROUTINE dom_grd ***436 !!437 !! ** Purpose : Read the NetCDF file(s) which contain(s) all the438 !! ocean domain informations (mesh and mask arrays). This (these)439 !! file(s) is (are) used for visualisation (SAXO software) and440 !! diagnostic computation.441 !!442 !! ** Method : Read in a file all the arrays generated in routines443 !! domhgr, domzgr, and dommsk. Note: the file contain depends on444 !! the vertical coord. used (z-coord, partial steps, s-coord)445 !! nmsh = 1 : 'mesh_mask.nc' file446 !! = 2 : 'mesh.nc' and mask.nc' files447 !! = 3 : 'mesh_hgr.nc', 'mesh_zgr.nc' and448 !! 'mask.nc' files449 !! For huge size domain, use option 2 or 3 depending on your450 !! vertical coordinate.451 !!452 !! ** input file :453 !! meshmask.nc : domain size, horizontal grid-point position,454 !! masks, depth and vertical scale factors455 !!----------------------------------------------------------------------456 USE iom457 !!458 INTEGER :: ji, jj, jk ! dummy loop indices459 INTEGER :: ik, inum0 , inum1 , inum2 , inum3 , inum4 ! local integers460 REAL(wp) :: zrefdep ! local real461 REAL(wp), POINTER, DIMENSION(:,:) :: zmbk, zprt, zprw462 !!----------------------------------------------------------------------463 464 IF(lwp) WRITE(numout,*)465 IF(lwp) WRITE(numout,*) 'dom_rea : read NetCDF mesh and mask information file(s)'466 IF(lwp) WRITE(numout,*) '~~~~~~~'467 468 CALL wrk_alloc( jpi, jpj, zmbk, zprt, zprw )469 470 zmbk(:,:) = 0._wp471 472 SELECT CASE (nmsh)473 ! ! ============================474 CASE ( 1 ) ! create 'mesh_mask.nc' file475 ! ! ============================476 477 IF(lwp) WRITE(numout,*) ' one file in "mesh_mask.nc" '478 CALL iom_open( 'mesh_mask', inum0 )479 480 inum2 = inum0 ! put all the informations481 inum3 = inum0 ! in unit inum0482 inum4 = inum0483 484 ! ! ============================485 CASE ( 2 ) ! create 'mesh.nc' and486 ! ! 'mask.nc' files487 ! ! ============================488 489 IF(lwp) WRITE(numout,*) ' two files in "mesh.nc" and "mask.nc" '490 CALL iom_open( 'mesh', inum1 )491 CALL iom_open( 'mask', inum2 )492 493 inum3 = inum1 ! put mesh informations494 inum4 = inum1 ! in unit inum1495 496 ! ! ============================497 CASE ( 3 ) ! create 'mesh_hgr.nc'498 ! ! 'mesh_zgr.nc' and499 ! ! 'mask.nc' files500 ! ! ============================501 502 IF(lwp) WRITE(numout,*) ' three files in "mesh_hgr.nc" , "mesh_zgr.nc" and "mask.nc" '503 CALL iom_open( 'mesh_hgr', inum3 ) ! create 'mesh_hgr.nc'504 CALL iom_open( 'mesh_zgr', inum4 ) ! create 'mesh_zgr.nc'505 CALL iom_open( 'mask' , inum2 ) ! create 'mask.nc'506 507 ! ! ===========================508 CASE DEFAULT ! return error509 ! ! mesh has to be provided510 ! ! ===========================511 CALL ctl_stop( ' OFFLINE mode requires the input mesh mask(s). ', &512 & ' Invalid nn_msh value in the namelist (0 is not allowed)' )513 514 END SELECT515 516 ! ! masks (inum2)517 CALL iom_get( inum2, jpdom_data, 'tmask', tmask )518 CALL iom_get( inum2, jpdom_data, 'umask', umask )519 CALL iom_get( inum2, jpdom_data, 'vmask', vmask )520 CALL iom_get( inum2, jpdom_data, 'fmask', fmask )521 522 CALL lbc_lnk( tmask, 'T', 1._wp ) ! Lateral boundary conditions523 CALL lbc_lnk( umask, 'U', 1._wp )524 CALL lbc_lnk( vmask, 'V', 1._wp )525 CALL lbc_lnk( fmask, 'F', 1._wp )526 527 #if defined key_c1d528 ! set umask and vmask equal tmask in 1D configuration529 IF(lwp) WRITE(numout,*)530 IF(lwp) WRITE(numout,*) '********** 1D configuration : set umask and vmask equal tmask ********'531 IF(lwp) WRITE(numout,*) '********** ********'532 533 umask(:,:,:) = tmask(:,:,:)534 vmask(:,:,:) = tmask(:,:,:)535 #endif536 537 #if defined key_degrad538 CALL iom_get( inum2, jpdom_data, 'facvolt', facvol )539 #endif540 541 ! ! horizontal mesh (inum3)542 CALL iom_get( inum3, jpdom_data, 'glamt', glamt )543 CALL iom_get( inum3, jpdom_data, 'glamu', glamu )544 CALL iom_get( inum3, jpdom_data, 'glamv', glamv )545 CALL iom_get( inum3, jpdom_data, 'glamf', glamf )546 547 CALL iom_get( inum3, jpdom_data, 'gphit', gphit )548 CALL iom_get( inum3, jpdom_data, 'gphiu', gphiu )549 CALL iom_get( inum3, jpdom_data, 'gphiv', gphiv )550 CALL iom_get( inum3, jpdom_data, 'gphif', gphif )551 552 CALL iom_get( inum3, jpdom_data, 'e1t', e1t )553 CALL iom_get( inum3, jpdom_data, 'e1u', e1u )554 CALL iom_get( inum3, jpdom_data, 'e1v', e1v )555 556 CALL iom_get( inum3, jpdom_data, 'e2t', e2t )557 CALL iom_get( inum3, jpdom_data, 'e2u', e2u )558 CALL iom_get( inum3, jpdom_data, 'e2v', e2v )559 560 CALL iom_get( inum3, jpdom_data, 'ff', ff )561 562 CALL iom_get( inum4, jpdom_data, 'mbathy', zmbk ) ! number of ocean t-points563 mbathy (:,:) = INT( zmbk(:,:) )564 misfdep(:,:) = 1 ! ice shelf case not yet done565 566 CALL zgr_bot_level ! mbk. arrays (deepest ocean t-, u- & v-points567 !568 654 IF( ln_sco ) THEN ! s-coordinate 569 CALL iom_get( inum 4, jpdom_data, 'hbatt', hbatt )570 CALL iom_get( inum 4, jpdom_data, 'hbatu', hbatu )571 CALL iom_get( inum 4, jpdom_data, 'hbatv', hbatv )572 CALL iom_get( inum 4, jpdom_data, 'hbatf', hbatf )655 CALL iom_get( inum, jpdom_data, 'hbatt', hbatt ) 656 CALL iom_get( inum, jpdom_data, 'hbatu', hbatu ) 657 CALL iom_get( inum, jpdom_data, 'hbatv', hbatv ) 658 CALL iom_get( inum, jpdom_data, 'hbatf', hbatf ) 573 659 574 CALL iom_get( inum4, jpdom_unknown, 'gsigt', gsigt ) ! scaling coef. 575 CALL iom_get( inum4, jpdom_unknown, 'gsigw', gsigw ) 576 CALL iom_get( inum4, jpdom_unknown, 'gsi3w', gsi3w ) 577 CALL iom_get( inum4, jpdom_unknown, 'esigt', esigt ) 578 CALL iom_get( inum4, jpdom_unknown, 'esigw', esigw ) 579 580 CALL iom_get( inum4, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) ! scale factors 581 CALL iom_get( inum4, jpdom_data, 'e3u_0', fse3u_n(:,:,:) ) 582 CALL iom_get( inum4, jpdom_data, 'e3v_0', fse3v_n(:,:,:) ) 583 CALL iom_get( inum4, jpdom_data, 'e3w_0', fse3w_n(:,:,:) ) 584 585 CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth 586 CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d ) 660 CALL iom_get( inum, jpdom_unknown, 'gsigt', gsigt ) ! scaling coef. 661 CALL iom_get( inum, jpdom_unknown, 'gsigw', gsigw ) 662 CALL iom_get( inum, jpdom_unknown, 'gsi3w', gsi3w ) 663 CALL iom_get( inum, jpdom_unknown, 'esigt', esigt ) 664 CALL iom_get( inum, jpdom_unknown, 'esigw', esigw ) 665 666 CALL iom_get( inum, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) ! scale factors 667 CALL iom_get( inum, jpdom_data, 'e3u_0', fse3u_n(:,:,:) ) 668 CALL iom_get( inum, jpdom_data, 'e3v_0', fse3v_n(:,:,:) ) 669 CALL iom_get( inum, jpdom_data, 'e3w_0', fse3w_n(:,:,:) ) 587 670 ENDIF 588 671 589 672 590 673 IF( ln_zps ) THEN ! z-coordinate - partial steps 591 CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! reference depth592 CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )593 CALL iom_get( inum4, jpdom_unknown, 'e3t_1d' , e3t_1d ) ! reference scale factors594 CALL iom_get( inum4, jpdom_unknown, 'e3w_1d' , e3w_1d )595 674 ! 596 IF( nmsh <= 6 ) THEN ! 3D vertical scale factors597 CALL iom_get( inum 4, jpdom_data, 'e3t_0', fse3t_n(:,:,:) )598 CALL iom_get( inum 4, jpdom_data, 'e3u_0', fse3u_n(:,:,:) )599 CALL iom_get( inum 4, jpdom_data, 'e3v_0', fse3v_n(:,:,:) )600 CALL iom_get( inum 4, jpdom_data, 'e3w_0', fse3w_n(:,:,:) )675 IF( iom_varid( inum, 'e3t_0', ldstop = .FALSE. ) > 0 ) THEN 676 CALL iom_get( inum, jpdom_data, 'e3t_0', fse3t_n(:,:,:) ) 677 CALL iom_get( inum, jpdom_data, 'e3u_0', fse3u_n(:,:,:) ) 678 CALL iom_get( inum, jpdom_data, 'e3v_0', fse3v_n(:,:,:) ) 679 CALL iom_get( inum, jpdom_data, 'e3w_0', fse3w_n(:,:,:) ) 601 680 ELSE ! 2D bottom scale factors 602 CALL iom_get( inum 4, jpdom_data, 'e3t_ps', e3tp )603 CALL iom_get( inum 4, jpdom_data, 'e3w_ps', e3wp )681 CALL iom_get( inum, jpdom_data, 'e3t_ps', e3tp ) 682 CALL iom_get( inum, jpdom_data, 'e3w_ps', e3wp ) 604 683 ! ! deduces the 3D scale factors 605 684 DO jk = 1, jpk … … 633 712 END IF 634 713 635 IF( iom_varid( inum 4, 'gdept_0', ldstop = .FALSE. ) > 0 ) THEN ! 3D depth of t- and w-level636 CALL iom_get( inum 4, jpdom_data, 'gdept_0', fsdept_n(:,:,:) )637 CALL iom_get( inum 4, jpdom_data, 'gdepw_0', fsdepw_n(:,:,:) )714 IF( iom_varid( inum, 'gdept_0', ldstop = .FALSE. ) > 0 ) THEN ! 3D depth of t- and w-level 715 CALL iom_get( inum, jpdom_data, 'gdept_0', fsdept_n(:,:,:) ) 716 CALL iom_get( inum, jpdom_data, 'gdepw_0', fsdepw_n(:,:,:) ) 638 717 ELSE ! 2D bottom depth 639 CALL iom_get( inum4, jpdom_data, 'hdept', zprt ) 640 CALL iom_get( inum4, jpdom_data, 'hdepw', zprw ) 718 CALL wrk_alloc( jpi, jpj, zprt, zprw ) 719 ! 720 CALL iom_get( inum, jpdom_data, 'hdept', zprt ) 721 CALL iom_get( inum, jpdom_data, 'hdepw', zprw ) 641 722 ! 642 723 DO jk = 1, jpk ! deduces the 3D depth … … 654 735 END DO 655 736 END DO 737 CALL wrk_dealloc( jpi, jpj, zprt, zprw ) 656 738 ENDIF 657 739 ! … … 659 741 660 742 IF( ln_zco ) THEN ! Vertical coordinates and scales factors 661 CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth662 CALL iom_get( inum4, jpdom_unknown, 'gdepw_1d', gdepw_1d )663 CALL iom_get( inum4, jpdom_unknown, 'e3t_1d' , e3t_1d )664 CALL iom_get( inum4, jpdom_unknown, 'e3w_1d' , e3w_1d )665 743 DO jk = 1, jpk 666 744 fse3t_n(:,:,jk) = e3t_1d(jk) ! set to the ref. factors … … 672 750 END DO 673 751 ENDIF 674 675 !!gm BUG in s-coordinate this does not work! 676 ! deepest/shallowest W level Above/Below ~10m 677 zrefdep = 10._wp - ( 0.1_wp * MINVAL(e3w_1d) ) ! ref. depth with tolerance (10% of minimum layer thickness) 678 nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 679 nla10 = nlb10 - 1 ! deepest W level Above ~10m 680 !!gm end bug 681 682 ! Control printing : Grid informations (if not restart) 683 ! ---------------- 684 685 IF(lwp .AND. .NOT.ln_rstart ) THEN 686 WRITE(numout,*) 687 WRITE(numout,*) ' longitude and e1 scale factors' 688 WRITE(numout,*) ' ------------------------------' 689 WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1), & 690 glamv(ji,1), glamf(ji,1), & 691 e1t(ji,1), e1u(ji,1), & 692 e1v(ji,1), ji = 1, jpi,10) 693 9300 FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x, & 694 f19.10, 1x, f19.10, 1x, f19.10 ) 695 696 WRITE(numout,*) 697 WRITE(numout,*) ' latitude and e2 scale factors' 698 WRITE(numout,*) ' -----------------------------' 699 WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj), & 700 & gphiv(1,jj), gphif(1,jj), & 701 & e2t (1,jj), e2u (1,jj), & 702 & e2v (1,jj), jj = 1, jpj, 10 ) 703 ENDIF 704 705 706 IF( nprint == 1 .AND. lwp ) THEN 707 WRITE(numout,*) ' e1u e2u ' 708 CALL prihre( e1u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 709 CALL prihre( e2u,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 710 WRITE(numout,*) ' e1v e2v ' 711 CALL prihre( e1v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 712 CALL prihre( e2v,jpi,jpj,jpi-5,jpi,1,jpj-5,jpj,1,0.,numout ) 713 ENDIF 714 715 IF(lwp) THEN 716 WRITE(numout,*) 717 WRITE(numout,*) ' Reference z-coordinate depth and scale factors:' 718 WRITE(numout, "(9x,' level gdept gdepw e3t e3w ')" ) 719 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk ) 720 ENDIF 721 722 DO jk = 1, jpk 723 IF( e3w_1d (jk) <= 0._wp .OR. e3t_1d (jk) <= 0._wp ) CALL ctl_stop( ' e3w_1d or e3t_1d =< 0 ' ) 724 IF( gdepw_1d(jk) < 0._wp .OR. gdept_1d(jk) < 0._wp ) CALL ctl_stop( ' gdepw_1d or gdept_1d < 0 ' ) 725 END DO 752 ! 753 ENDIF 726 754 ! ! ============================ 727 755 ! ! close the files 728 756 ! ! ============================ 729 SELECT CASE ( nmsh ) 730 CASE ( 1 ) 731 CALL iom_close( inum0 ) 732 CASE ( 2 ) 733 CALL iom_close( inum1 ) 734 CALL iom_close( inum2 ) 735 CASE ( 3 ) 736 CALL iom_close( inum2 ) 737 CALL iom_close( inum3 ) 738 CALL iom_close( inum4 ) 739 END SELECT 740 ! 741 CALL wrk_dealloc( jpi, jpj, zmbk, zprt, zprw ) 742 ! 743 END SUBROUTINE dom_grd 744 745 746 SUBROUTINE zgr_bot_level 747 !!---------------------------------------------------------------------- 748 !! *** ROUTINE zgr_bot_level *** 749 !! 750 !! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays) 751 !! 752 !! ** Method : computes from mbathy with a minimum value of 1 over land 753 !! 754 !! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest 755 !! ocean level at t-, u- & v-points 756 !! (min value = 1 over land) 757 !!---------------------------------------------------------------------- 758 ! 759 INTEGER :: ji, jj ! dummy loop indices 760 REAL(wp), POINTER, DIMENSION(:,:) :: zmbk 761 !!---------------------------------------------------------------------- 762 763 ! 764 IF(lwp) WRITE(numout,*) 765 IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels ' 766 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~' 767 ! 768 CALL wrk_alloc( jpi, jpj, zmbk ) 769 ! 770 mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land) 771 mikt(:,:) = 1 ; miku(:,:) = 1; mikv(:,:) = 1; ! top k-index of T-level (=1 over open ocean; >1 beneath ice shelf) 772 ! ! bottom k-index of W-level = mbkt+1 773 DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level 774 DO ji = 1, jpim1 775 mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) ) 776 mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) ) 777 END DO 778 END DO 779 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 780 zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 781 zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 782 ! 783 CALL wrk_dealloc( jpi, jpj, zmbk ) 784 ! 785 END SUBROUTINE zgr_bot_level 786 787 SUBROUTINE dom_msk 788 !!--------------------------------------------------------------------- 789 !! *** ROUTINE dom_msk *** 790 !! 791 !! ** Purpose : Off-line case: defines the interior domain T-mask. 792 !! 793 !! ** Method : The interior ocean/land mask is computed from tmask 794 !! setting to zero the duplicated row and lines due to 795 !! MPP exchange halos, est-west cyclic and north fold 796 !! boundary conditions. 797 !! 798 !! ** Action : tmask_i : interiorland/ocean mask at t-point 799 !! tpol : ??? 800 !!---------------------------------------------------------------------- 801 ! 802 INTEGER :: ji, jj, jk ! dummy loop indices 803 INTEGER :: iif, iil, ijf, ijl ! local integers 804 INTEGER, POINTER, DIMENSION(:,:) :: imsk 805 ! 806 !!--------------------------------------------------------------------- 807 808 CALL wrk_alloc( jpi, jpj, imsk ) 809 ! 810 ! Interior domain mask (used for global sum) 811 ! -------------------- 812 ssmask(:,:) = tmask(:,:,1) 813 tmask_i(:,:) = tmask(:,:,1) 814 iif = jpreci ! thickness of exchange halos in i-axis 815 iil = nlci - jpreci + 1 816 ijf = jprecj ! thickness of exchange halos in j-axis 817 ijl = nlcj - jprecj + 1 818 ! 819 tmask_i( 1 :iif, : ) = 0._wp ! first columns 820 tmask_i(iil:jpi, : ) = 0._wp ! last columns (including mpp extra columns) 821 tmask_i( : , 1 :ijf) = 0._wp ! first rows 822 tmask_i( : ,ijl:jpj) = 0._wp ! last rows (including mpp extra rows) 823 ! 824 ! ! north fold mask 825 tpol(1:jpiglo) = 1._wp 826 ! 827 IF( jperio == 3 .OR. jperio == 4 ) tpol(jpiglo/2+1:jpiglo) = 0._wp ! T-point pivot 828 IF( jperio == 5 .OR. jperio == 6 ) tpol( 1 :jpiglo) = 0._wp ! F-point pivot 829 IF( jperio == 3 .OR. jperio == 4 ) THEN ! T-point pivot: only half of the nlcj-1 row 830 IF( mjg(ijl-1) == jpjglo-1 ) THEN 831 DO ji = iif+1, iil-1 832 tmask_i(ji,ijl-1) = tmask_i(ji,ijl-1) * tpol(mig(ji)) 833 END DO 834 ENDIF 835 ENDIF 836 ! 837 ! (ISF) MIN(1,SUM(umask)) is here to check if you have effectively at 838 ! least 1 wet u point 839 DO jj = 1, jpjm1 840 DO ji = 1, fs_jpim1 ! vector loop 841 umask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji+1,jj ) * MIN(1._wp,SUM(umask(ji,jj,:))) 842 vmask_i(ji,jj) = ssmask(ji,jj) * ssmask(ji ,jj+1) * MIN(1._wp,SUM(vmask(ji,jj,:))) 843 END DO 844 DO ji = 1, jpim1 ! NO vector opt. 845 fmask_i(ji,jj) = ssmask(ji,jj ) * ssmask(ji+1,jj ) & 846 & * ssmask(ji,jj+1) * ssmask(ji+1,jj+1) * MIN(1._wp,SUM(fmask(ji,jj,:))) 847 END DO 848 END DO 849 CALL lbc_lnk( umask_i, 'U', 1._wp ) ! Lateral boundary conditions 850 CALL lbc_lnk( vmask_i, 'V', 1._wp ) 851 CALL lbc_lnk( fmask_i, 'F', 1._wp ) 852 853 ! 3. Ocean/land mask at wu-, wv- and w points 854 !---------------------------------------------- 855 wmask (:,:,1) = tmask(:,:,1) ! ???????? 856 wumask(:,:,1) = umask(:,:,1) ! ???????? 857 wvmask(:,:,1) = vmask(:,:,1) ! ???????? 858 DO jk=2,jpk 859 wmask (:,:,jk)=tmask(:,:,jk) * tmask(:,:,jk-1) 860 wumask(:,:,jk)=umask(:,:,jk) * umask(:,:,jk-1) 861 wvmask(:,:,jk)=vmask(:,:,jk) * vmask(:,:,jk-1) 862 END DO 863 ! 864 IF( nprint == 1 .AND. lwp ) THEN ! Control print 865 imsk(:,:) = INT( tmask_i(:,:) ) 866 WRITE(numout,*) ' tmask_i : ' 867 CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout) 868 WRITE (numout,*) 869 WRITE (numout,*) ' dommsk: tmask for each level' 870 WRITE (numout,*) ' ----------------------------' 871 DO jk = 1, jpk 872 imsk(:,:) = INT( tmask(:,:,jk) ) 873 WRITE(numout,*) 874 WRITE(numout,*) ' level = ',jk 875 CALL prihin( imsk(:,:), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 1, numout) 876 END DO 877 ENDIF 878 ! 879 CALL wrk_dealloc( jpi, jpj, imsk ) 880 ! 881 END SUBROUTINE dom_msk 757 CALL iom_close( inum ) 758 ! 759 ! 760 END SUBROUTINE dom_zgr 761 762 SUBROUTINE dom_ctl 763 !!---------------------------------------------------------------------- 764 !! *** ROUTINE dom_ctl *** 765 !! 766 !! ** Purpose : Domain control. 767 !! 768 !! ** Method : compute and print extrema of masked scale factors 769 !! 770 !! History : 771 !! 8.5 ! 02-08 (G. Madec) Original code 772 !!---------------------------------------------------------------------- 773 !! * Local declarations 774 INTEGER :: iimi1, ijmi1, iimi2, ijmi2, iima1, ijma1, iima2, ijma2 775 INTEGER, DIMENSION(2) :: iloc ! 776 REAL(wp) :: ze1min, ze1max, ze2min, ze2max 777 !!---------------------------------------------------------------------- 778 779 ! Extrema of the scale factors 780 781 IF(lwp)WRITE(numout,*) 782 IF(lwp)WRITE(numout,*) 'dom_ctl : extrema of the masked scale factors' 783 IF(lwp)WRITE(numout,*) '~~~~~~~' 784 785 IF (lk_mpp) THEN 786 CALL mpp_minloc( e1t(:,:), tmask(:,:,1), ze1min, iimi1,ijmi1 ) 787 CALL mpp_minloc( e2t(:,:), tmask(:,:,1), ze2min, iimi2,ijmi2 ) 788 CALL mpp_maxloc( e1t(:,:), tmask(:,:,1), ze1max, iima1,ijma1 ) 789 CALL mpp_maxloc( e2t(:,:), tmask(:,:,1), ze2max, iima2,ijma2 ) 790 ELSE 791 ze1min = MINVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 ) 792 ze2min = MINVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 ) 793 ze1max = MAXVAL( e1t(:,:), mask = tmask(:,:,1) == 1.e0 ) 794 ze2max = MAXVAL( e2t(:,:), mask = tmask(:,:,1) == 1.e0 ) 795 796 iloc = MINLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 ) 797 iimi1 = iloc(1) + nimpp - 1 798 ijmi1 = iloc(2) + njmpp - 1 799 iloc = MINLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 ) 800 iimi2 = iloc(1) + nimpp - 1 801 ijmi2 = iloc(2) + njmpp - 1 802 iloc = MAXLOC( e1t(:,:), mask = tmask(:,:,1) == 1.e0 ) 803 iima1 = iloc(1) + nimpp - 1 804 ijma1 = iloc(2) + njmpp - 1 805 iloc = MAXLOC( e2t(:,:), mask = tmask(:,:,1) == 1.e0 ) 806 iima2 = iloc(1) + nimpp - 1 807 ijma2 = iloc(2) + njmpp - 1 808 ENDIF 809 810 IF(lwp) THEN 811 WRITE(numout,"(14x,'e1t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1max, iima1, ijma1 812 WRITE(numout,"(14x,'e1t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze1min, iimi1, ijmi1 813 WRITE(numout,"(14x,'e2t maxi: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2max, iima2, ijma2 814 WRITE(numout,"(14x,'e2t mini: ',1f10.2,' at i = ',i5,' j= ',i5)") ze2min, iimi2, ijmi2 815 ENDIF 816 817 END SUBROUTINE dom_ctl 882 818 883 819 !!====================================================================== -
branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90
r7256 r7806 22 22 USE c1d ! 1D configuration: lk_c1d 23 23 USE dom_oce ! ocean domain: variables 24 USE domvvl ! variable volume 24 25 USE zdf_oce ! ocean vertical physics: variables 25 26 USE sbc_oce ! surface module: variables … … 28 29 USE trabbl ! active tracer: bottom boundary layer 29 30 USE ldfslp ! lateral diffusion: iso-neutral slopes 31 USE sbcrnf ! river runoffs 30 32 USE ldfeiv ! eddy induced velocity coef. 31 33 USE ldftra_oce ! ocean tracer lateral physics … … 39 41 USE prtctl ! print control 40 42 USE fldread ! read input fields 43 USE wrk_nemo ! Memory allocation 41 44 USE timing ! Timing 45 USE trc, ONLY : ln_rsttr, numrtr, numrtw, lrst_trc 42 46 43 47 IMPLICIT NONE … … 46 50 PUBLIC dta_dyn_init ! called by opa.F90 47 51 PUBLIC dta_dyn ! called by step.F90 48 49 CHARACTER(len=100) :: cn_dir !: Root directory for location of ssr files 50 LOGICAL :: ln_dynwzv !: vertical velocity read in a file (T) or computed from u/v (F) 51 LOGICAL :: ln_dynbbl !: bbl coef read in a file (T) or computed (F) 52 LOGICAL :: ln_degrad !: degradation option enabled or not 53 LOGICAL :: ln_dynrnf !: read runoff data in file (T) or set to zero (F) 54 55 INTEGER , PARAMETER :: jpfld = 21 ! maximum number of fields to read 52 PUBLIC dta_dyn_swp ! called by step.F90 53 54 CHARACTER(len=100) :: cn_dir !: Root directory for location of ssr files 55 LOGICAL :: ln_ssh_ini !: initial ssh from dyn file (T) or not (F) - ssh is then read from passive tracer restart 56 LOGICAL :: ln_dynrnf !: read runoff data in file (T) or set to zero (F) 57 LOGICAL :: ln_dynrnf_depth !: read runoff data in file (T) or set to zero (F) 58 REAL(wp) :: fwbcorr 59 60 61 INTEGER , PARAMETER :: jpfld = 20 ! maximum number of fields to read 56 62 INTEGER , SAVE :: jf_tem ! index of temperature 57 63 INTEGER , SAVE :: jf_sal ! index of salinity 58 INTEGER , SAVE :: jf_uwd ! index of u- wind59 INTEGER , SAVE :: jf_vwd ! index of v- wind60 INTEGER , SAVE :: jf_wwd ! index of w-wind64 INTEGER , SAVE :: jf_uwd ! index of u-transport 65 INTEGER , SAVE :: jf_vwd ! index of v-transport 66 INTEGER , SAVE :: jf_wwd ! index of v-transport 61 67 INTEGER , SAVE :: jf_avt ! index of Kz 62 68 INTEGER , SAVE :: jf_mld ! index of mixed layer deptht 63 69 INTEGER , SAVE :: jf_emp ! index of water flux 70 INTEGER , SAVE :: jf_empb ! index of water flux 64 71 INTEGER , SAVE :: jf_qsr ! index of solar radiation 65 72 INTEGER , SAVE :: jf_wnd ! index of wind speed 66 73 INTEGER , SAVE :: jf_ice ! index of sea ice cover 67 74 INTEGER , SAVE :: jf_rnf ! index of river runoff 75 INTEGER , SAVE :: jf_fmf ! index of downward salt flux 68 76 INTEGER , SAVE :: jf_ubl ! index of u-bbl coef 69 77 INTEGER , SAVE :: jf_vbl ! index of v-bbl coef 70 INTEGER , SAVE :: jf_ahu ! index of u-diffusivity coef 71 INTEGER , SAVE :: jf_ahv ! index of v-diffusivity coef 72 INTEGER , SAVE :: jf_ahw ! index of w-diffusivity coef 73 INTEGER , SAVE :: jf_eiu ! index of u-eiv 74 INTEGER , SAVE :: jf_eiv ! index of v-eiv 75 INTEGER , SAVE :: jf_eiw ! index of w-eiv 76 INTEGER , SAVE :: jf_fmf ! index of downward salt flux 77 78 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_dyn ! structure of input fields (file informations, fields read) 78 INTEGER , SAVE :: jf_div ! index of e3t 79 80 81 TYPE(FLD), ALLOCATABLE, SAVE, DIMENSION(:) :: sf_dyn ! structure of input fields (file informations, fields read) 79 82 ! ! 80 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wdta ! vertical velocity at 2 time step81 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,: ) :: wnow ! vertical velocity at 2 time step82 83 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uslpdta ! zonal isopycnal slopes 83 84 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vslpdta ! meridional isopycnal slopes 84 85 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpidta ! zonal diapycnal slopes 85 86 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpjdta ! meridional diapycnal slopes 86 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: uslpnow ! zonal isopycnal slopes 87 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: vslpnow ! meridional isopycnal slopes 88 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wslpinow ! zonal diapycnal slopes 89 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wslpjnow ! meridional diapycnal slopes 90 91 INTEGER :: nrecprev_tem , nrecprev_uwd 87 88 INTEGER, SAVE :: nprevrec, nsecdyn 92 89 93 90 !! * Substitutions … … 113 110 !!---------------------------------------------------------------------- 114 111 ! 115 USE oce, ONLY: zts => tsa 116 USE oce, ONLY: zuslp => ua , zvslp => va 117 USE oce, ONLY: zwslpi => rotb , zwslpj => rotn 118 USE oce, ONLY: zu => ub , zv => vb, zw => hdivb 119 ! 112 USE oce, ONLY: zhdivtr => ua 120 113 INTEGER, INTENT(in) :: kt ! ocean time-step index 121 ! 122 INTEGER :: ji, jj ! dummy loop indices 123 INTEGER :: isecsbc ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step 124 REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation 125 REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation 126 INTEGER :: iswap_tem, iswap_uwd ! 114 INTEGER :: ji, jj, jk 115 REAL(wp), POINTER, DIMENSION(:,:) :: zemp 116 ! 127 117 !!---------------------------------------------------------------------- 128 118 … … 130 120 IF( nn_timing == 1 ) CALL timing_start( 'dta_dyn') 131 121 ! 132 isecsbc = nsec_year + nsec1jan000 133 ! 134 IF( kt == nit000 ) THEN 135 nrecprev_tem = 0 136 nrecprev_uwd = 0 137 ! 138 CALL fld_read( kt, 1, sf_dyn ) !== read data at kt time step ==! 139 ! 140 IF( lk_ldfslp .AND. .NOT.lk_c1d .AND. sf_dyn(jf_tem)%ln_tint ) THEN ! Computes slopes (here avt is used as workspace) 141 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:) ! temperature 142 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:) ! salinity 143 avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:) ! vertical diffusive coef. 144 CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 145 uslpdta (:,:,:,1) = zuslp (:,:,:) 146 vslpdta (:,:,:,1) = zvslp (:,:,:) 147 wslpidta(:,:,:,1) = zwslpi(:,:,:) 148 wslpjdta(:,:,:,1) = zwslpj(:,:,:) 149 ENDIF 150 IF( ln_dynwzv .AND. sf_dyn(jf_uwd)%ln_tint ) THEN ! compute vertical velocity from u/v 151 zu(:,:,:) = sf_dyn(jf_uwd)%fdta(:,:,:,1) 152 zv(:,:,:) = sf_dyn(jf_vwd)%fdta(:,:,:,1) 153 CALL dta_dyn_wzv( zu, zv, zw ) 154 wdta(:,:,:,1) = zw(:,:,:) * tmask(:,:,:) 155 ENDIF 156 ELSE 157 nrecprev_tem = sf_dyn(jf_tem)%nrec_a(2) 158 nrecprev_uwd = sf_dyn(jf_uwd)%nrec_a(2) 159 ! 160 CALL fld_read( kt, 1, sf_dyn ) !== read data at kt time step ==! 161 ! 162 ENDIF 163 ! 164 IF( lk_ldfslp .AND. .NOT.lk_c1d ) THEN ! Computes slopes (here avt is used as workspace) 165 iswap_tem = 0 166 IF( kt /= nit000 .AND. ( sf_dyn(jf_tem)%nrec_a(2) - nrecprev_tem ) /= 0 ) iswap_tem = 1 167 IF( ( isecsbc > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap_tem == 1 ) .OR. kt == nit000 ) THEN ! read/update the after data 168 IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt 169 IF( sf_dyn(jf_tem)%ln_tint ) THEN ! time interpolation of data 170 IF( kt /= nit000 ) THEN 171 uslpdta (:,:,:,1) = uslpdta (:,:,:,2) ! swap the data 172 vslpdta (:,:,:,1) = vslpdta (:,:,:,2) 173 wslpidta(:,:,:,1) = wslpidta(:,:,:,2) 174 wslpjdta(:,:,:,1) = wslpjdta(:,:,:,2) 175 ENDIF 176 ! 177 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:) ! temperature 178 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity 179 avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef. 180 CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 181 ! 182 uslpdta (:,:,:,2) = zuslp (:,:,:) 183 vslpdta (:,:,:,2) = zvslp (:,:,:) 184 wslpidta(:,:,:,2) = zwslpi(:,:,:) 185 wslpjdta(:,:,:,2) = zwslpj(:,:,:) 186 ELSE 187 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) 188 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) 189 avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) 190 CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 191 uslpnow (:,:,:) = zuslp (:,:,:) 192 vslpnow (:,:,:) = zvslp (:,:,:) 193 wslpinow(:,:,:) = zwslpi(:,:,:) 194 wslpjnow(:,:,:) = zwslpj(:,:,:) 195 ENDIF 196 ENDIF 197 IF( sf_dyn(jf_tem)%ln_tint ) THEN 198 ztinta = REAL( isecsbc - sf_dyn(jf_tem)%nrec_b(2), wp ) & 199 & / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp ) 200 ztintb = 1. - ztinta 201 uslp (:,:,:) = ztintb * uslpdta (:,:,:,1) + ztinta * uslpdta (:,:,:,2) 202 vslp (:,:,:) = ztintb * vslpdta (:,:,:,1) + ztinta * vslpdta (:,:,:,2) 203 wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1) + ztinta * wslpidta(:,:,:,2) 204 wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1) + ztinta * wslpjdta(:,:,:,2) 205 ELSE 206 uslp (:,:,:) = uslpnow (:,:,:) 207 vslp (:,:,:) = vslpnow (:,:,:) 208 wslpi(:,:,:) = wslpinow(:,:,:) 209 wslpj(:,:,:) = wslpjnow(:,:,:) 210 ENDIF 211 ENDIF 212 ! 213 IF( ln_dynwzv ) THEN ! compute vertical velocity from u/v 214 iswap_uwd = 0 215 IF( kt /= nit000 .AND. ( sf_dyn(jf_uwd)%nrec_a(2) - nrecprev_uwd ) /= 0 ) iswap_uwd = 1 216 IF( ( isecsbc > sf_dyn(jf_uwd)%nrec_b(2) .AND. iswap_uwd == 1 ) .OR. kt == nit000 ) THEN ! read/update the after data 217 IF(lwp) WRITE(numout,*) ' Compute new vertical velocity at kt = ', kt 218 IF(lwp) WRITE(numout,*) 219 IF( sf_dyn(jf_uwd)%ln_tint ) THEN ! time interpolation of data 220 IF( kt /= nit000 ) THEN 221 wdta(:,:,:,1) = wdta(:,:,:,2) ! swap the data for initialisation 222 ENDIF 223 zu(:,:,:) = sf_dyn(jf_uwd)%fdta(:,:,:,2) 224 zv(:,:,:) = sf_dyn(jf_vwd)%fdta(:,:,:,2) 225 CALL dta_dyn_wzv( zu, zv, zw ) 226 wdta(:,:,:,2) = zw(:,:,:) * tmask(:,:,:) 227 ELSE 228 zu(:,:,:) = sf_dyn(jf_uwd)%fnow(:,:,:) 229 zv(:,:,:) = sf_dyn(jf_vwd)%fnow(:,:,:) 230 CALL dta_dyn_wzv( zu, zv, zw ) 231 wnow(:,:,:) = zw(:,:,:) * tmask(:,:,:) 232 ENDIF 233 ENDIF 234 IF( sf_dyn(jf_uwd)%ln_tint ) THEN 235 ztinta = REAL( isecsbc - sf_dyn(jf_uwd)%nrec_b(2), wp ) & 236 & / REAL( sf_dyn(jf_uwd)%nrec_a(2) - sf_dyn(jf_uwd)%nrec_b(2), wp ) 237 ztintb = 1. - ztinta 238 wn(:,:,:) = ztintb * wdta(:,:,:,1) + ztinta * wdta(:,:,:,2) 239 ELSE 240 wn(:,:,:) = wnow(:,:,:) 241 ENDIF 242 ENDIF 122 ! 123 nsecdyn = nsec_year + nsec1jan000 ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step 124 ! 125 IF( kt == nit000 ) THEN ; nprevrec = 0 126 ELSE ; nprevrec = sf_dyn(jf_tem)%nrec_a(2) 127 ENDIF 128 ! 129 CALL fld_read( kt, 1, sf_dyn ) != read data at kt time step ==! 130 ! 131 IF( lk_ldfslp .AND. .NOT.lk_c1d ) CALL dta_dyn_slp( kt ) ! Computation of slopes 243 132 ! 244 133 tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature 245 134 tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity 246 ! 135 wndm(:,:) = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1) ! wind speed - needed for gas exchange 136 fmmflx(:,:) = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1) ! downward salt flux (v3.5+) 137 fr_i(:,:) = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1) ! Sea-ice fraction 138 qsr (:,:) = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1) ! solar radiation 139 emp (:,:) = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1) ! E-P 140 IF( ln_dynrnf ) THEN 141 rnf (:,:) = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1) ! E-P 142 IF( ln_dynrnf_depth .AND. lk_vvl ) CALL dta_dyn_hrnf 143 ENDIF 144 ! 145 un(:,:,:) = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:) ! effective u-transport 146 vn(:,:,:) = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:) ! effective v-transport 147 wn(:,:,:) = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:) ! effective v-transport 148 ! 149 IF( lk_vvl ) THEN 150 CALL wrk_alloc(jpi, jpj, zemp ) 151 zhdivtr(:,:,:) = sf_dyn(jf_div)%fnow(:,:,:) * tmask(:,:,:) ! effective u-transport 152 emp_b (:,:) = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1) ! E-P 153 zemp(:,:) = 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr * tmask(:,:,1) 154 CALL dta_dyn_ssh( kt, zhdivtr, sshb, zemp, ssha, fse3t_a(:,:,:) ) != ssh, vertical scale factor & vertical transport 155 CALL wrk_dealloc(jpi, jpj, zemp ) 156 ! Write in the tracer restart file 157 ! ******************************* 158 IF( lrst_trc ) THEN 159 IF(lwp) WRITE(numout,*) 160 IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file ', & 161 & 'at it= ', kt,' date= ', ndastp 162 IF(lwp) WRITE(numout,*) '~~~~' 163 CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssha ) 164 CALL iom_rstput( kt, nitrst, numrtw, 'sshb', sshn ) 165 ENDIF 166 ENDIF 247 167 ! 248 168 CALL eos ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop … … 251 171 252 172 rn2b(:,:,:) = rn2(:,:,:) ! need for zdfmxl 253 CALL zdf_mxl( kt ) ! In any case, we need mxl 254 ! 255 avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coefficient 256 un (:,:,:) = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:) ! u-velocity 257 vn (:,:,:) = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:) ! v-velocity 258 IF( .NOT.ln_dynwzv ) & ! w-velocity read in file 259 wn (:,:,:) = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:) 260 hmld(:,:) = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1) ! mixed layer depht 261 wndm(:,:) = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1) ! wind speed - needed for gas exchange 262 emp (:,:) = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1) ! E-P 263 fmmflx(:,:) = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1) ! downward salt flux (v3.5+) 264 fr_i(:,:) = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1) ! Sea-ice fraction 265 qsr (:,:) = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1) ! solar radiation 266 IF( ln_dynrnf ) & 267 rnf (:,:) = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1) ! river runoffs 268 269 ! ! bbl diffusive coef 173 CALL zdf_mxl( kt ) ! In any case, we need mxl 174 ! 175 hmld(:,:) = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1) ! mixed layer depht 176 avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coefficient 177 ! 270 178 #if defined key_trabbl && ! defined key_c1d 271 IF( ln_dynbbl ) THEN ! read in a file 272 ahu_bbl(:,:) = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1) 273 ahv_bbl(:,:) = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1) 274 ELSE ! Compute bbl coefficients if needed 275 tsb(:,:,:,:) = tsn(:,:,:,:) 276 CALL bbl( kt, nit000, 'TRC') 277 END IF 179 ahu_bbl(:,:) = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1) ! bbl diffusive coef 180 ahv_bbl(:,:) = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1) 278 181 #endif 279 #if ( ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv ) && ! defined key_c1d 280 aeiw(:,:) = sf_dyn(jf_eiw)%fnow(:,:,1) * tmask(:,:,1) ! w-eiv 281 ! ! Computes the horizontal values from the vertical value 282 DO jj = 2, jpjm1 283 DO ji = fs_2, fs_jpim1 ! vector opt. 284 aeiu(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji+1,jj ) ) ! Average the diffusive coefficient at u- v- points 285 aeiv(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji ,jj+1) ) ! at u- v- points 286 END DO 287 END DO 288 CALL lbc_lnk( aeiu, 'U', 1. ) ; CALL lbc_lnk( aeiv, 'V', 1. ) ! lateral boundary condition 289 #endif 290 291 #if defined key_degrad && ! defined key_c1d 292 ! ! degrad option : diffusive and eiv coef are 3D 293 ahtu(:,:,:) = sf_dyn(jf_ahu)%fnow(:,:,:) * umask(:,:,:) 294 ahtv(:,:,:) = sf_dyn(jf_ahv)%fnow(:,:,:) * vmask(:,:,:) 295 ahtw(:,:,:) = sf_dyn(jf_ahw)%fnow(:,:,:) * tmask(:,:,:) 296 # if defined key_traldf_eiv 297 aeiu(:,:,:) = sf_dyn(jf_eiu)%fnow(:,:,:) * umask(:,:,:) 298 aeiv(:,:,:) = sf_dyn(jf_eiv)%fnow(:,:,:) * vmask(:,:,:) 299 aeiw(:,:,:) = sf_dyn(jf_eiw)%fnow(:,:,:) * tmask(:,:,:) 300 # endif 301 #endif 182 ! 183 ! 184 CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 302 185 ! 303 186 IF(ln_ctl) THEN ! print control … … 308 191 CALL prt_ctl(tab3d_1=wn , clinfo1=' wn - : ', mask1=tmask, ovlap=1, kdim=jpk ) 309 192 CALL prt_ctl(tab3d_1=avt , clinfo1=' kz - : ', mask1=tmask, ovlap=1, kdim=jpk ) 310 CALL prt_ctl(tab2d_1=fr_i , clinfo1=' fr_i - : ', mask1=tmask, ovlap=1 )311 CALL prt_ctl(tab2d_1=hmld , clinfo1=' hmld - : ', mask1=tmask, ovlap=1 )312 CALL prt_ctl(tab2d_1=fmmflx , clinfo1=' fmmflx - : ', mask1=tmask, ovlap=1 )313 CALL prt_ctl(tab2d_1=emp , clinfo1=' emp - : ', mask1=tmask, ovlap=1 )314 CALL prt_ctl(tab2d_1=wndm , clinfo1=' wspd - : ', mask1=tmask, ovlap=1 )315 CALL prt_ctl(tab2d_1=qsr , clinfo1=' qsr - : ', mask1=tmask, ovlap=1 )193 ! CALL prt_ctl(tab2d_1=fr_i , clinfo1=' fr_i - : ', mask1=tmask, ovlap=1 ) 194 ! CALL prt_ctl(tab2d_1=hmld , clinfo1=' hmld - : ', mask1=tmask, ovlap=1 ) 195 ! CALL prt_ctl(tab2d_1=fmmflx , clinfo1=' fmmflx - : ', mask1=tmask, ovlap=1 ) 196 ! CALL prt_ctl(tab2d_1=emp , clinfo1=' emp - : ', mask1=tmask, ovlap=1 ) 197 ! CALL prt_ctl(tab2d_1=wndm , clinfo1=' wspd - : ', mask1=tmask, ovlap=1 ) 198 ! CALL prt_ctl(tab2d_1=qsr , clinfo1=' qsr - : ', mask1=tmask, ovlap=1 ) 316 199 ENDIF 317 200 ! … … 335 218 INTEGER :: inum, idv, idimv ! local integer 336 219 INTEGER :: ios ! Local integer output status for namelist read 337 !! 338 CHARACTER(len=100) :: cn_dir ! Root directory for location of core files 339 TYPE(FLD_N), DIMENSION(jpfld) :: slf_d ! array of namelist informations on the fields to read 340 TYPE(FLD_N) :: sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd, sn_rnf ! informations about the fields to be read 341 TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl ! " " 342 TYPE(FLD_N) :: sn_ahu, sn_ahv, sn_ahw, sn_eiu, sn_eiv, sn_eiw, sn_fmf ! " " 343 !!---------------------------------------------------------------------- 344 ! 345 NAMELIST/namdta_dyn/cn_dir, ln_dynwzv, ln_dynbbl, ln_degrad, ln_dynrnf, & 346 & sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd, sn_rnf, & 347 & sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl, & 348 & sn_ahu, sn_ahv, sn_ahw, sn_eiu, sn_eiv, sn_eiw, sn_fmf 220 INTEGER :: ji, jj, jk 221 REAL(wp) :: zcoef 222 INTEGER :: nkrnf_max 223 REAL(wp) :: hrnf_max 224 !! 225 CHARACTER(len=100) :: cn_dir ! Root directory for location of core files 226 TYPE(FLD_N), DIMENSION(jpfld) :: slf_d ! array of namelist informations on the fields to read 227 TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_empb, sn_emp ! informations about the fields to be read 228 TYPE(FLD_N) :: sn_tem , sn_sal , sn_avt ! " " 229 TYPE(FLD_N) :: sn_mld, sn_qsr, sn_wnd , sn_ice , sn_fmf ! " " 230 TYPE(FLD_N) :: sn_ubl, sn_vbl, sn_rnf ! " " 231 TYPE(FLD_N) :: sn_div ! informations about the fields to be read 232 !!---------------------------------------------------------------------- 233 234 NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth, ln_ssh_ini, fwbcorr, & 235 & sn_uwd, sn_vwd, sn_wwd, sn_emp, & 236 & sn_avt, sn_tem, sn_sal, sn_mld , sn_qsr , & 237 & sn_wnd, sn_ice, sn_fmf, & 238 & sn_ubl, sn_vbl, sn_rnf, & 239 & sn_empb, sn_div 349 240 ! 350 241 REWIND( numnam_ref ) ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data … … 363 254 WRITE(numout,*) '~~~~~~~ ' 364 255 WRITE(numout,*) ' Namelist namdta_dyn' 365 WRITE(numout,*) ' vertical velocity read from file (T) or computed (F) ln_dynwzv = ', ln_dynwzv366 WRITE(numout,*) ' bbl coef read from file (T) or computed (F) ln_dynbbl = ', ln_dynbbl367 WRITE(numout,*) ' degradation option enabled (T) or not (F) ln_degrad = ', ln_degrad368 WRITE(numout,*) ' river runoff option enabled (T) or not (F) ln_dynrnf = ', ln_dynrnf256 WRITE(numout,*) ' ssh initialised from dyn file (T) or not (F) ln_ssh_ini = ', ln_ssh_ini 257 WRITE(numout,*) ' runoffs option enabled (T) or not (F) ln_dynrnf = ', ln_dynrnf 258 WRITE(numout,*) ' runoffs is spread in vertical ln_dynrnf_depth = ', ln_dynrnf_depth 259 WRITE(numout,*) ' annual global mean of empmr for ssh correction fwbcorr = ', fwbcorr 369 260 WRITE(numout,*) 370 261 ENDIF 371 262 ! 372 IF( ln_degrad .AND. .NOT.lk_degrad ) THEN 373 CALL ctl_warn( 'dta_dyn_init: degradation option requires key_degrad activated ; force ln_degrad to false' ) 374 ln_degrad = .FALSE. 375 ENDIF 376 IF( ln_dynbbl .AND. ( .NOT.lk_trabbl .OR. lk_c1d ) ) THEN 377 CALL ctl_warn( 'dta_dyn_init: bbl option requires key_trabbl activated ; force ln_dynbbl to false' ) 378 ln_dynbbl = .FALSE. 379 ENDIF 380 381 jf_tem = 1 ; jf_sal = 2 ; jf_mld = 3 ; jf_emp = 4 ; jf_fmf = 5 ; jf_ice = 6 ; jf_qsr = 7 382 jf_wnd = 8 ; jf_uwd = 9 ; jf_vwd = 10 ; jf_wwd = 11 ; jf_avt = 12 ; jfld = jf_avt 383 ! 384 slf_d(jf_tem) = sn_tem ; slf_d(jf_sal) = sn_sal ; slf_d(jf_mld) = sn_mld 385 slf_d(jf_emp) = sn_emp ; slf_d(jf_fmf ) = sn_fmf ; slf_d(jf_ice) = sn_ice 386 slf_d(jf_qsr) = sn_qsr ; slf_d(jf_wnd) = sn_wnd ; slf_d(jf_avt) = sn_avt 387 slf_d(jf_uwd) = sn_uwd ; slf_d(jf_vwd) = sn_vwd ; slf_d(jf_wwd) = sn_wwd 388 263 264 jf_uwd = 1 ; jf_vwd = 2 ; jf_wwd = 3 ; jf_emp = 4 ; jf_avt = 5 265 jf_tem = 6 ; jf_sal = 7 ; jf_mld = 8 ; jf_qsr = 9 266 jf_wnd = 10 ; jf_ice = 11 ; jf_fmf = 12 ; jfld = jf_fmf 267 268 ! 269 slf_d(jf_uwd) = sn_uwd ; slf_d(jf_vwd) = sn_vwd ; slf_d(jf_wwd) = sn_wwd 270 slf_d(jf_emp) = sn_emp ; slf_d(jf_avt) = sn_avt 271 slf_d(jf_tem) = sn_tem ; slf_d(jf_sal) = sn_sal ; slf_d(jf_mld) = sn_mld 272 slf_d(jf_qsr) = sn_qsr ; slf_d(jf_wnd) = sn_wnd ; slf_d(jf_ice) = sn_ice 273 slf_d(jf_fmf) = sn_fmf 274 275 276 ! 277 IF( lk_vvl ) THEN 278 jf_div = jfld + 1 ; jf_empb = jfld + 2 ; jfld = jf_empb 279 slf_d(jf_div) = sn_div ; slf_d(jf_empb) = sn_empb 280 ENDIF 281 ! 282 IF( lk_trabbl ) THEN 283 jf_ubl = jfld + 1 ; jf_vbl = jfld + 2 ; jfld = jf_vbl 284 slf_d(jf_ubl) = sn_ubl ; slf_d(jf_vbl) = sn_vbl 285 ENDIF 389 286 ! 390 287 IF( ln_dynrnf ) THEN 391 jf_rnf = jfld + 1 ;jfld = jf_rnf392 slf_d(jf_rnf) = sn_rnf288 jf_rnf = jfld + 1 ; jfld = jf_rnf 289 slf_d(jf_rnf) = sn_rnf 393 290 ELSE 394 rnf (:,:) = 0._wp 395 ENDIF 396 397 ! 398 IF( .NOT.ln_degrad ) THEN ! no degrad option 399 IF( lk_traldf_eiv .AND. ln_dynbbl ) THEN ! eiv & bbl 400 jf_ubl = jfld + 1 ; jf_vbl = jfld + 2 ; jf_eiw = jfld + 3 ; jfld = jf_eiw 401 slf_d(jf_ubl) = sn_ubl ; slf_d(jf_vbl) = sn_vbl ; slf_d(jf_eiw) = sn_eiw 402 ENDIF 403 IF( .NOT.lk_traldf_eiv .AND. ln_dynbbl ) THEN ! no eiv & bbl 404 jf_ubl = jfld + 1 ; jf_vbl = jfld + 2 ; jfld = jf_vbl 405 slf_d(jf_ubl) = sn_ubl ; slf_d(jf_vbl) = sn_vbl 406 ENDIF 407 IF( lk_traldf_eiv .AND. .NOT.ln_dynbbl ) THEN ! eiv & no bbl 408 jf_eiw = jfld + 1 ; jfld = jf_eiw ; slf_d(jf_eiw) = sn_eiw 409 ENDIF 410 ELSE 411 jf_ahu = jfld + 1 ; jf_ahv = jfld + 2 ; jf_ahw = jfld + 3 ; jfld = jf_ahw 412 slf_d(jf_ahu) = sn_ahu ; slf_d(jf_ahv) = sn_ahv ; slf_d(jf_ahw) = sn_ahw 413 IF( lk_traldf_eiv .AND. ln_dynbbl ) THEN ! eiv & bbl 414 jf_ubl = jfld + 1 ; jf_vbl = jfld + 2 ; 415 slf_d(jf_ubl) = sn_ubl ; slf_d(jf_vbl) = sn_vbl 416 jf_eiu = jfld + 3 ; jf_eiv = jfld + 4 ; jf_eiw = jfld + 5 ; jfld = jf_eiw 417 slf_d(jf_eiu) = sn_eiu ; slf_d(jf_eiv) = sn_eiv ; slf_d(jf_eiw) = sn_eiw 418 ENDIF 419 IF( .NOT.lk_traldf_eiv .AND. ln_dynbbl ) THEN ! no eiv & bbl 420 jf_ubl = jfld + 1 ; jf_vbl = jfld + 2 ; jfld = jf_vbl 421 slf_d(jf_ubl) = sn_ubl ; slf_d(jf_vbl) = sn_vbl 422 ENDIF 423 IF( lk_traldf_eiv .AND. .NOT.ln_dynbbl ) THEN ! eiv & no bbl 424 jf_eiu = jfld + 1 ; jf_eiv = jfld + 2 ; jf_eiw = jfld + 3 ; jfld = jf_eiw 425 slf_d(jf_eiu) = sn_eiu ; slf_d(jf_eiv) = sn_eiv ; slf_d(jf_eiw) = sn_eiw 426 ENDIF 427 ENDIF 291 rnf(:,:) = 0._wp 292 ENDIF 293 428 294 429 295 ALLOCATE( sf_dyn(jfld), STAT=ierr ) ! set sf structure 430 IF( ierr > 0 ) THEN296 IF( ierr > 0 ) THEN 431 297 CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' ) ; RETURN 432 298 ENDIF 433 299 ! ! fill sf with slf_i and control print 434 300 CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' ) 301 ! 435 302 ! Open file for each variable to get his number of dimension 436 303 DO ifpr = 1, jfld … … 456 323 ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2), & 457 324 & wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2), STAT=ierr2 ) 458 ELSE 459 ALLOCATE( uslpnow (jpi,jpj,jpk) , vslpnow (jpi,jpj,jpk) , & 460 & wslpinow(jpi,jpj,jpk) , wslpjnow(jpi,jpj,jpk) , STAT=ierr2 ) 461 ENDIF 462 IF( ierr2 > 0 ) THEN 463 CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' ) ; RETURN 325 ! 326 IF( ierr2 > 0 ) THEN 327 CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' ) ; RETURN 328 ENDIF 464 329 ENDIF 465 330 ENDIF 466 IF( ln_dynwzv ) THEN ! slopes 467 IF( sf_dyn(jf_uwd)%ln_tint ) THEN ! time interpolation 468 ALLOCATE( wdta(jpi,jpj,jpk,2), STAT=ierr3 ) 469 ELSE 470 ALLOCATE( wnow(jpi,jpj,jpk) , STAT=ierr3 ) 471 ENDIF 472 IF( ierr3 > 0 ) THEN 473 CALL ctl_stop( 'dta_dyn_init : unable to allocate wdta arrays' ) ; RETURN 474 ENDIF 475 ENDIF 476 ! 477 CALL dta_dyn( nit000 ) 478 ! 479 END SUBROUTINE dta_dyn_init 480 481 SUBROUTINE dta_dyn_wzv( pu, pv, pw ) 482 !!---------------------------------------------------------------------- 483 !! *** ROUTINE wzv *** 484 !! 485 !! ** Purpose : Compute the now vertical velocity after the array swap 486 !! 487 !! ** Method : - compute the now divergence given by : 488 !! * z-coordinate ONLY !!!! 489 !! hdiv = 1/(e1t*e2t) [ di(e2u u) + dj(e1v v) ] 490 !! - Using the incompressibility hypothesis, the vertical 491 !! velocity is computed by integrating the horizontal divergence 492 !! from the bottom to the surface. 493 !! The boundary conditions are w=0 at the bottom (no flux). 494 !!---------------------------------------------------------------------- 495 USE oce, ONLY: zhdiv => hdivn 496 ! 497 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu, pv !: horizontal velocities 498 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: pw !: vertical velocity 499 !! 500 INTEGER :: ji, jj, jk 501 REAL(wp) :: zu, zu1, zv, zv1, zet 502 !!---------------------------------------------------------------------- 503 ! 504 ! Computation of vertical velocity using horizontal divergence 505 zhdiv(:,:,:) = 0._wp 506 DO jk = 1, jpkm1 507 DO jj = 2, jpjm1 508 DO ji = fs_2, fs_jpim1 ! vector opt. 509 zu = pu(ji ,jj ,jk) * umask(ji ,jj ,jk) * e2u(ji ,jj ) * fse3u(ji ,jj ,jk) 510 zu1 = pu(ji-1,jj ,jk) * umask(ji-1,jj ,jk) * e2u(ji-1,jj ) * fse3u(ji-1,jj ,jk) 511 zv = pv(ji ,jj ,jk) * vmask(ji ,jj ,jk) * e1v(ji ,jj ) * fse3v(ji ,jj ,jk) 512 zv1 = pv(ji ,jj-1,jk) * vmask(ji ,jj-1,jk) * e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) 513 zet = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 514 zhdiv(ji,jj,jk) = ( zu - zu1 + zv - zv1 ) * zet 331 ! 332 IF( lk_vvl ) THEN 333 IF( ln_ssh_ini ) THEN ! Restart: read in restart file 334 IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the dynamics restart file for initialisation' 335 CALL iom_open( 'restart', inum ) 336 CALL iom_get( inum, jpdom_autoglo, 'sshn', sshn(:,:) ) 337 CALL iom_get( inum, jpdom_autoglo, 'sshb', sshb(:,:) ) 338 CALL iom_close( inum ) ! close file 339 ELSE 340 IF(lwp) WRITE(numout,*) ' sshn forcing fields read in passive tracers restart file for initialisation' 341 CALL iom_get( numrtr, jpdom_autoglo, 'sshn', sshn(:,:) ) 342 CALL iom_get( numrtr, jpdom_autoglo, 'sshb', sshb(:,:) ) 343 ENDIF 344 ! 345 DO jk = 1, jpkm1 346 fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 347 ENDDO 348 fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 349 350 ! Horizontal scale factor interpolations 351 ! -------------------------------------- 352 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 353 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 354 355 ! Vertical scale factor interpolations 356 ! ------------------------------------ 357 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n(:,:,:), 'W' ) 358 359 fse3t_b(:,:,:) = fse3t_n(:,:,:) 360 fse3u_b(:,:,:) = fse3u_n(:,:,:) 361 fse3v_b(:,:,:) = fse3v_n(:,:,:) 362 363 ! t- and w- points depth 364 ! ---------------------- 365 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 366 fsdepw_n(:,:,1) = 0.0_wp 367 368 DO jk = 2, jpk 369 DO jj = 1,jpj 370 DO ji = 1,jpi 371 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere 372 ! tmask = wmask, ie everywhere expect at jk = mikt 373 ! 1 for jk = 374 ! mikt 375 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 376 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 377 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 378 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 379 END DO 380 END DO 381 END DO 382 383 fsdept_b(:,:,:) = fsdept_n(:,:,:) 384 fsdepw_b(:,:,:) = fsdepw_n(:,:,:) 385 ! 386 ENDIF 387 ! 388 IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN ! read depht over which runoffs are distributed 389 IF(lwp) WRITE(numout,*) 390 IF(lwp) WRITE(numout,*) ' read in the file depht over which runoffs are distributed' 391 CALL iom_open ( "runoffs", inum ) ! open file 392 CALL iom_get ( inum, jpdom_data, 'rodepth', h_rnf ) ! read the river mouth array 393 CALL iom_close( inum ) ! close file 394 ! 395 nk_rnf(:,:) = 0 ! set the number of level over which river runoffs are applied 396 DO jj = 1, jpj 397 DO ji = 1, jpi 398 IF( h_rnf(ji,jj) > 0._wp ) THEN 399 jk = 2 400 DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1 401 END DO 402 nk_rnf(ji,jj) = jk 403 ELSEIF( h_rnf(ji,jj) == -1._wp ) THEN ; nk_rnf(ji,jj) = 1 404 ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN ; nk_rnf(ji,jj) = mbkt(ji,jj) 405 ELSE 406 CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999' ) 407 WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj) 408 ENDIF 515 409 END DO 516 410 END DO 411 DO jj = 1, jpj ! set the associated depth 412 DO ji = 1, jpi 413 h_rnf(ji,jj) = 0._wp 414 DO jk = 1, nk_rnf(ji,jj) 415 h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk) 416 END DO 417 END DO 418 END DO 419 ELSE ! runoffs applied at the surface 420 nk_rnf(:,:) = 1 421 h_rnf (:,:) = fse3t(:,:,1) 422 ENDIF 423 nkrnf_max = MAXVAL( nk_rnf(:,:) ) 424 hrnf_max = MAXVAL( h_rnf(:,:) ) 425 IF( lk_mpp ) THEN 426 CALL mpp_max( nkrnf_max ) ! max over the global domain 427 CALL mpp_max( hrnf_max ) ! max over the global domain 428 ENDIF 429 IF(lwp) WRITE(numout,*) ' ' 430 IF(lwp) WRITE(numout,*) ' max depht of runoff : ', hrnf_max,' max level : ', nkrnf_max 431 IF(lwp) WRITE(numout,*) ' ' 432 ! 433 CALL dta_dyn( nit000 ) 434 ! 435 END SUBROUTINE dta_dyn_init 436 437 SUBROUTINE dta_dyn_swp( kt ) 438 !!--------------------------------------------------------------------- 439 !! *** ROUTINE dta_dyn_swp *** 440 !! 441 !! ** Purpose : Swap and the data and compute the vertical scale factor at U/V/W point 442 !! and the depht 443 !! 444 !!--------------------------------------------------------------------- 445 INTEGER, INTENT(in) :: kt ! time step 446 INTEGER :: ji, jj, jk 447 REAL(wp) :: zcoef 448 ! 449 !!--------------------------------------------------------------------- 450 451 IF( kt == nit000 ) THEN 452 IF(lwp) WRITE(numout,*) 453 IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height' 454 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 455 ENDIF 456 457 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:)) ! before <-- now filtered 458 sshn(:,:) = ssha(:,:) 459 460 fse3t_n(:,:,:) = fse3t_a(:,:,:) 461 462 ! Reconstruction of all vertical scale factors at now and before time steps 463 ! ============================================================================= 464 465 ! Horizontal scale factor interpolations 466 ! -------------------------------------- 467 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 468 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 469 470 ! Vertical scale factor interpolations 471 ! ------------------------------------ 472 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 473 474 fse3t_b(:,:,:) = fse3t_n(:,:,:) 475 fse3u_b(:,:,:) = fse3u_n(:,:,:) 476 fse3v_b(:,:,:) = fse3v_n(:,:,:) 477 478 ! t- and w- points depth 479 ! ---------------------- 480 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 481 fsdepw_n(:,:,1) = 0.0_wp 482 483 DO jk = 2, jpk 484 DO jj = 1,jpj 485 DO ji = 1,jpi 486 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 487 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 488 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 489 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 490 END DO 491 END DO 492 END DO 493 494 fsdept_b(:,:,:) = fsdept_n(:,:,:) 495 fsdepw_b(:,:,:) = fsdepw_n(:,:,:) 496 497 ! 498 END SUBROUTINE dta_dyn_swp 499 500 SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb, pemp, pssha, pe3ta ) 501 !!---------------------------------------------------------------------- 502 !! *** ROUTINE dta_dyn_wzv *** 503 !! 504 !! ** Purpose : compute the after ssh (ssha) and the now vertical velocity 505 !! 506 !! ** Method : Using the incompressibility hypothesis, 507 !! - the ssh increment is computed by integrating the horizontal divergence 508 !! and multiply by the time step. 509 !! 510 !! - compute the after scale factor : repartition of ssh INCREMENT proportionnaly 511 !! to the level thickness ( z-star case ) 512 !! 513 !! - the vertical velocity is computed by integrating the horizontal divergence 514 !! from the bottom to the surface minus the scale factor evolution. 515 !! The boundary conditions are w=0 at the bottom (no flux) 516 !! 517 !! ** action : ssha / e3t_a / wn 518 !! 519 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 520 !!---------------------------------------------------------------------- 521 !! * Arguments 522 INTEGER, INTENT(in ) :: kt ! time-step 523 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: phdivtr ! horizontal divergence transport 524 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: psshb ! now ssh 525 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: pemp ! evaporation minus precipitation 526 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(inout) :: pssha ! after ssh 527 REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out) :: pe3ta ! after vertical scale factor 528 !! * Local declarations 529 INTEGER :: jk 530 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv 531 REAL(wp) :: z2dt 532 !!---------------------------------------------------------------------- 533 534 ! 535 z2dt = 2._wp * rdt 536 ! 537 zhdiv(:,:) = 0._wp 538 DO jk = 1, jpkm1 539 zhdiv(:,:) = zhdiv(:,:) + phdivtr(:,:,jk) * tmask(:,:,jk) 517 540 END DO 518 CALL lbc_lnk( zhdiv, 'T', 1. ) ! Lateral boundary conditions on zhdiv519 !520 ! computation of vertical velocity from the bottom521 pw(:,:,jpk) = 0._wp522 DO jk = jpkm1, 1, -1523 pw(:,:,jk) = pw(:,:,jk+1) - fse3t(:,:,jk) * zhdiv(:,:,jk)541 ! ! Sea surface elevation time-stepping 542 pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rau0 * pemp(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 543 ! ! 544 ! ! After acale factors at t-points ( z_star coordinate ) 545 DO jk = 1, jpkm1 546 pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 524 547 END DO 525 548 ! 526 END SUBROUTINE dta_dyn_wzv 527 528 SUBROUTINE dta_dyn_slp( kt, pts, puslp, pvslp, pwslpi, pwslpj ) 549 END SUBROUTINE dta_dyn_ssh 550 551 552 SUBROUTINE dta_dyn_hrnf 553 !!---------------------------------------------------------------------- 554 !! *** ROUTINE sbc_rnf *** 555 !! 556 !! ** Purpose : update the horizontal divergence with the runoff inflow 557 !! 558 !! ** Method : 559 !! CAUTION : rnf is positive (inflow) decreasing the 560 !! divergence and expressed in m/s 561 !! 562 !! ** Action : phdivn decreased by the runoff inflow 563 !!---------------------------------------------------------------------- 564 !! 565 INTEGER :: ji, jj, jk ! dummy loop indices 566 !!---------------------------------------------------------------------- 567 ! 568 DO jj = 1, jpj ! update the depth over which runoffs are distributed 569 DO ji = 1, jpi 570 h_rnf(ji,jj) = 0._wp 571 DO jk = 1, nk_rnf(ji,jj) ! recalculates h_rnf to be the depth in metres 572 h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk) ! to the bottom of the relevant grid box 573 END DO 574 END DO 575 END DO 576 ! 577 END SUBROUTINE dta_dyn_hrnf 578 579 580 581 SUBROUTINE dta_dyn_slp( kt ) 582 !!--------------------------------------------------------------------- 583 !! *** ROUTINE dta_dyn_slp *** 584 !! 585 !! ** Purpose : Computation of slope 586 !! 587 !!--------------------------------------------------------------------- 588 USE oce, ONLY: zts => tsa 589 ! 590 INTEGER, INTENT(in) :: kt ! time step 591 ! 592 INTEGER :: ji, jj ! dummy loop indices 593 REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation 594 REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation 595 INTEGER :: iswap 596 REAL(wp), POINTER, DIMENSION(:,:,:) :: zuslp, zvslp, zwslpi, zwslpj 597 !!--------------------------------------------------------------------- 598 ! 599 CALL wrk_alloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj ) 600 ! 601 IF( sf_dyn(jf_tem)%ln_tint ) THEN ! Computes slopes (here avt is used as workspace) 602 IF( kt == nit000 ) THEN 603 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:) ! temperature 604 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:) ! salinity 605 avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:) ! vertical diffusive coef. 606 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 607 uslpdta (:,:,:,1) = zuslp (:,:,:) 608 vslpdta (:,:,:,1) = zvslp (:,:,:) 609 wslpidta(:,:,:,1) = zwslpi(:,:,:) 610 wslpjdta(:,:,:,1) = zwslpj(:,:,:) 611 ! 612 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:) ! temperature 613 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity 614 avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef. 615 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 616 uslpdta (:,:,:,2) = zuslp (:,:,:) 617 vslpdta (:,:,:,2) = zvslp (:,:,:) 618 wslpidta(:,:,:,2) = zwslpi(:,:,:) 619 wslpjdta(:,:,:,2) = zwslpj(:,:,:) 620 ELSE 621 ! 622 iswap = 0 623 IF( sf_dyn(jf_tem)%nrec_a(2) - nprevrec /= 0 ) iswap = 1 624 IF( nsecdyn > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap == 1 ) THEN ! read/update the after data 625 IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt 626 uslpdta (:,:,:,1) = uslpdta (:,:,:,2) ! swap the data 627 vslpdta (:,:,:,1) = vslpdta (:,:,:,2) 628 wslpidta(:,:,:,1) = wslpidta(:,:,:,2) 629 wslpjdta(:,:,:,1) = wslpjdta(:,:,:,2) 630 ! 631 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:) ! temperature 632 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity 633 avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef. 634 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 635 ! 636 uslpdta (:,:,:,2) = zuslp (:,:,:) 637 vslpdta (:,:,:,2) = zvslp (:,:,:) 638 wslpidta(:,:,:,2) = zwslpi(:,:,:) 639 wslpjdta(:,:,:,2) = zwslpj(:,:,:) 640 ENDIF 641 ENDIF 642 ENDIF 643 ! 644 IF( sf_dyn(jf_tem)%ln_tint ) THEN 645 ztinta = REAL( nsecdyn - sf_dyn(jf_tem)%nrec_b(2), wp ) & 646 & / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp ) 647 ztintb = 1. - ztinta 648 #if defined key_ldfslp && ! defined key_c1d 649 uslp (:,:,:) = ztintb * uslpdta (:,:,:,1) + ztinta * uslpdta (:,:,:,2) 650 vslp (:,:,:) = ztintb * vslpdta (:,:,:,1) + ztinta * vslpdta (:,:,:,2) 651 wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1) + ztinta * wslpidta(:,:,:,2) 652 wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1) + ztinta * wslpjdta(:,:,:,2) 653 #endif 654 ELSE 655 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature 656 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity 657 avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coef. 658 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 659 ! 660 #if defined key_ldfslp && ! defined key_c1d 661 uslp (:,:,:) = zuslp (:,:,:) 662 vslp (:,:,:) = zvslp (:,:,:) 663 wslpi(:,:,:) = zwslpi(:,:,:) 664 wslpj(:,:,:) = zwslpj(:,:,:) 665 #endif 666 ENDIF 667 ! 668 CALL wrk_dealloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj ) 669 ! 670 END SUBROUTINE dta_dyn_slp 671 672 SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj ) 529 673 !!--------------------------------------------------------------------- 530 674 !! *** ROUTINE dta_dyn_slp *** … … 568 712 #endif 569 713 ! 570 END SUBROUTINE dta_dyn_slp 714 END SUBROUTINE compute_slopes 715 571 716 !!====================================================================== 572 717 END MODULE dtadyn -
branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90
r5602 r7806 34 34 USE trcstp ! passive tracer time-stepping (trc_stp routine) 35 35 USE dtadyn ! Lecture and interpolation of the dynamical fields 36 ! ! I/O & MPP36 ! ! I/O & MPP 37 37 USE iom ! I/O library 38 38 USE in_out_manager ! I/O manager … … 50 50 USE trcnam 51 51 USE trcrst 52 USE diaptr ! Need to initialise this as some variables are used in if statements later53 52 54 53 IMPLICIT NONE … … 94 93 istp = nit000 95 94 ! 96 CALL iom_init( cxios_context ) ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS)97 !98 95 DO WHILE ( istp <= nitend .AND. nstop == 0 ) ! time stepping 99 96 ! 100 IF( istp /= nit000 ) CALL day ( istp ) ! Calendar (day was already called at nit000 in day_init) 101 CALL iom_setkt( istp - nit000 + 1, "nemo" ) ! say to iom that we are at time step kstp 102 CALL dta_dyn ( istp ) ! Interpolation of the dynamical fields 103 CALL trc_stp ( istp ) ! time-stepping 104 CALL stp_ctl ( istp, indic ) ! Time loop: control and print 97 IF( istp == nit000 ) CALL iom_init( cxios_context ) ! iom_put initialization 98 IF( istp /= nit000 ) CALL day ( istp ) ! Calendar (day was already called at nit000 in day_init) 99 CALL iom_setkt ( istp - nit000 + 1, cxios_context ) ! say to iom that we are at time step kstp 100 CALL trc_rst_opn( istp ) ! Open tracer ! restart file 101 CALL dta_dyn ( istp ) ! Interpolation of the dynamical fields 102 CALL trc_stp ( istp ) ! time-stepping 103 IF( lk_vvl ) CALL dta_dyn_swp( istp ) ! swap of sea surface height and vertical scale factors 104 CALL stp_ctl ( istp, indic ) ! Time loop: control and print 105 105 istp = istp + 1 106 106 IF( lk_mpp ) CALL mpp_max( nstop ) … … 265 265 IF( nn_timing == 1 ) CALL timing_start( 'nemo_init') 266 266 ! 267 CALL phy_cst ! Physical constants 268 CALL eos_init ! Equation of state 269 IF( lk_c1d ) CALL c1d_init ! 1D column configuration 270 CALL dom_cfg ! Domain configuration 267 CALL phy_cst ! Physical constants 268 CALL eos_init ! Equation of state 269 IF( lk_c1d ) CALL c1d_init ! 1D column configuration 270 CALL dom_cfg ! Domain configuration 271 ! 271 272 ! 272 273 INQUIRE( FILE='coordinates.nc', EXIST = llexist ) ! Check if coordinate file exist 273 274 ! 274 IF( llexist ) THEN 275 ELSE 275 IF( llexist ) THEN ; CALL dom_init ! compute the grid from coordinates and bathymetry 276 ELSE ; CALL dom_rea ! read grid from the meskmask 276 277 ENDIF 277 278 CALL istate_init ! ocean initial state (Dynamics and tracers) 278 279 279 IF( ln_nnogather ) CALL nemo_northcomms ! Initialise the northfold neighbour lists (must be done after the masks are defined)280 281 IF( ln_ctl ) CALLprt_ctl_init ! Print control280 IF( ln_nnogather ) CALL nemo_northcomms ! Initialise the northfold neighbour lists (must be done after the masks are defined) 281 282 IF( ln_ctl ) CALL prt_ctl_init ! Print control 282 283 283 284 CALL sbc_init ! Forcings : surface module … … 289 290 290 291 CALL tra_qsr_init ! penetrative solar radiation qsr 291 IF( lk_trabbl )CALL tra_bbl_init ! advective (and/or diffusive) bottom boundary layer scheme292 IF( lk_trabbl ) CALL tra_bbl_init ! advective (and/or diffusive) bottom boundary layer scheme 292 293 293 294 CALL trc_nam_run ! Needed to get restart parameters for passive tracers 294 295 CALL trc_rst_cal( nit000, 'READ' ) ! calendar 295 296 CALL dta_dyn_init ! Initialization for the dynamics 296 297 297 CALL trc_init ! Passive tracers initialization 298 CALL dia_ptr_init ! Initialise diaptr as some variables are used299 298 ! ! in various advection and diffusion routines 300 299 IF(lwp) WRITE(numout,cform_aaa) ! Flag AAAAAAA
Note: See TracChangeset
for help on using the changeset viewer.