Changeset 6952
- Timestamp:
- 2016-09-23T16:02:11+02:00 (8 years ago)
- Location:
- branches/2015/dev_r6946_Offline_vvl/NEMOGCM
- Files:
-
- 13 added
- 11 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/OFF_SRC/domrea.F90
r5504 r6952 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_r6946_Offline_vvl/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90
r5767 r6952 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_dynrnf !: read runoff data in file (T) or set to zero (F) 56 LOGICAL :: ln_dynrnf_depth !: read runoff data in file (T) or set to zero (F) 57 REAL(wp) :: fwbcorr 58 59 60 INTEGER , PARAMETER :: jpfld = 20 ! maximum number of fields to read 56 61 INTEGER , SAVE :: jf_tem ! index of temperature 57 62 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-wind63 INTEGER , SAVE :: jf_uwd ! index of u-transport 64 INTEGER , SAVE :: jf_vwd ! index of v-transport 65 INTEGER , SAVE :: jf_wwd ! index of v-transport 61 66 INTEGER , SAVE :: jf_avt ! index of Kz 62 67 INTEGER , SAVE :: jf_mld ! index of mixed layer deptht 63 INTEGER , SAVE :: jf_emp ! index of water flux 68 INTEGER , SAVE :: jf_empn ! index of water flux 69 INTEGER , SAVE :: jf_empb ! index of water flux 64 70 INTEGER , SAVE :: jf_qsr ! index of solar radiation 65 71 INTEGER , SAVE :: jf_wnd ! index of wind speed 66 72 INTEGER , SAVE :: jf_ice ! index of sea ice cover 67 73 INTEGER , SAVE :: jf_rnf ! index of river runoff 74 INTEGER , SAVE :: jf_fmf ! index of downward salt flux 68 75 INTEGER , SAVE :: jf_ubl ! index of u-bbl coef 69 76 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) 77 INTEGER , SAVE :: jf_div ! index of e3t 78 79 80 TYPE(FLD), ALLOCATABLE, SAVE, DIMENSION(:) :: sf_dyn ! structure of input fields (file informations, fields read) 79 81 ! ! 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 82 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uslpdta ! zonal isopycnal slopes 83 83 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vslpdta ! meridional isopycnal slopes 84 84 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpidta ! zonal diapycnal slopes 85 85 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 86 92 87 93 88 !! * Substitutions … … 113 108 !!---------------------------------------------------------------------- 114 109 ! 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 ! 110 USE oce, ONLY: zhdivtr => ua 120 111 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 ! 112 INTEGER :: ji, jj, jk 113 REAL(wp), POINTER, DIMENSION(:,:) :: zemp 114 ! 127 115 !!---------------------------------------------------------------------- 128 116 … … 130 118 IF( nn_timing == 1 ) CALL timing_start( 'dta_dyn') 131 119 ! 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 120 CALL fld_read( kt, 1, sf_dyn ) != read data at kt time step ==! 243 121 ! 244 122 tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature 245 123 tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity 246 ! 124 wndm(:,:) = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1) ! wind speed - needed for gas exchange 125 fmmflx(:,:) = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1) ! downward salt flux (v3.5+) 126 fr_i(:,:) = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1) ! Sea-ice fraction 127 qsr (:,:) = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1) ! solar radiation 128 emp (:,:) = sf_dyn(jf_empn)%fnow(:,:,1) * tmask(:,:,1) ! E-P 129 emp_b (:,:) = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1) ! E-P 130 IF( ln_dynrnf ) THEN 131 rnf (:,:) = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1) ! E-P 132 IF( ln_dynrnf_depth .AND. lk_vvl ) CALL dta_dyn_hrnf 133 ENDIF 134 ! 135 un(:,:,:) = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:) ! effective u-transport 136 vn(:,:,:) = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:) ! effective v-transport 137 wn(:,:,:) = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:) ! effective v-transport 138 zhdivtr(:,:,:) = sf_dyn(jf_div)%fnow(:,:,:) * tmask(:,:,:) ! effective u-transport 139 ! 140 IF( lk_vvl ) THEN 141 CALL wrk_alloc(jpi, jpj, zemp ) 142 zemp(:,:) = 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr * tmask(:,:,1) 143 CALL dta_dyn_ssh( kt, zhdivtr, sshb, zemp, ssha, fse3t_a(:,:,:) ) != ssh, vertical scale factor & vertical transport 144 CALL wrk_dealloc(jpi, jpj, zemp ) 145 ! Write in the tracer restart file 146 ! ******************************* 147 IF( lrst_trc ) THEN 148 IF(lwp) WRITE(numout,*) 149 IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file ', & 150 & 'at it= ', kt,' date= ', ndastp 151 IF(lwp) WRITE(numout,*) '~~~~' 152 CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssha ) 153 CALL iom_rstput( kt, nitrst, numrtw, 'sshb', sshn ) 154 ENDIF 155 ENDIF 247 156 ! 248 157 CALL eos ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop … … 251 160 252 161 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 162 CALL zdf_mxl( kt ) ! In any case, we need mxl 163 ! 164 IF( lk_ldfslp .AND. .NOT.lk_c1d ) CALL dta_dyn_slp( kt ) ! Computation of slopes 165 ! 166 hmld(:,:) = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1) ! mixed layer depht 167 avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coefficient 168 ! 270 169 #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 170 ahu_bbl(:,:) = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1) ! bbl diffusive coef 171 ahv_bbl(:,:) = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1) 278 172 #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 173 ! 174 ! 175 CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 302 176 ! 303 177 IF(ln_ctl) THEN ! print control … … 308 182 CALL prt_ctl(tab3d_1=wn , clinfo1=' wn - : ', mask1=tmask, ovlap=1, kdim=jpk ) 309 183 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 )184 ! CALL prt_ctl(tab2d_1=fr_i , clinfo1=' fr_i - : ', mask1=tmask, ovlap=1 ) 185 ! CALL prt_ctl(tab2d_1=hmld , clinfo1=' hmld - : ', mask1=tmask, ovlap=1 ) 186 ! CALL prt_ctl(tab2d_1=fmmflx , clinfo1=' fmmflx - : ', mask1=tmask, ovlap=1 ) 187 ! CALL prt_ctl(tab2d_1=emp , clinfo1=' emp - : ', mask1=tmask, ovlap=1 ) 188 ! CALL prt_ctl(tab2d_1=wndm , clinfo1=' wspd - : ', mask1=tmask, ovlap=1 ) 189 ! CALL prt_ctl(tab2d_1=qsr , clinfo1=' qsr - : ', mask1=tmask, ovlap=1 ) 316 190 ENDIF 317 191 ! … … 335 209 INTEGER :: inum, idv, idimv ! local integer 336 210 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 211 INTEGER :: ji, jj, jk 212 REAL(wp) :: zcoef 213 INTEGER :: nkrnf_max 214 REAL(wp) :: hrnf_max 215 !! 216 CHARACTER(len=100) :: cn_dir ! Root directory for location of core files 217 TYPE(FLD_N), DIMENSION(jpfld) :: slf_d ! array of namelist informations on the fields to read 218 TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_empb, sn_empn ! informations about the fields to be read 219 TYPE(FLD_N) :: sn_tem , sn_sal , sn_avt ! " " 220 TYPE(FLD_N) :: sn_mld, sn_qsr, sn_wnd , sn_ice , sn_fmf ! " " 221 TYPE(FLD_N) :: sn_ubl, sn_vbl, sn_rnf ! " " 222 TYPE(FLD_N) :: sn_div ! informations about the fields to be read 223 224 !!---------------------------------------------------------------------- 225 ! 226 NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth, fwbcorr, & 227 & sn_uwd, sn_vwd, sn_wwd, sn_empb, sn_empn, & 228 & sn_avt, sn_tem, sn_sal, sn_mld , sn_qsr , & 229 & sn_wnd, sn_ice, sn_fmf, & 230 & sn_ubl, sn_vbl, sn_rnf, & 231 & sn_div 349 232 ! 350 233 REWIND( numnam_ref ) ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data … … 363 246 WRITE(numout,*) '~~~~~~~ ' 364 247 WRITE(numout,*) ' Namelist namdta_dyn' 365 WRITE(numout,*) ' vertical velocity read from file (T) or computed (F) ln_dynwzv = ', ln_dynwzv 366 WRITE(numout,*) ' bbl coef read from file (T) or computed (F) ln_dynbbl = ', ln_dynbbl 367 WRITE(numout,*) ' degradation option enabled (T) or not (F) ln_degrad = ', ln_degrad 368 WRITE(numout,*) ' river runoff option enabled (T) or not (F) ln_dynrnf = ', ln_dynrnf 248 WRITE(numout,*) ' runoffs option enabled (T) or not (F) ln_dynrnf = ', ln_dynrnf 249 WRITE(numout,*) ' runoffs is spread in vertical ln_dynrnf_depth = ', ln_dynrnf_depth 250 WRITE(numout,*) ' annual global mean of empmr for ssh correction fwbcorr = ', fwbcorr 369 251 WRITE(numout,*) 370 252 ENDIF 371 253 ! 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 254 255 jf_uwd = 1 ; jf_vwd = 2 ; jf_wwd = 3 256 jf_empb = 4 ; jf_empn = 5 ; jf_avt = 6 257 jf_tem = 7 ; jf_sal = 8 ; jf_mld = 9 ; jf_qsr = 10 258 jf_wnd = 11 ; jf_ice = 12 ; jf_fmf = 13 ; jfld = jf_fmf 259 260 ! 261 slf_d(jf_uwd) = sn_uwd ; slf_d(jf_vwd) = sn_vwd ; slf_d(jf_wwd) = sn_wwd 262 slf_d(jf_empb) = sn_empb ; slf_d(jf_empn) = sn_empn ; slf_d(jf_avt) = sn_avt 263 slf_d(jf_tem) = sn_tem ; slf_d(jf_sal) = sn_sal ; slf_d(jf_mld) = sn_mld 264 slf_d(jf_qsr) = sn_qsr ; slf_d(jf_wnd) = sn_wnd ; slf_d(jf_ice) = sn_ice 265 slf_d(jf_fmf) = sn_fmf 266 267 ! 268 IF( lk_vvl ) THEN 269 jf_div = jfld + 1 ; jfld = jf_div 270 slf_d(jf_div) = sn_div 271 ENDIF 272 ! 273 IF( lk_trabbl ) THEN 274 jf_ubl = jfld + 1 ; jf_vbl = jfld + 2 ; jfld = jf_vbl 275 slf_d(jf_ubl) = sn_ubl ; slf_d(jf_vbl) = sn_vbl 276 ENDIF 389 277 ! 390 278 IF( ln_dynrnf ) THEN 391 jf_rnf = jfld + 1 ;jfld = jf_rnf392 slf_d(jf_rnf) = sn_rnf279 jf_rnf = jfld + 1 ; jfld = jf_rnf 280 slf_d(jf_rnf) = sn_rnf 393 281 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 282 rnf(:,:) = 0._wp 283 ENDIF 284 428 285 429 286 ALLOCATE( sf_dyn(jfld), STAT=ierr ) ! set sf structure 430 IF( ierr > 0 ) THEN287 IF( ierr > 0 ) THEN 431 288 CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' ) ; RETURN 432 289 ENDIF 433 290 ! ! fill sf with slf_i and control print 434 291 CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' ) 292 ! 435 293 ! Open file for each variable to get his number of dimension 436 294 DO ifpr = 1, jfld … … 456 314 ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2), & 457 315 & 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 316 ! 317 IF( ierr2 > 0 ) THEN 318 CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' ) ; RETURN 319 ENDIF 464 320 ENDIF 465 321 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 322 ! 323 IF( lk_vvl ) THEN 324 IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND. & ! Restart: read in restart file 325 iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN 326 IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation' 327 CALL iom_get( numrtr, jpdom_autoglo, 'sshn', sshn(:,:) ) 328 CALL iom_get( numrtr, jpdom_autoglo, 'sshb', sshb(:,:) ) 329 ELSE 330 IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation' 331 CALL iom_open( 'restart', inum ) 332 CALL iom_get( inum, jpdom_autoglo, 'sshn', sshn(:,:) ) 333 CALL iom_get( inum, jpdom_autoglo, 'sshb', sshb(:,:) ) 334 CALL iom_close( inum ) ! close file 335 ENDIF 336 ! 337 DO jk = 1, jpkm1 338 fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 339 ENDDO 340 fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 341 342 ! Horizontal scale factor interpolations 343 ! -------------------------------------- 344 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 345 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 346 347 ! Vertical scale factor interpolations 348 ! ------------------------------------ 349 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n(:,:,:), 'W' ) 350 351 fse3t_b(:,:,:) = fse3t_n(:,:,:) 352 fse3u_b(:,:,:) = fse3u_n(:,:,:) 353 fse3v_b(:,:,:) = fse3v_n(:,:,:) 354 355 ! t- and w- points depth 356 ! ---------------------- 357 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 358 fsdepw_n(:,:,1) = 0.0_wp 359 360 DO jk = 2, jpk 361 DO jj = 1,jpj 362 DO ji = 1,jpi 363 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere 364 ! tmask = wmask, ie everywhere expect at jk = mikt 365 ! 1 for jk = 366 ! mikt 367 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 368 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 369 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 370 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 371 END DO 372 END DO 373 END DO 374 375 fsdept_b(:,:,:) = fsdept_n(:,:,:) 376 fsdepw_b(:,:,:) = fsdepw_n(:,:,:) 377 ! 378 ENDIF 379 ! 380 IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN ! read depht over which runoffs are distributed 381 IF(lwp) WRITE(numout,*) 382 IF(lwp) WRITE(numout,*) ' read in the file depht over which runoffs are distributed' 383 CALL iom_open ( "runoffs", inum ) ! open file 384 CALL iom_get ( inum, jpdom_data, 'rodepth', h_rnf ) ! read the river mouth array 385 CALL iom_close( inum ) ! close file 386 ! 387 nk_rnf(:,:) = 0 ! set the number of level over which river runoffs are applied 388 DO jj = 1, jpj 389 DO ji = 1, jpi 390 IF( h_rnf(ji,jj) > 0._wp ) THEN 391 jk = 2 392 DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1 393 END DO 394 nk_rnf(ji,jj) = jk 395 ELSEIF( h_rnf(ji,jj) == -1._wp ) THEN ; nk_rnf(ji,jj) = 1 396 ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN ; nk_rnf(ji,jj) = mbkt(ji,jj) 397 ELSE 398 CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999' ) 399 WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj) 400 ENDIF 515 401 END DO 516 402 END DO 403 DO jj = 1, jpj ! set the associated depth 404 DO ji = 1, jpi 405 h_rnf(ji,jj) = 0._wp 406 DO jk = 1, nk_rnf(ji,jj) 407 h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk) 408 END DO 409 END DO 410 END DO 411 ELSE ! runoffs applied at the surface 412 nk_rnf(:,:) = 1 413 h_rnf (:,:) = fse3t(:,:,1) 414 ENDIF 415 nkrnf_max = MAXVAL( nk_rnf(:,:) ) 416 hrnf_max = MAXVAL( h_rnf(:,:) ) 417 IF( lk_mpp ) THEN 418 CALL mpp_max( nkrnf_max ) ! max over the global domain 419 CALL mpp_max( hrnf_max ) ! max over the global domain 420 ENDIF 421 IF(lwp) WRITE(numout,*) ' ' 422 IF(lwp) WRITE(numout,*) ' max depht of runoff : ', hrnf_max,' max level : ', nkrnf_max 423 IF(lwp) WRITE(numout,*) ' ' 424 ! 425 CALL dta_dyn( nit000 ) 426 ! 427 END SUBROUTINE dta_dyn_init 428 429 SUBROUTINE dta_dyn_swp( kt ) 430 !!--------------------------------------------------------------------- 431 !! *** ROUTINE dta_dyn_swp *** 432 !! 433 !! ** Purpose : Swap and the data and compute the vertical scale factor at U/V/W point 434 !! and the depht 435 !! 436 !!--------------------------------------------------------------------- 437 INTEGER, INTENT(in) :: kt ! time step 438 INTEGER :: ji, jj, jk 439 REAL(wp) :: zcoef 440 ! 441 !!--------------------------------------------------------------------- 442 443 IF( kt == nit000 ) THEN 444 IF(lwp) WRITE(numout,*) 445 IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height' 446 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 447 ENDIF 448 449 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:)) ! before <-- now filtered 450 sshn(:,:) = ssha(:,:) 451 452 fse3t_n(:,:,:) = fse3t_a(:,:,:) 453 454 ! Reconstruction of all vertical scale factors at now and before time steps 455 ! ============================================================================= 456 457 ! Horizontal scale factor interpolations 458 ! -------------------------------------- 459 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 460 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 461 462 ! Vertical scale factor interpolations 463 ! ------------------------------------ 464 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 465 466 fse3t_b(:,:,:) = fse3t_n(:,:,:) 467 fse3u_b(:,:,:) = fse3u_n(:,:,:) 468 fse3v_b(:,:,:) = fse3v_n(:,:,:) 469 470 ! t- and w- points depth 471 ! ---------------------- 472 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 473 fsdepw_n(:,:,1) = 0.0_wp 474 475 DO jk = 2, jpk 476 DO jj = 1,jpj 477 DO ji = 1,jpi 478 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 479 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 480 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 481 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 482 END DO 483 END DO 484 END DO 485 486 fsdept_b(:,:,:) = fsdept_n(:,:,:) 487 fsdepw_b(:,:,:) = fsdepw_n(:,:,:) 488 489 ! 490 END SUBROUTINE dta_dyn_swp 491 492 SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb, pemp, pssha, pe3ta ) 493 !!---------------------------------------------------------------------- 494 !! *** ROUTINE dta_dyn_wzv *** 495 !! 496 !! ** Purpose : compute the after ssh (ssha) and the now vertical velocity 497 !! 498 !! ** Method : Using the incompressibility hypothesis, 499 !! - the ssh increment is computed by integrating the horizontal divergence 500 !! and multiply by the time step. 501 !! 502 !! - compute the after scale factor : repartition of ssh INCREMENT proportionnaly 503 !! to the level thickness ( z-star case ) 504 !! 505 !! - the vertical velocity is computed by integrating the horizontal divergence 506 !! from the bottom to the surface minus the scale factor evolution. 507 !! The boundary conditions are w=0 at the bottom (no flux) 508 !! 509 !! ** action : ssha / e3t_a / wn 510 !! 511 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 512 !!---------------------------------------------------------------------- 513 !! * Arguments 514 INTEGER, INTENT(in ) :: kt ! time-step 515 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: phdivtr ! horizontal divergence transport 516 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: psshb ! now ssh 517 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: pemp ! evaporation minus precipitation 518 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(inout) :: pssha ! after ssh 519 REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out) :: pe3ta ! after vertical scale factor 520 !! * Local declarations 521 INTEGER :: jk 522 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv 523 REAL(wp) :: z2dt 524 !!---------------------------------------------------------------------- 525 526 ! 527 z2dt = 2._wp * rdt 528 ! 529 zhdiv(:,:) = 0._wp 530 DO jk = 1, jpkm1 531 zhdiv(:,:) = zhdiv(:,:) + phdivtr(:,:,jk) * tmask(:,:,jk) 517 532 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)533 ! ! Sea surface elevation time-stepping 534 pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rau0 * pemp(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 535 ! ! 536 ! ! After acale factors at t-points ( z_star coordinate ) 537 DO jk = 1, jpkm1 538 pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 524 539 END DO 525 540 ! 526 END SUBROUTINE dta_dyn_wzv 527 528 SUBROUTINE dta_dyn_slp( kt, pts, puslp, pvslp, pwslpi, pwslpj ) 541 END SUBROUTINE dta_dyn_ssh 542 543 544 SUBROUTINE dta_dyn_hrnf 545 !!---------------------------------------------------------------------- 546 !! *** ROUTINE sbc_rnf *** 547 !! 548 !! ** Purpose : update the horizontal divergence with the runoff inflow 549 !! 550 !! ** Method : 551 !! CAUTION : rnf is positive (inflow) decreasing the 552 !! divergence and expressed in m/s 553 !! 554 !! ** Action : phdivn decreased by the runoff inflow 555 !!---------------------------------------------------------------------- 556 !! 557 INTEGER :: ji, jj, jk ! dummy loop indices 558 !!---------------------------------------------------------------------- 559 ! 560 DO jj = 1, jpj ! update the depth over which runoffs are distributed 561 DO ji = 1, jpi 562 h_rnf(ji,jj) = 0._wp 563 DO jk = 1, nk_rnf(ji,jj) ! recalculates h_rnf to be the depth in metres 564 h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk) ! to the bottom of the relevant grid box 565 END DO 566 END DO 567 END DO 568 ! 569 END SUBROUTINE dta_dyn_hrnf 570 571 572 573 SUBROUTINE dta_dyn_slp( kt ) 574 !!--------------------------------------------------------------------- 575 !! *** ROUTINE dta_dyn_slp *** 576 !! 577 !! ** Purpose : Computation of slope 578 !! 579 !!--------------------------------------------------------------------- 580 USE oce, ONLY: zts => tsa 581 ! 582 INTEGER, INTENT(in) :: kt ! time step 583 ! 584 INTEGER :: ji, jj ! dummy loop indices 585 REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation 586 REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation 587 INTEGER :: iswap, isecdyn, iprevrec 588 REAL(wp), POINTER, DIMENSION(:,:,:) :: zuslp, zvslp, zwslpi, zwslpj 589 !!--------------------------------------------------------------------- 590 ! 591 CALL wrk_alloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj ) 592 ! 593 isecdyn = nsec_year + nsec1jan000 ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step 594 ! 595 IF( kt == nit000 ) THEN ; iprevrec = 0 596 ELSE ; iprevrec = sf_dyn(jf_tem)%nrec_a(2) 597 ENDIF 598 ! 599 IF( sf_dyn(jf_tem)%ln_tint ) THEN ! Computes slopes (here avt is used as workspace) 600 IF( kt == nit000 ) THEN 601 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:) ! temperature 602 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:) ! salinity 603 avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:) ! vertical diffusive coef. 604 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 605 uslpdta (:,:,:,1) = zuslp (:,:,:) 606 vslpdta (:,:,:,1) = zvslp (:,:,:) 607 wslpidta(:,:,:,1) = zwslpi(:,:,:) 608 wslpjdta(:,:,:,1) = zwslpj(:,:,:) 609 ! 610 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:) ! temperature 611 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity 612 avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef. 613 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 614 uslpdta (:,:,:,2) = zuslp (:,:,:) 615 vslpdta (:,:,:,2) = zvslp (:,:,:) 616 wslpidta(:,:,:,2) = zwslpi(:,:,:) 617 wslpjdta(:,:,:,2) = zwslpj(:,:,:) 618 ELSE 619 ! 620 iswap = 0 621 IF( sf_dyn(jf_tem)%nrec_a(2) - iprevrec /= 0 ) iswap = 1 622 IF( isecdyn > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap == 1 ) THEN ! read/update the after data 623 IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt 624 uslpdta (:,:,:,1) = uslpdta (:,:,:,2) ! swap the data 625 vslpdta (:,:,:,1) = vslpdta (:,:,:,2) 626 wslpidta(:,:,:,1) = wslpidta(:,:,:,2) 627 wslpjdta(:,:,:,1) = wslpjdta(:,:,:,2) 628 ! 629 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:) ! temperature 630 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity 631 avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef. 632 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 633 ! 634 uslpdta (:,:,:,2) = zuslp (:,:,:) 635 vslpdta (:,:,:,2) = zvslp (:,:,:) 636 wslpidta(:,:,:,2) = zwslpi(:,:,:) 637 wslpjdta(:,:,:,2) = zwslpj(:,:,:) 638 ENDIF 639 ENDIF 640 ENDIF 641 ! 642 IF( sf_dyn(jf_tem)%ln_tint ) THEN 643 ztinta = REAL( isecdyn - sf_dyn(jf_tem)%nrec_b(2), wp ) & 644 & / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp ) 645 ztintb = 1. - ztinta 646 #if defined key_ldfslp && ! defined key_c1d 647 uslp (:,:,:) = ztintb * uslpdta (:,:,:,1) + ztinta * uslpdta (:,:,:,2) 648 vslp (:,:,:) = ztintb * vslpdta (:,:,:,1) + ztinta * vslpdta (:,:,:,2) 649 wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1) + ztinta * wslpidta(:,:,:,2) 650 wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1) + ztinta * wslpjdta(:,:,:,2) 651 #endif 652 ELSE 653 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature 654 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity 655 avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coef. 656 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 657 ! 658 #if defined key_ldfslp && ! defined key_c1d 659 uslp (:,:,:) = zuslp (:,:,:) 660 vslp (:,:,:) = zvslp (:,:,:) 661 wslpi(:,:,:) = zwslpi(:,:,:) 662 wslpj(:,:,:) = zwslpj(:,:,:) 663 #endif 664 ENDIF 665 ! 666 CALL wrk_dealloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj ) 667 ! 668 END SUBROUTINE dta_dyn_slp 669 670 SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj ) 529 671 !!--------------------------------------------------------------------- 530 672 !! *** ROUTINE dta_dyn_slp *** … … 568 710 #endif 569 711 ! 570 END SUBROUTINE dta_dyn_slp 712 END SUBROUTINE compute_slopes 713 571 714 !!====================================================================== 572 715 END MODULE dtadyn -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90
r5504 r6952 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 -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r5363 r6952 23 23 USE dom_oce ! domain: ocean 24 24 USE sbc_oce ! surface boundary condition: ocean 25 USE trc_oce, ONLY : lk_offline ! 25 26 USE phycst ! physical constants 26 27 USE closea ! closed seas … … 97 98 END DO 98 99 ! 99 IF( lk_vvl ) CALL dom_vvl_init ! Vertical variable mesh 100 ! 101 IF( lk_c1d ) CALL cor_c1d ! 1D configuration: Coriolis set at T-point 102 ! 103 ! 104 hu(:,:) = 0._wp ! Ocean depth at U-points 105 hv(:,:) = 0._wp ! Ocean depth at V-points 106 ht(:,:) = 0._wp ! Ocean depth at T-points 107 DO jk = 1, jpkm1 108 hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 109 hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 110 ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 111 END DO 112 ! ! Inverse of the local depth 113 hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 114 hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 115 100 IF( .NOT. lk_offline ) THEN 101 ! 102 IF( lk_vvl ) CALL dom_vvl_init ! Vertical variable mesh 103 ! 104 IF( lk_c1d ) CALL cor_c1d ! 1D configuration: Coriolis set at T-point 105 ! 106 ! 107 hu(:,:) = 0._wp ! Ocean depth at U-points 108 hv(:,:) = 0._wp ! Ocean depth at V-points 109 ht(:,:) = 0._wp ! Ocean depth at T-points 110 DO jk = 1, jpkm1 111 hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 112 hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 113 ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 114 END DO 115 ! ! Inverse of the local depth 116 hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 117 hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 118 ! 119 ENDIF 116 120 CALL dom_stp ! time step 117 121 IF( nmsh /= 0 ) CALL dom_wri ! Create a domain file -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90
r6943 r6952 261 261 IF( ln_dust .AND. ln_solub ) THEN ; ll_solub = .TRUE. 262 262 ELSE ; ll_solub = .FALSE. 263 ENDIF264 265 ! set the number of level over which river runoffs are applied266 ! online configuration : computed in sbcrnf267 IF( lk_offline ) THEN268 nk_rnf(:,:) = 1269 h_rnf (:,:) = fsdept(:,:,1)270 263 ENDIF 271 264 -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90
r5385 r6952 88 88 r2dt(:) = 2. * rdttrc(:) ! = 2 rdttrc (leapfrog) 89 89 ENDIF 90 ! ! effective transport 91 DO jk = 1, jpkm1 92 ! ! eulerian transport only 93 zun(:,:,jk) = e2u (:,:) * fse3u(:,:,jk) * un(:,:,jk) 94 zvn(:,:,jk) = e1v (:,:) * fse3v(:,:,jk) * vn(:,:,jk) 95 zwn(:,:,jk) = e1e2t(:,:) * wn(:,:,jk) 96 ! 97 END DO 98 ! 99 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 100 zun(:,:,:) = zun(:,:,:) + un_td(:,:,:) 101 zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:) 90 ! 91 IF( lk_offline ) THEN 92 zun(:,:,:) = un(:,:,:) ! effective transport already in un/vn/wn 93 zvn(:,:,:) = vn(:,:,:) 94 zwn(:,:,:) = wn(:,:,:) 95 ELSE 96 ! ! effective transport 97 DO jk = 1, jpkm1 98 ! ! eulerian transport only 99 zun(:,:,jk) = e2u (:,:) * fse3u(:,:,jk) * un(:,:,jk) 100 zvn(:,:,jk) = e1v (:,:) * fse3v(:,:,jk) * vn(:,:,jk) 101 zwn(:,:,jk) = e1e2t(:,:) * wn(:,:,jk) 102 ! 103 END DO 104 ! 105 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 106 zun(:,:,:) = zun(:,:,:) + un_td(:,:,:) 107 zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:) 108 ENDIF 109 ! 110 zun(:,:,jpk) = 0._wp ! no transport trough the bottom 111 zvn(:,:,jpk) = 0._wp ! no transport trough the bottom 112 zwn(:,:,jpk) = 0._wp ! no transport trough the bottom 113 114 IF( lk_traldf_eiv .AND. .NOT. ln_traldf_grif ) & ! add the eiv transport (if necessary) 115 & CALL tra_adv_eiv( kt, nittrc000, zun, zvn, zwn, 'TRC' ) 116 ! 117 IF( ln_mle ) CALL tra_adv_mle( kt, nittrc000, zun, zvn, zwn, 'TRC' ) ! add the mle transport (if necessary) 118 ! 102 119 ENDIF 103 !104 zun(:,:,jpk) = 0._wp ! no transport trough the bottom105 zvn(:,:,jpk) = 0._wp ! no transport trough the bottom106 zwn(:,:,jpk) = 0._wp ! no transport trough the bottom107 108 IF( lk_traldf_eiv .AND. .NOT. ln_traldf_grif ) & ! add the eiv transport (if necessary)109 & CALL tra_adv_eiv( kt, nittrc000, zun, zvn, zwn, 'TRC' )110 !111 IF( ln_mle ) CALL tra_adv_mle( kt, nittrc000, zun, zvn, zwn, 'TRC' ) ! add the mle transport (if necessary)112 120 ! 113 121 SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==! -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/TRP/trcnxt.F90
r6204 r6952 45 45 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: r2dt 46 46 47 !! * Substitutions 48 # include "domzgr_substitute.h90" 49 # include "vectopt_loop_substitute.h90" 47 50 !!---------------------------------------------------------------------- 48 51 !! NEMO/TOP 3.3 , NEMO Consortium (2010) … … 50 53 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 51 54 !!---------------------------------------------------------------------- 55 52 56 CONTAINS 53 57 … … 136 140 ! 137 141 ELSE 138 ! Leap-Frog + Asselin filter time stepping 139 IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt, nittrc000, rdttrc, 'TRC', trb, trn, tra, & 140 & sbc_trc, sbc_trc_b, jptra ) ! variable volume level (vvl) 141 ELSE ; CALL tra_nxt_fix( kt, nittrc000, 'TRC', trb, trn, tra, jptra ) ! fixed volume level 142 IF( .NOT. lk_offline ) THEN ! Leap-Frog + Asselin filter time stepping 143 IF( lk_vvl ) THEN ; CALL tra_nxt_vvl( kt, nittrc000, rdttrc, 'TRC', trb, trn, tra, & 144 & sbc_trc, sbc_trc_b, jptra ) ! variable volume level (vvl) 145 ELSE ; CALL tra_nxt_fix( kt, nittrc000, 'TRC', trb, trn, tra, jptra ) ! fixed volume level 146 ENDIF 147 ELSE 148 CALL trc_nxt_off( kt ) ! offline 142 149 ENDIF 143 150 ENDIF … … 165 172 END SUBROUTINE trc_nxt 166 173 174 SUBROUTINE trc_nxt_off( kt ) 175 !!---------------------------------------------------------------------- 176 !! *** ROUTINE tra_nxt_vvl *** 177 !! 178 !! ** Purpose : Time varying volume: apply the Asselin time filter 179 !! and swap the tracer fields. 180 !! 181 !! ** Method : - Apply a thickness weighted Asselin time filter on now fields. 182 !! - save in (ta,sa) a thickness weighted average over the three 183 !! time levels which will be used to compute rdn and thus the semi- 184 !! implicit hydrostatic pressure gradient (ln_dynhpg_imp = T) 185 !! - swap tracer fields to prepare the next time_step. 186 !! This can be summurized for tempearture as: 187 !! ztm = ( e3t_n*tn + rbcp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) ln_dynhpg_imp = T 188 !! /( e3t_n + rbcp*[ e3t_b - 2 e3t_n + e3t_a ] ) 189 !! ztm = 0 otherwise 190 !! tb = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] ) 191 !! /( e3t_n + atfp*[ e3t_b - 2 e3t_n + e3t_a ] ) 192 !! tn = ta 193 !! ta = zt (NB: reset to 0 after eos_bn2 call) 194 !! 195 !! ** Action : - (tb,sb) and (tn,sn) ready for the next time step 196 !! - (ta,sa) time averaged (t,s) (ln_dynhpg_imp = T) 197 !!---------------------------------------------------------------------- 198 INTEGER , INTENT(in ) :: kt ! ocean time-step index 199 !! 200 INTEGER :: ji, jj, jk, jn ! dummy loop indices 201 REAL(wp) :: zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d ! local scalar 202 REAL(wp) :: zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d ! - - 203 !!---------------------------------------------------------------------- 204 ! 205 IF( kt == nittrc000 ) THEN 206 IF(lwp) WRITE(numout,*) 207 IF(lwp) WRITE(numout,*) 'trc_nxt_off : time stepping' 208 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 209 ENDIF 210 ! 211 IF(lwp) WRITE(numout,*) 'trc_nxt_off : time stepping kt = ', kt 212 IF(lwp) WRITE(numout,*) 213 214 DO jn = 1, jptra 215 DO jk = 1, jpkm1 216 zfact1 = atfp * rdttrc(jk) 217 zfact2 = zfact1 / rau0 218 DO jj = 1, jpj 219 DO ji = 1, jpi 220 ze3t_b = fse3t_b(ji,jj,jk) 221 ze3t_n = fse3t_n(ji,jj,jk) 222 ze3t_a = fse3t_a(ji,jj,jk) 223 ! ! tracer content at Before, now and after 224 ztc_b = trb(ji,jj,jk,jn) * ze3t_b 225 ztc_n = trn(ji,jj,jk,jn) * ze3t_n 226 ztc_a = tra(ji,jj,jk,jn) * ze3t_a 227 ! 228 ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b 229 ztc_d = ztc_a - 2. * ztc_n + ztc_b 230 ! 231 ze3t_f = ze3t_n + atfp * ze3t_d 232 ztc_f = ztc_n + atfp * ztc_d 233 ! 234 IF( lk_vvl .AND. jk == mikt(ji,jj) ) THEN ! first level 235 ze3t_f = ze3t_f - zfact2 * ( emp_b(ji,jj) - emp(ji,jj) ) 236 ztc_f = ztc_f - zfact1 * ( sbc_trc(ji,jj,jn) - sbc_trc_b(ji,jj,jn) ) 237 ENDIF 238 239 ze3t_f = 1.e0 / ze3t_f 240 trb(ji,jj,jk,jn) = ztc_f * ze3t_f ! ptb <-- ptn filtered 241 trn(ji,jj,jk,jn) = tra(ji,jj,jk,jn) ! ptn <-- pta 242 ! 243 END DO 244 END DO 245 END DO 246 ! 247 END DO 248 ! 249 END SUBROUTINE trc_nxt_off 167 250 #else 168 251 !!---------------------------------------------------------------------- -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/TRP/trcrad.F90
r4990 r6952 64 64 IF( lk_c14b ) CALL trc_rad_sms( kt, trb, trn, jp_c14b0, jp_c14b1 ) ! bomb C14 65 65 IF( lk_pisces ) CALL trc_rad_sms( kt, trb, trn, jp_pcs0 , jp_pcs1, cpreserv='Y' ) ! PISCES model 66 IF( lk_my_trc ) CALL trc_rad_sms( kt, trb, trn, jp_myt0 , jp_myt1 66 IF( lk_my_trc ) CALL trc_rad_sms( kt, trb, trn, jp_myt0 , jp_myt1, cpreserv='Y' ) ! MY_TRC model 67 67 68 68 ! -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90
r6941 r6952 128 128 ! Coupling offline : runoff are in emp which contains E-P-R 129 129 ! 130 IF( .NOT. lk_offline .AND.lk_vvl ) THEN ! online coupling with vvl130 IF( lk_vvl ) THEN ! online coupling with vvl 131 131 zsfx(:,:) = 0._wp 132 132 ELSE ! online coupling free surface or offline with free surface -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/trcnam.F90
r6204 r6952 331 331 !!--------------------------------------------------------------------- 332 332 333 IF(lwp) WRITE(numout,*)334 IF(lwp) WRITE(numout,*) 'trc_nam_dia : read the passive tracer diagnostics options'335 IF(lwp) WRITE(numout,*) '~~~~~~~'336 337 333 IF(lwp) WRITE(numout,*) 338 334 IF(lwp) WRITE(numout,*) 'trc_nam_dia : read the passive tracer diagnostics options' -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/trcstp.F90
r6941 r6952 86 86 tra(:,:,:,:) = 0.e0 87 87 ! 88 88 IF( .NOT.lk_offline ) CALL trc_rst_opn ( kt ) ! Open tracer restart file 89 89 IF( lrst_trc ) CALL trc_rst_cal ( kt, 'WRITE' ) ! calendar 90 90 IF( lk_iomput ) THEN ; CALL trc_wri ( kt ) ! output of passive tracers with iom I/O manager
Note: See TracChangeset
for help on using the changeset viewer.