Changeset 7367
- Timestamp:
- 2016-11-29T11:52:31+01:00 (8 years ago)
- Location:
- branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 14 added
- 94 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r7363 r7367 10 10 !! NEMO 3.3 ! 2010-05 (D. Lea) Update to work with NEMO v3.2 11 11 !! - ! 2010-05 (D. Lea) add calc_month_len routine based on day_init 12 !! 3.4 ! 2012-10 (A. Weaver and K. Mogensen) Fix for direct initialization 12 13 !!---------------------------------------------------------------------- 13 14 … … 20 21 !! dyn_asm_inc : Apply the dynamic (u and v) increments 21 22 !! ssh_asm_inc : Apply the SSH increment 23 !! seaice_asm_inc : Apply the seaice increment 22 24 !!---------------------------------------------------------------------- 23 25 USE wrk_nemo ! Memory Allocation … … 25 27 USE dom_oce ! Ocean space and time domain 26 28 USE oce ! Dynamics and active tracers defined in memory 27 USE divcur ! Horizontal divergence and relative vorticity28 29 USE ldfdyn_oce ! ocean dynamics: lateral physics 29 30 USE eosbn2 ! Equation of state - in situ and potential density … … 33 34 USE c1d ! 1D initialization 34 35 USE in_out_manager ! I/O manager 35 USE lib_mpp ! MPP library 36 USE lib_mpp ! MPP library 37 #if defined key_lim3 || defined key_lim2 || defined key_cice 38 #if defined key_lim3 39 USE ice_3, ONLY : & ! LIM Ice model variables () 40 & frld, pfrld, hicif, hsnif, phicif 41 USE sbc_oce, ONLY : & 42 & fr_i ! ice fraction 43 #endif 44 #if defined key_lim2 45 USE ice_2, ONLY : & ! LIM Ice model variables 46 & frld, pfrld, hicif, hsnif, phicif 47 USE sbc_oce, ONLY : & 48 & fr_i ! ice fraction 49 #endif 50 #if defined key_cice 51 USE sbc_oce, ONLY : & 52 & fr_i ! ice fraction 53 USE sbc_ice, ONLY : & ! CICE Ice model variables 54 & naicet, ndaice_da, nfresh_da, nfsalt_da, nTf 55 USE ice_constants, only: Lfresh, rhoi,rhos ! for updating ice and snow enthalphy 56 ! USE ice_therm_itd, only: hfrazilmin ! thickness at new ice points 57 USE ice_domain_size, only: ncat,ntilyr,ntslyr 58 #endif 59 USE phycst, ONLY : & ! Physical Ice variables 60 & soce, sice, rhoic, rhosn, rday 61 #endif 62 USE sbc_oce ! Surface boundary condition variables. 63 64 USE eosbn2, only: tfreez 65 66 USE zdfmxl, ONLY : & 67 & hmld_tref, & 68 #if defined key_karaml 69 & hmld_kara, & 70 #endif 71 & hmld, & 72 & hmlp 73 74 #if defined key_bdy 75 USE bdy_oce, ONLY: bdytmask 76 #endif 77 USE histcom 36 78 37 79 IMPLICIT NONE … … 43 85 PUBLIC dyn_asm_inc !: Apply the dynamic (u and v) increments 44 86 PUBLIC ssh_asm_inc !: Apply the SSH increment 87 PUBLIC seaice_asm_inc !: Apply the seaice increment 45 88 46 89 #if defined key_asminc … … 56 99 LOGICAL, PUBLIC :: ln_dyninc = .FALSE. !: No dynamics (u and v) assimilation increments 57 100 LOGICAL, PUBLIC :: ln_sshinc = .FALSE. !: No sea surface height assimilation increment 101 LOGICAL, PUBLIC :: ln_seaiceinc = .FALSE. !: No sea ice concentration increment 58 102 LOGICAL, PUBLIC :: ln_salfix = .FALSE. !: Apply minimum salinity check 103 LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing 59 104 INTEGER, PUBLIC :: nn_divdmp = 0 !: Apply divergence damping filter nn_divdmp times 60 105 … … 78 123 79 124 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ssh_bkg, ssh_bkginc ! Background sea surface height and its increment 80 125 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: seaice_bkginc ! Increment to the background sea ice conc 126 127 INTEGER :: mld_choice = 4 !: choice of mld criteria to use 128 !: 1) turbocline depth 129 !: 2) surface to 0.001 kg/m^3 change 130 !: 3) Kara MLD 131 !: 4) Temperature criteria. 132 81 133 !! * Substitutions 82 134 # include "domzgr_substitute.h90" … … 122 174 REAL(wp) :: zdate_bkg ! Date in background state file for DI 123 175 REAL(wp) :: zdate_inc ! Time axis in increments file 124 176 177 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: & 178 & t_bkginc_2d !file for reading in 2D 179 !temperature increments 180 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: & 181 & z_mld !Mixed layer depth 182 125 183 REAL(wp), POINTER, DIMENSION(:,:) :: hdiv 184 126 185 !! 127 186 NAMELIST/nam_asminc/ ln_bkgwri, ln_trjwri, & … … 130 189 & nitbkg, nitdin, nitiaustr, nitiaufin, niaufn, & 131 190 & nittrjfrq, ln_salfix, salfixmin, & 132 & nn_divdmp 191 & nn_divdmp, mld_choice 133 192 !!---------------------------------------------------------------------- 134 193 … … 143 202 ln_dyninc = .FALSE. 144 203 ln_sshinc = .FALSE. 204 ln_seaiceinc = .FALSE. 145 205 ln_asmdin = .FALSE. 146 206 ln_asmiau = .TRUE. 147 207 ln_salfix = .FALSE. 208 ln_temnofreeze = .FALSE. 148 209 salfixmin = -9999 149 210 nitbkg = 0 … … 156 217 REWIND ( numnam ) 157 218 READ ( numnam, nam_asminc ) 158 219 220 ! Set the data time for diagnostics to the end of the IAU period 221 ! and multiply by the timestep to get seconds from start of run 222 data_time = rdt * nitiaufin 223 224 IF( ln_sco .AND. (ln_sshinc .OR. ln_seaiceinc .OR. ln_asmdin & 225 & .OR. ln_dyninc ) )THEN 226 CALL ctl_warn("Only SST assimilation currently supported in "//& 227 & "s-coordinates") 228 ln_sshinc = .FALSE. 229 ln_seaiceinc = .FALSE. 230 ln_asmdin = .FALSE. 231 ln_dyninc = .FALSE. 232 ENDIF 233 159 234 ! Control print 160 235 IF(lwp) THEN … … 169 244 WRITE(numout,*) ' Logical switch for applying SSH increments ln_sshinc = ', ln_sshinc 170 245 WRITE(numout,*) ' Logical switch for Direct Initialization (DI) ln_asmdin = ', ln_asmdin 246 WRITE(numout,*) ' Logical switch for applying sea ice increments ln_seaiceinc = ', ln_seaiceinc 171 247 WRITE(numout,*) ' Logical switch for Incremental Analysis Updating (IAU) ln_asmiau = ', ln_asmiau 172 248 WRITE(numout,*) ' Timestep of background in [0,nitend-nit000-1] nitbkg = ', nitbkg … … 235 311 236 312 IF ( ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) & 237 .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) )) &238 & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc and ln_sshinc is set to .true.', &313 .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) & 314 & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', & 239 315 & ' but ln_asmdin and ln_asmiau are both set to .false. :', & 240 316 & ' Inconsistent options') … … 248 324 & ' Type IAU weighting function is invalid') 249 325 250 IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ) &326 IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) & 251 327 & ) & 252 & CALL ctl_warn( ' ln_trainc, ln_dyninc and ln_sshinc are set to .false. :', &328 & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', & 253 329 & ' The assimilation increments are not applied') 254 330 … … 353 429 ALLOCATE( v_bkginc(jpi,jpj,jpk) ) 354 430 ALLOCATE( ssh_bkginc(jpi,jpj) ) 431 ALLOCATE( seaice_bkginc(jpi,jpj)) 355 432 #if defined key_asminc 356 433 ALLOCATE( ssh_iau(jpi,jpj) ) … … 361 438 v_bkginc(:,:,:) = 0.0 362 439 ssh_bkginc(:,:) = 0.0 440 seaice_bkginc(:,:) = 0.0 363 441 #if defined key_asminc 364 442 ssh_iau(:,:) = 0.0 365 443 #endif 366 IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) ) THEN444 IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN 367 445 368 446 !-------------------------------------------------------------------- … … 397 475 398 476 IF ( ln_trainc ) THEN 399 CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 ) 400 CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 ) 401 ! Apply the masks 402 t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:) 403 s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:) 404 ! Set missing increments to 0.0 rather than 1e+20 405 ! to allow for differences in masks 406 WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0 407 WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0 477 478 IF (ln_sco) THEN 479 480 ALLOCATE(z_mld(jpi,jpj)) 481 SELECT CASE(mld_choice) 482 CASE(1) 483 z_mld = hmld 484 CASE(2) 485 z_mld = hmlp 486 CASE(3) 487 #if defined key_karaml 488 z_mld = hmld_kara 489 #endif 490 CALL ctl_stop("Kara mixed layer not defined in current version of NEMO") ! JW: Safty feature, should be removed 491 ! once the kara mixed layer is availible 492 CASE(4) 493 z_mld = hmld_tref 494 END SELECT 495 496 ALLOCATE( t_bkginc_2d(jpi,jpj) ) 497 CALL iom_get( inum, jpdom_autoglo, 'bckinsurft', t_bkginc_2d, 1) 498 #if defined key_bdy 499 DO jk = 1,jpkm1 500 WHERE( z_mld(:,:) > fsdepw(:,:,jk) ) 501 t_bkginc(:,:,jk) = t_bkginc_2d(:,:) * bdytmask(:,:) 502 ELSEWHERE 503 t_bkginc(:,:,jk) = 0. 504 ENDWHERE 505 ENDDO 506 #else 507 t_bkginc(:,:,:) = 0. 508 #endif 509 s_bkginc(:,:,:) = 0. 510 511 !DEALLOCATE temporary arrays 512 DEALLOCATE(z_mld, t_bkginc_2d) 513 514 ELSE 515 516 CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 ) 517 CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 ) 518 ! Apply the masks 519 t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:) 520 s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:) 521 ! Set missing increments to 0.0 rather than 1e+20 522 ! to allow for differences in masks 523 WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0 524 WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0 525 526 ENDIF 527 408 528 ENDIF 409 529 … … 429 549 ENDIF 430 550 551 IF ( ln_seaiceinc ) THEN 552 CALL iom_get( inum, jpdom_autoglo, 'bckinseaice', seaice_bkginc, 1 ) 553 ! Apply the masks 554 seaice_bkginc(:,:) = seaice_bkginc(:,:) * tmask(:,:,1) 555 ! Set missing increments to 0.0 rather than 1e+20 556 ! to allow for differences in masks 557 WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0 558 ENDIF 559 431 560 CALL iom_close( inum ) 432 561 … … 437 566 !----------------------------------------------------------------------- 438 567 439 440 568 IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN 441 569 442 CALL wrk_alloc(jpi,jpj,hdiv) 443 444 DO jt = 1, nn_divdmp 445 446 DO jk = 1, jpkm1 447 448 hdiv(:,:) = 0._wp 449 450 DO jj = 2, jpjm1 451 DO ji = fs_2, fs_jpim1 ! vector opt. 452 hdiv(ji,jj) = & 453 ( e2u(ji ,jj)*fse3u(ji ,jj,jk) * u_bkginc(ji ,jj,jk) & 454 - e2u(ji-1,jj)*fse3u(ji-1,jj,jk) * u_bkginc(ji-1,jj,jk) & 455 + e1v(ji,jj )*fse3v(ji,jj ,jk) * v_bkginc(ji,jj ,jk) & 456 - e1v(ji,jj-1)*fse3v(ji,jj-1,jk) * v_bkginc(ji,jj-1,jk) ) & 457 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 570 CALL wrk_alloc(jpi,jpj,hdiv) 571 572 DO jt = 1, nn_divdmp 573 574 DO jk = 1, jpkm1 575 576 hdiv(:,:) = 0._wp 577 578 DO jj = 2, jpjm1 579 DO ji = fs_2, fs_jpim1 ! vector opt. 580 hdiv(ji,jj) = & 581 ( e2u(ji ,jj ) * fse3u(ji ,jj ,jk) * u_bkginc(ji ,jj ,jk) & 582 - e2u(ji-1,jj ) * fse3u(ji-1,jj ,jk) * u_bkginc(ji-1,jj ,jk) & 583 + e1v(ji ,jj ) * fse3v(ji ,jj ,jk) * v_bkginc(ji ,jj ,jk) & 584 - e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) * v_bkginc(ji ,jj-1,jk) ) & 585 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 586 END DO 458 587 END DO 588 589 CALL lbc_lnk( hdiv, 'T', 1. ) ! lateral boundary cond. (no sign change) 590 591 DO jj = 2, jpjm1 592 DO ji = fs_2, fs_jpim1 ! vector opt. 593 u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji+1,jj)*e2t(ji+1,jj) * hdiv(ji+1,jj) & 594 - e1t(ji ,jj)*e2t(ji ,jj) * hdiv(ji ,jj) ) & 595 / e1u(ji,jj) * umask(ji,jj,jk) 596 v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji,jj+1)*e2t(ji,jj+1) * hdiv(ji,jj+1) & 597 - e1t(ji,jj )*e2t(ji,jj ) * hdiv(ji,jj ) ) & 598 / e2v(ji,jj) * vmask(ji,jj,jk) 599 END DO 600 END DO 601 459 602 END DO 460 603 461 CALL lbc_lnk( hdiv, 'T', 1. ) ! lateral boundary cond. (no sign change) 462 463 DO jj = 2, jpjm1 464 DO ji = fs_2, fs_jpim1 ! vector opt. 465 u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2 * ( e1t(ji+1,jj)*e2t(ji+1,jj) * hdiv(ji+1,jj) & 466 - e1t(ji ,jj)*e2t(ji ,jj) * hdiv(ji ,jj) ) & 467 / e1u(ji,jj) * umask(ji,jj,jk) 468 v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2 * ( e1t(ji,jj+1)*e2t(ji,jj+1) * hdiv(ji,jj+1) & 469 - e1t(ji,jj )*e2t(ji,jj ) * hdiv(ji,jj ) ) & 470 / e2v(ji,jj) * vmask(ji,jj,jk) 471 END DO 472 END DO 473 474 END DO 475 476 END DO 477 478 CALL wrk_dealloc(jpi,jpj,hdiv) 604 END DO 605 606 CALL wrk_dealloc(jpi,jpj,hdiv) 479 607 480 608 ENDIF … … 506 634 CALL iom_open( c_asmdin, inum ) 507 635 508 CALL iom_get( inum, ' zdate', zdate_bkg )636 CALL iom_get( inum, 'rdastp', zdate_bkg ) 509 637 510 638 IF(lwp) THEN … … 662 790 INTEGER :: it 663 791 REAL(wp) :: zincwgt ! IAU weight for current time step 664 !!---------------------------------------------------------------------- 792 REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values 793 !!---------------------------------------------------------------------- 794 795 ! freezing point calculation taken from oc_fz_pt (but calculated for all depths) 796 ! used to prevent the applied increments taking the temperature below the local freezing point 797 798 ! Note: For NEMO/CICE this will be identical to nTf (for the surface), but defined at the now point. 799 800 DO jk=1, jpkm1 801 fzptnz (:,:,jk) = tfreez(tsn(:,:,jk,jp_sal )) - 7.53e-4_wp * fsdepw(:,:,jk) 802 ENDDO 665 803 666 804 IF ( ln_asmiau ) THEN … … 684 822 ! Update the tracer tendencies 685 823 DO jk = 1, jpkm1 686 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 687 tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt 824 IF (ln_temnofreeze) THEN 825 ! Do not apply negative increments if the temperature will fall below freezing 826 WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. & 827 & tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 828 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 829 END WHERE 830 ELSE 831 tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 832 ENDIF 833 IF (ln_salfix) THEN 834 ! Do not apply negative increments if the salinity will fall below a specified 835 ! minimum value salfixmin 836 WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. & 837 & tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 838 tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt 839 END WHERE 840 ELSE 841 tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt 842 ENDIF 688 843 END DO 689 690 ! Salinity fix691 IF (ln_salfix) THEN692 DO jk = 1, jpkm1693 DO jj = 1, jpj694 DO ji= 1, jpi695 tsa(ji,jj,jk,jp_sal) = MAX( tsa(ji,jj,jk,jp_sal), salfixmin )696 END DO697 END DO698 END DO699 ENDIF700 844 701 845 ENDIF … … 718 862 719 863 ! Initialize the now fields with the background + increment 720 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:) 721 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:) 722 723 ! Optional salinity fix 864 IF (ln_temnofreeze) THEN 865 ! Do not apply negative increments if the temperature will fall below freezing 866 WHERE(t_bkginc(:,:,:) > 0.0_wp .OR. & 867 & tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 868 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:) 869 END WHERE 870 ELSE 871 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:) 872 ENDIF 724 873 IF (ln_salfix) THEN 725 DO jk = 1, jpkm1 726 DO jj = 1, jpj 727 DO ji= 1, jpi 728 tsn(ji,jj,jk,jp_sal) = MAX( tsn(ji,jj,jk,jp_sal), salfixmin ) 729 END DO 730 END DO 731 END DO 874 ! Do not apply negative increments if the salinity will fall below a specified 875 ! minimum value salfixmin 876 WHERE(s_bkginc(:,:,:) > 0.0_wp .OR. & 877 & tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 878 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:) 879 END WHERE 880 ELSE 881 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:) 732 882 ENDIF 733 883 734 tsb(:,:,:,:) = tsn(:,:,:,:) 884 tsb(:,:,:,:) = tsn(:,:,:,:) ! Update before fields 735 885 736 886 CALL eos( tsb, rhd, rhop ) ! Before potential and in situ densities 737 887 738 888 IF( ln_zps .AND. .NOT. lk_c1d ) & 739 & CALL zps_hde( nit000, jpts, tsb, 740 & gtsu, gtsv, rhd, 889 & CALL zps_hde( nit000, jpts, tsb, & ! Partial steps: before horizontal derivative 890 & gtsu, gtsv, rhd, & ! of T, S, rd at the bottom ocean level 741 891 & gru , grv ) 892 893 #if defined key_zdfkpp 894 CALL eos( tsn, rhd ) ! Compute rhd 895 #endif 742 896 743 897 DEALLOCATE( t_bkginc ) … … 748 902 ! 749 903 ENDIF 904 ! Perhaps the following call should be in step 905 IF ( ln_seaiceinc ) CALL seaice_asm_inc ( kt ) ! apply sea ice concentration increment 750 906 ! 751 907 END SUBROUTINE tra_asm_inc … … 817 973 vb(:,:,:) = vn(:,:,:) 818 974 819 CALL div_cur( kt ) ! Compute divergence and curl for now fields820 821 rotb (:,:,:) = rotn (:,:,:) ! Update before fields822 hdivb(:,:,:) = hdivn(:,:,:)823 824 975 DEALLOCATE( u_bkg ) 825 976 DEALLOCATE( v_bkg ) … … 846 997 ! 847 998 INTEGER :: it 999 INTEGER :: jk 848 1000 REAL(wp) :: zincwgt ! IAU weight for current time step 849 1001 !!---------------------------------------------------------------------- … … 891 1043 sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) 892 1044 893 sshb(:,:) = sshn(:,:) ! Update before fields 1045 ! Update before fields 1046 sshb(:,:) = sshn(:,:) 894 1047 895 1048 DEALLOCATE( ssh_bkg ) … … 902 1055 END SUBROUTINE ssh_asm_inc 903 1056 1057 SUBROUTINE seaice_asm_inc( kt, kindic ) 1058 !!---------------------------------------------------------------------- 1059 !! *** ROUTINE seaice_asm_inc *** 1060 !! 1061 !! ** Purpose : Apply the sea ice assimilation increment. 1062 !! 1063 !! ** Method : Direct initialization or Incremental Analysis Updating. 1064 !! 1065 !! ** Action : 1066 !! 1067 !! History : 1068 !! ! 07-2011 (D. Lea) Initial version based on ssh_asm_inc 1069 !!---------------------------------------------------------------------- 1070 1071 IMPLICIT NONE 1072 1073 !! * Arguments 1074 INTEGER, INTENT(IN) :: kt ! Current time step 1075 INTEGER, OPTIONAL, INTENT(IN) :: kindic ! flag for disabling the deallocation 1076 1077 !! * Local declarations 1078 INTEGER :: it 1079 REAL(wp) :: zincwgt ! IAU weight for current time step 1080 1081 #if defined key_lim3 || defined key_lim2 1082 REAL(wp), DIMENSION(jpi,jpj) :: zofrld, zohicif, zseaicendg, zhicifinc ! LIM 1083 REAL(wp) :: zhicifmin=0.5_wp ! ice minimum depth in metres 1084 1085 #endif 1086 1087 1088 IF ( ln_asmiau ) THEN 1089 1090 !-------------------------------------------------------------------- 1091 ! Incremental Analysis Updating 1092 !-------------------------------------------------------------------- 1093 1094 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 1095 1096 it = kt - nit000 + 1 1097 zincwgt = wgtiau(it) ! IAU weight for the current time step 1098 ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments) 1099 1100 IF(lwp) THEN 1101 WRITE(numout,*) 1102 WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', & 1103 & kt,' with IAU weight = ', wgtiau(it) 1104 WRITE(numout,*) '~~~~~~~~~~~~' 1105 ENDIF 1106 1107 #if defined key_lim3 || defined key_lim2 1108 1109 zofrld(:,:)=frld(:,:) 1110 zohicif(:,:)=hicif(:,:) 1111 1112 frld = MIN( MAX( frld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 1113 pfrld = MIN( MAX( pfrld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 1114 fr_i(:,:) = 1.0_wp - frld(:,:) ! adjust ice fraction 1115 1116 zseaicendg(:,:)=zofrld(:,:) - frld(:,:) ! find out actual sea ice nudge applied 1117 1118 ! Nudge sea ice depth to bring it up to a required minimum depth 1119 1120 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 1121 zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt 1122 ELSEWHERE 1123 zhicifinc(:,:) = 0.0_wp 1124 END WHERE 1125 1126 ! nudge ice depth 1127 hicif(:,:)=hicif(:,:) + zhicifinc(:,:) 1128 phicif(:,:)=phicif(:,:) + zhicifinc(:,:) 1129 1130 ! seaice salinity balancing (to add) 1131 1132 #endif 1133 1134 #if defined key_cice 1135 1136 ! Pass ice increment tendency into CICE 1137 ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt 1138 1139 #endif 1140 1141 IF ( kt == nitiaufin_r ) THEN 1142 DEALLOCATE( seaice_bkginc ) 1143 ENDIF 1144 1145 ELSE 1146 1147 #if defined key_cice 1148 1149 ! Zero ice increment tendency into CICE 1150 ndaice_da(:,:) = 0.0_wp 1151 1152 #endif 1153 1154 ENDIF 1155 1156 ELSEIF ( ln_asmdin ) THEN 1157 1158 !-------------------------------------------------------------------- 1159 ! Direct Initialization 1160 !-------------------------------------------------------------------- 1161 1162 IF ( kt == nitdin_r ) THEN 1163 1164 neuler = 0 ! Force Euler forward step 1165 1166 #if defined key_lim3 || defined key_lim2 1167 1168 zofrld(:,:)=frld(:,:) 1169 zohicif(:,:)=hicif(:,:) 1170 1171 ! Initialize the now fields the background + increment 1172 1173 frld(:,:) = MIN( MAX( frld(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp) 1174 pfrld(:,:) = frld(:,:) 1175 fr_i(:,:) = 1.0_wp - frld(:,:) ! adjust ice fraction 1176 1177 zseaicendg(:,:)=zofrld(:,:) - frld(:,:) ! find out actual sea ice nudge applied 1178 1179 ! Nudge sea ice depth to bring it up to a required minimum depth 1180 1181 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 1182 zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt 1183 ELSEWHERE 1184 zhicifinc(:,:) = 0.0_wp 1185 END WHERE 1186 1187 ! nudge ice depth 1188 hicif(:,:)=hicif(:,:) + zhicifinc(:,:) 1189 phicif(:,:)=phicif(:,:) 1190 1191 ! seaice salinity balancing (to add) 1192 1193 #endif 1194 1195 #if defined key_cice 1196 1197 ! Pass ice increment tendency into CICE - is this correct? 1198 ndaice_da(:,:) = seaice_bkginc(:,:) / rdt 1199 1200 #endif 1201 IF ( .NOT. PRESENT(kindic) ) THEN 1202 DEALLOCATE( seaice_bkginc ) 1203 END IF 1204 1205 ELSE 1206 1207 #if defined key_cice 1208 1209 ! Zero ice increment tendency into CICE 1210 ndaice_da(:,:) = 0.0_wp 1211 1212 #endif 1213 1214 ENDIF 1215 1216 !#if defined key_lim3 || defined key_lim2 || defined key_cice 1217 ! 1218 ! IF (ln_seaicebal ) THEN 1219 ! !! balancing salinity increments 1220 ! !! simple case from limflx.F90 (doesn't include a mass flux) 1221 ! !! assumption is that as ice concentration is reduced or increased 1222 ! !! the snow and ice depths remain constant 1223 ! !! note that snow is being created where ice concentration is being increased 1224 ! !! - could be more sophisticated and 1225 ! !! not do this (but would need to alter h_snow) 1226 ! 1227 ! usave(:,:,:)=sb(:,:,:) ! use array as a temporary store 1228 ! 1229 ! DO jj = 1, jpj 1230 ! DO ji = 1, jpi 1231 ! ! calculate change in ice and snow mass per unit area 1232 ! ! positive values imply adding salt to the ocean (results from ice formation) 1233 ! ! fwf : ice formation and melting 1234 ! 1235 ! zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rdt 1236 ! 1237 ! ! change salinity down to mixed layer depth 1238 ! mld=hmld_kara(ji,jj) 1239 ! 1240 ! ! prevent small mld 1241 ! ! less than 10m can cause salinity instability 1242 ! IF (mld < 10) mld=10 1243 ! 1244 ! ! set to bottom of a level 1245 ! DO jk = jpk-1, 2, -1 1246 ! IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN 1247 ! mld=gdepw(ji,jj,jk+1) 1248 ! jkmax=jk 1249 ! ENDIF 1250 ! ENDDO 1251 ! 1252 ! ! avoid applying salinity balancing in shallow water or on land 1253 ! ! 1254 ! 1255 ! ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m) 1256 ! 1257 ! dsal_ocn=0.0_wp 1258 ! sal_thresh=5.0_wp ! minimum salinity threshold for salinity balancing 1259 ! 1260 ! if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) & 1261 ! dsal_ocn = zfons / (rhop(ji,jj,1) * mld) 1262 ! 1263 ! ! put increments in for levels in the mixed layer 1264 ! ! but prevent salinity below a threshold value 1265 ! 1266 ! DO jk = 1, jkmax 1267 ! 1268 ! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 1269 ! sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn 1270 ! sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn 1271 ! ENDIF 1272 ! 1273 ! ENDDO 1274 ! 1275 ! ! ! salt exchanges at the ice/ocean interface 1276 ! ! zpmess = zfons / rdt_ice ! rdt_ice is ice timestep 1277 ! ! 1278 ! !! Adjust fsalt. A +ve fsalt means adding salt to ocean 1279 ! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt 1280 ! !! 1281 ! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d) 1282 ! !! ! E-P (kg m-2 s-2) 1283 ! ! emp(ji,jj) = emp(ji,jj) + zpmess ! E-P (kg m-2 s-2) 1284 ! ENDDO !ji 1285 ! ENDDO !jj! 1286 ! 1287 ! ENDIF !ln_seaicebal 1288 ! 1289 !#endif 1290 1291 1292 ENDIF 1293 1294 END SUBROUTINE seaice_asm_inc 904 1295 !!====================================================================== 905 1296 END MODULE asminc -
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/ASM/asmtrj.F90
r7363 r7367 14 14 !! ! 2009-06 (F. Vigilant) asm_trj_wri: special case when kt=nit000-1 15 15 !! ! 2009-07 (F. Vigilant) asm_trj_wri: add computation of eiv at restart 16 17 !! ! 2012-11 (A. Weaver) Save avt_bkg for mixing layer computation, remove en_bkg 16 18 !!---------------------------------------------------------------------- 17 19 … … 35 37 USE zdfmxl ! Mixed layer depth 36 38 USE dom_oce, ONLY : ndastp 37 USE sol_oce, ONLY : gcx ! Solver variables defined in memory38 39 USE in_out_manager ! I/O manager 39 40 USE iom ! I/O module … … 43 44 USE ldfeiv ! eddy induced velocity coef. (ldf_eiv routine) 44 45 #endif 45 46 #if defined key_lim2 47 USE ice_2 48 #endif 49 #if defined key_lim3 50 USE ice 51 #endif 46 52 IMPLICIT NONE 47 53 PRIVATE … … 110 116 CALL iom_rstput( kt, nitbkg_r, inum, 'tn' , tsn(:,:,:,jp_tem) ) 111 117 CALL iom_rstput( kt, nitbkg_r, inum, 'sn' , tsn(:,:,:,jp_sal) ) 118 CALL iom_rstput( kt, nitbkg_r, inum, 'avt' , avt ) 112 119 CALL iom_rstput( kt, nitbkg_r, inum, 'sshn' , sshn ) 113 #if defined key_zdftke114 CALL iom_rstput( kt, nitbkg_r, inum, 'en' , en )115 #endif116 CALL iom_rstput( kt, nitbkg_r, inum, 'gcx' , gcx )117 120 ! 118 121 CALL iom_close( inum ) … … 148 151 CALL iom_rstput( kt, nitdin_r, inum, 'tn' , tsn(:,:,:,jp_tem) ) 149 152 CALL iom_rstput( kt, nitdin_r, inum, 'sn' , tsn(:,:,:,jp_sal) ) 153 CALL iom_rstput( kt, nitbkg_r, inum, 'avt' , avt ) 150 154 CALL iom_rstput( kt, nitdin_r, inum, 'sshn' , sshn ) 155 #if defined key_lim2 || defined key_lim3 156 IF(( nn_ice == 2 ) .OR. ( nn_ice == 3 )) THEN 157 CALL iom_rstput( kt, nitdin_r, inum, 'iceconc', 1.0 - frld(:,:) ) 158 ENDIF 159 #endif 151 160 ! 152 161 CALL iom_close( inum ) … … 222 231 CALL iom_rstput( it, it, inum, 'avt' , avt ) 223 232 #if defined key_ldfslp 224 CALL iom_rstput( it, it, inum, 'uslp' , uslp ) 225 CALL iom_rstput( it, it, inum, 'vslp' , vslp ) 226 CALL iom_rstput( it, it, inum, 'wslpi' , wslpi ) 227 CALL iom_rstput( it, it, inum, 'wslpj' , wslpj ) 233 CALL iom_rstput( it, it, inum, 'uslp_hor' , uslp_hor ) 234 CALL iom_rstput( it, it, inum, 'vslp_hor' , vslp_hor ) 235 CALL iom_rstput( it, it, inum, 'wslpi_hor' , wslpi_hor ) 236 CALL iom_rstput( it, it, inum, 'wslpj_hor' , wslpj_hor ) 237 CALL iom_rstput( it, it, inum, 'uslp_iso' , uslp_iso ) 238 CALL iom_rstput( it, it, inum, 'vslp_iso' , vslp_iso ) 239 CALL iom_rstput( it, it, inum, 'wslpi_iso' , wslpi_iso ) 240 CALL iom_rstput( it, it, inum, 'wslpj_iso' , wslpj_iso ) 228 241 #endif 229 242 #if defined key_zdfddm -
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
r7363 r7367 57 57 LOGICAL :: ln_mask_file !: =T read bdymask from file 58 58 LOGICAL :: ln_vol !: =T volume correction 59 LOGICAL, DIMENSION(jp_bdy) :: ln_sponge !: =T use sponge layer 59 60 ! 60 61 INTEGER :: nb_bdy !: number of open boundary sets … … 62 63 INTEGER :: nn_volctl !: = 0 the total volume will have the variability of the surface Flux E-P 63 64 ! ! = 1 the volume will be constant during all the integration. 65 REAL(wp) :: rn_sponge !: multiplier of diffusion for sponge layer 66 64 67 INTEGER, DIMENSION(jp_bdy) :: nn_dyn2d ! Choice of boundary condition for barotropic variables (U,V,SSH) 65 68 INTEGER, DIMENSION(jp_bdy) :: nn_dyn2d_dta !: = 0 use the initial state as bdy dta ; … … 83 86 !! Global variables 84 87 !!---------------------------------------------------------------------- 85 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: bdytmask !: Mask defining computational domain at T-points 86 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: bdyumask !: Mask defining computational domain at U-points 87 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: bdyvmask !: Mask defining computational domain at V-points 88 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: bdytmask !: Mask defining computational domain at T-points 89 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: bdyumask !: Mask defining computational domain at U-points 90 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: bdyvmask !: Mask defining computational domain at V-points 91 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sponge_factor !: Multiplier for diffusion for sponge layer 88 92 89 93 REAL(wp) :: bdysurftot !: Lateral surface of unstructured open boundary … … 120 124 ! 121 125 ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj), & 122 & STAT=bdy_oce_alloc )126 & sponge_factor(jpi,jpj), STAT=bdy_oce_alloc ) 123 127 ! 124 128 IF( lk_mpp ) CALL mpp_sum ( bdy_oce_alloc ) -
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
r7363 r7367 32 32 USE ice_2 33 33 #endif 34 USE sbcapr 34 35 35 36 IMPLICIT NONE … … 238 239 ENDIF 239 240 ENDIF 240 jstart = j end+1241 jstart = jstart + nb_bdy_fld(ib_bdy) 241 242 242 243 ! If full velocities in boundary data then split into barotropic and baroclinic data … … 281 282 END IF ! nn_dta(ib_bdy) = 1 282 283 END DO ! ib_bdy 284 285 IF ( ln_apr_obc ) THEN 286 DO ib_bdy = 1, nb_bdy 287 IF (nn_tra(ib_bdy).NE.4)THEN 288 igrd = 1 ! meridional velocity 289 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 290 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 291 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 292 dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + ssh_ib(ii,ij) 293 ENDDO 294 ENDIF 295 ENDDO 296 ENDIF 283 297 284 298 IF( nn_timing == 1 ) CALL timing_stop('bdy_dta') … … 317 331 TYPE(FLD_N) :: bn_frld, bn_hicif, bn_hsnif ! 318 332 #endif 319 NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d 320 #if defined key_lim2 321 NAMELIST/nambdy_dta/ bn_frld, bn_hicif, bn_hsnif 322 #endif 323 NAMELIST/nambdy_dta/ ln_full_vel 333 NAMELIST/nambdy_dta_1/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d 334 NAMELIST/nambdy_dta_2/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d 335 #if defined key_lim2 336 NAMELIST/nambdy_dta_1/ bn_frld, bn_hicif, bn_hsnif 337 NAMELIST/nambdy_dta_2/ bn_frld, bn_hicif, bn_hsnif 338 #endif 339 NAMELIST/nambdy_dta_1/ ln_full_vel 340 NAMELIST/nambdy_dta_2/ ln_full_vel 324 341 !!--------------------------------------------------------------------------- 325 342 … … 403 420 404 421 ! Important NOT to rewind here. 405 READ( numnam, nambdy_dta ) 422 if ( ib_bdy == 1 ) READ( numnam, nambdy_dta_1 ) 423 if ( ib_bdy == 2 ) READ( numnam, nambdy_dta_2 ) 406 424 407 425 cn_dir_array(ib_bdy) = cn_dir … … 554 572 ! Recalculate field counts 555 573 !------------------------- 556 nb_bdy_fld_sum = 0557 574 IF( ib_bdy .eq. 1 ) THEN 575 nb_bdy_fld_sum = 0 558 576 nb_bdy_fld(ib_bdy) = jfld 559 577 nb_bdy_fld_sum = jfld -
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim2.F90
r7363 r7367 53 53 CYCLE 54 54 CASE(jp_frs) 55 CALL bdy_ice_frs( idx_bdy(ib_bdy), dta_ idx(ib_bdy) )55 CALL bdy_ice_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy) ) 56 56 CASE DEFAULT 57 57 CALL ctl_stop( 'bdy_ice_lim_2 : unrecognised option for open boundaries for ice fields' ) -
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
r7363 r7367 83 83 & nn_ice_lim2, nn_ice_lim2_dta, & 84 84 #endif 85 & ln_vol, nn_volctl, nn_rimwidth85 & ln_vol, nn_volctl, ln_sponge, rn_sponge, nn_rimwidth 86 86 !! 87 87 NAMELIST/nambdy_index/ nbdysege, jpieob, jpjedt, jpjeft, & … … 127 127 ln_vol = .false. 128 128 nn_volctl = -1 ! uninitialised flag 129 ln_sponge(:) = .false. 130 rn_sponge = 0.0 129 131 nn_rimwidth(:) = -1 ! uninitialised flag 130 132 … … 224 226 IF(lwp) WRITE(numout,*) 225 227 228 IF( ln_sponge(ib_bdy) ) THEN ! check sponge layer choice 229 IF(lwp) WRITE(numout,*) 'Sponge layer applied at open boundaries' 230 IF(lwp) WRITE(numout,*) 'Multiplier for diffusion in sponge layer : ', rn_sponge 231 IF(lwp) WRITE(numout,*) 232 ELSE 233 IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 234 IF(lwp) WRITE(numout,*) 235 ENDIF 236 226 237 ENDDO 227 238 228 IF( ln_vol ) THEN ! check volume conservation (nn_volctl value) 229 IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 230 IF(lwp) WRITE(numout,*) 231 SELECT CASE ( nn_volctl ) 232 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' The total volume will be constant' 233 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' The total volume will vary according to the surface E-P flux' 234 CASE DEFAULT ; CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 235 END SELECT 236 IF(lwp) WRITE(numout,*) 237 ELSE 238 IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 239 IF(lwp) WRITE(numout,*) 240 ENDIF 239 IF( ln_vol ) THEN ! check volume conservation (nn_volctl value) 240 IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries' 241 IF(lwp) WRITE(numout,*) 242 SELECT CASE ( nn_volctl ) 243 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' The total volume will be constant' 244 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' The total volume will vary according to the surface E-P flux' 245 CASE DEFAULT ; CALL ctl_stop( 'nn_volctl must be 0 or 1' ) 246 END SELECT 247 IF(lwp) WRITE(numout,*) 248 ELSE 249 IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries' 250 IF(lwp) WRITE(numout,*) 251 ENDIF 252 253 sponge_factor(:,:) = 1.0 241 254 242 255 ! ------------------------------------------------- … … 247 260 ! --------------------------------------------- 248 261 REWIND( numnam ) 262 jpbdta = 1 249 263 DO ib_bdy = 1, nb_bdy 250 264 251 jpbdta = 1252 265 IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters 253 266 … … 317 330 318 331 nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy) 319 jpbdta = MAX VAL(nblendta(:,ib_bdy))332 jpbdta = MAX( jpbdta, MAXVAL(nblendta(:,ib_bdy)) ) 320 333 321 334 … … 324 337 325 338 CALL iom_open( cn_coords_file(ib_bdy), inum ) 326 jpbdta = 1 339 327 340 DO igrd = 1, jpbgrd 328 341 id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) … … 330 343 jpbdta = MAX(jpbdta, kdimsz(1)) 331 344 ENDDO 345 CALL iom_close( inum ) 332 346 333 347 ENDIF … … 507 521 ELSE ! Read global index arrays from boundary coordinates file. 508 522 523 CALL iom_open( cn_coords_file(ib_bdy), inum ) 509 524 DO igrd = 1, jpbgrd 510 525 CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) ) … … 616 631 END DO 617 632 618 ENDDO 633 ! Compute multiplier for diffusion for sponge layer 634 ! ------------------------------------------------- 635 IF( ln_sponge(ib_bdy) ) THEN 636 igrd=1 637 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 638 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 639 nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 640 nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 641 sponge_factor(nbi,nbj) = 1.0 + (rn_sponge-1.0) * ( 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) ) 642 END DO 643 ENDIF 644 645 ENDDO ! ib_bdy 619 646 620 647 ! ------------------------------------------------------ … … 773 800 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 774 801 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 775 nbj => idx_bdy(ib_bdy)%nb i(ib,igrd)802 nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 776 803 flagu => idx_bdy(ib_bdy)%flagu(ib) 777 804 bdysurftot = bdysurftot + hu (nbi , nbj) & … … 786 813 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 787 814 nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 788 nbj => idx_bdy(ib_bdy)%nb i(ib,igrd)815 nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 789 816 flagv => idx_bdy(ib_bdy)%flagv(ib) 790 817 bdysurftot = bdysurftot + hv (nbi, nbj ) & -
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90
r3294 r7367 38 38 USE dianam ! build name of file 39 39 USE lib_mpp ! distributed memory computing library 40 #if defined key_lim2 || defined key_lim3 41 USE ice 40 #if defined key_lim2 41 USE ice_2 42 #endif 43 #if defined key_lim3 44 USE ice_3 42 45 #endif 43 46 USE domvvl … … 49 52 50 53 !! * Routine accessibility 51 PUBLIC dia_dct ! routine called by step.F90 52 PUBLIC dia_dct_init! routine called by opa.F90 54 PUBLIC dia_dct ! routine called by step.F90 55 PUBLIC dia_dct_init ! routine called by opa.F90 56 PUBLIC diadct_alloc ! routine called by nemo_init in nemogcm.F90 53 57 PRIVATE readsec 54 58 PRIVATE removepoints 55 59 PRIVATE transport 56 60 PRIVATE dia_dct_wri 61 PRIVATE dia_dct_wri_NOOS 57 62 58 63 #include "domzgr_substitute.h90" … … 65 70 INTEGER :: nn_dctwri = 1 ! Frequency of output 66 71 INTEGER :: nn_secdebug = 0 ! Number of the section to debug 72 INTEGER :: nn_dct_h = 1 ! Frequency of computation for NOOS hourly files 73 INTEGER :: nn_dctwri_h = 1 ! Frequency of output for NOOS hourly files 67 74 68 INTEGER, PARAMETER :: nb_class_max = 10 69 INTEGER, PARAMETER :: nb_sec_max = 150 70 INTEGER, PARAMETER :: nb_point_max = 2000 71 INTEGER, PARAMETER :: nb_type_class = 14 75 INTEGER, PARAMETER :: nb_class_max = 11 ! maximum number of classes, i.e. depth levels or density classes 76 INTEGER, PARAMETER :: nb_sec_max = 30 ! maximum number of sections 77 INTEGER, PARAMETER :: nb_point_max = 375 ! maximum number of points in a single section 78 INTEGER, PARAMETER :: nb_type = 14 ! types of calculations, i.e. pos transport, neg transport, heat transport, salt transport 79 INTEGER, PARAMETER :: nb_3d_vars = 5 80 INTEGER, PARAMETER :: nb_2d_vars = 2 72 81 INTEGER :: nb_sec 73 82 … … 82 91 TYPE SECTION 83 92 CHARACTER(len=60) :: name ! name of the sec 84 LOGICAL :: llstrpond ! true if you want the computation of salt and 85 ! heat transports 93 LOGICAL :: llstrpond ! true if you want the computation of salinity and heat transports 86 94 LOGICAL :: ll_ice_section ! ice surface and ice volume computation 87 95 LOGICAL :: ll_date_line ! = T if the section crosses the date-line … … 95 103 ztem ,&! temperature classes(99 if you don't want) 96 104 zlay ! level classes (99 if you don't want) 97 REAL(wp), DIMENSION(nb_type_class,nb_class_max) :: transport ! transport output 105 REAL(wp), DIMENSION(nb_type,nb_class_max) :: transport ! transport output 106 REAL(wp), DIMENSION(nb_type,nb_class_max) :: transport_h ! transport output 98 107 REAL(wp) :: slopeSection ! slope of the section 99 108 INTEGER :: nb_point ! number of points in the section … … 102 111 103 112 TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections 104 113 114 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: transports_3d 115 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: transports_2d 116 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: transports_3d_h 117 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: transports_2d_h 118 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_hr_output 105 119 106 120 CONTAINS 121 122 INTEGER FUNCTION diadct_alloc() 123 !!---------------------------------------------------------------------- 124 !! *** FUNCTION diadct_alloc *** 125 !!---------------------------------------------------------------------- 126 INTEGER :: ierr(2) 127 !!---------------------------------------------------------------------- 128 ! 129 ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(1) ) 130 ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max) , STAT=ierr(2) ) 131 ALLOCATE(transports_3d_h(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(3) ) 132 ALLOCATE(transports_2d_h(nb_2d_vars,nb_sec_max,nb_point_max) , STAT=ierr(4) ) 133 ALLOCATE(z_hr_output(nb_sec_max,24,nb_class_max) , STAT=ierr(5) ) 134 ! 135 diadct_alloc = MAXVAL( ierr ) 136 IF( diadct_alloc /= 0 ) CALL ctl_warn('diadct_alloc: failed to allocate arrays') 137 ! 138 END FUNCTION diadct_alloc 139 107 140 108 141 SUBROUTINE dia_dct_init … … 110 143 !! *** ROUTINE diadct *** 111 144 !! 112 !! ** Purpose: Read the namelist paramet res145 !! ** Purpose: Read the namelist parameters 113 146 !! Open output files 114 147 !! … … 121 154 REWIND( numnam ) 122 155 READ ( numnam, namdct ) 156 157 IF( ln_NOOS ) THEN 158 nn_dct=3600./rdt ! hard coded for NOOS transects, to give 25 hour means 159 nn_dctwri=86400./rdt 160 nn_dct_h=1 ! hard coded for NOOS transects, to give hourly data 161 nn_dctwri_h=3600./rdt 162 ENDIF 123 163 124 164 IF( lwp ) THEN … … 126 166 WRITE(numout,*) "diadct_init: compute transports through sections " 127 167 WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~" 128 WRITE(numout,*) " Frequency of computation: nn_dct = ",nn_dct 129 WRITE(numout,*) " Frequency of write: nn_dctwri = ",nn_dctwri 168 IF( ln_NOOS ) THEN 169 WRITE(numout,*) " Frequency of computation hard coded to be every hour: nn_dct = ",nn_dct 170 WRITE(numout,*) " Frequency of write hard coded to average 25 instantaneous hour values: nn_dctwri = ",nn_dctwri 171 WRITE(numout,*) " Frequency of hourly computation hard coded to be every timestep: nn_dct_h = ",nn_dct_h 172 WRITE(numout,*) " Frequency of hourly write hard coded to every hour: nn_dctwri_h = ",nn_dctwri_h 173 ELSE 174 WRITE(numout,*) " Frequency of computation: nn_dct = ",nn_dct 175 WRITE(numout,*) " Frequency of write: nn_dctwri = ",nn_dctwri 176 ENDIF 130 177 131 178 IF ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN … … 146 193 !open output file 147 194 IF( lwp ) THEN 148 CALL ctl_opn( numdct_vol, 'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 149 CALL ctl_opn( numdct_heat, 'heat_transport' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 150 CALL ctl_opn( numdct_salt, 'salt_transport' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 195 IF( ln_NOOS ) THEN 196 CALL ctl_opn( numdct_NOOS ,'NOOS_transport' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 197 CALL ctl_opn( numdct_NOOS_h,'NOOS_transport_h', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 198 ELSE 199 CALL ctl_opn( numdct_vol , 'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 200 CALL ctl_opn( numdct_temp, 'heat_transport' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 201 CALL ctl_opn( numdct_sal , 'salt_transport' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 202 ENDIF 151 203 ENDIF 204 205 ! Initialise arrays to zero 206 transports_3d(:,:,:,:) =0._wp 207 transports_2d(:,:,:) =0._wp 208 transports_3d_h(:,:,:,:)=0._wp 209 transports_2d_h(:,:,:) =0._wp 210 z_hr_output(:,:,:) =0._wp 152 211 153 212 IF( nn_timing == 1 ) CALL timing_stop('dia_dct_init') … … 160 219 !! *** ROUTINE diadct *** 161 220 !! 162 !! ** Purpose: Compute sections tranport and write it in numdct file 221 !! Purpose :: Compute section transports and write it in numdct files 222 !! 223 !! Method :: All arrays initialised to zero in dct_init 224 !! Each nn_dct time step call subroutine 'transports' for 225 !! each section to sum the transports. 226 !! Each nn_dctwri time step: 227 !! Divide the arrays by the number of summations to gain 228 !! an average value 229 !! Call dia_dct_sum to sum relevant grid boxes to obtain 230 !! totals for each class (density, depth, temp or sal) 231 !! Call dia_dct_wri to write the transports into file 232 !! Reinitialise all relevant arrays to zero 163 233 !!--------------------------------------------------------------------- 164 234 !! * Arguments … … 167 237 !! * Local variables 168 238 INTEGER :: jsec, &! loop on sections 169 iost, &! error for opening fileout 170 itotal ! nb_sec_max*nb_type_class*nb_class_max 239 itotal ! nb_sec_max*nb_type*nb_class_max 171 240 LOGICAL :: lldebug =.FALSE. ! debug a section 172 CHARACTER(len=160) :: clfileout ! fileout name173 241 174 242 175 INTEGER , DIMENSION(1) :: ish ! tmp array for mpp_sum176 INTEGER , DIMENSION(3) :: ish2 ! "177 REAL(wp), POINTER, DIMENSION(:) :: zwork ! "178 REAL(wp), POINTER, DIMENSION(:,:,:):: zsum ! "243 INTEGER , DIMENSION(1) :: ish ! tmp array for mpp_sum 244 INTEGER , DIMENSION(3) :: ish2 ! " 245 REAL(wp), POINTER, DIMENSION(:) :: zwork ! " 246 REAL(wp), POINTER, DIMENSION(:,:,:):: zsum ! " 179 247 180 248 !!--------------------------------------------------------------------- … … 182 250 183 251 IF( lk_mpp )THEN 184 itotal = nb_sec_max*nb_type _class*nb_class_max185 CALL wrk_alloc( itotal 186 CALL wrk_alloc( nb_sec_max,nb_type _class,nb_class_max , zsum )252 itotal = nb_sec_max*nb_type*nb_class_max 253 CALL wrk_alloc( itotal , zwork ) 254 CALL wrk_alloc( nb_sec_max,nb_type,nb_class_max , zsum ) 187 255 ENDIF 188 256 257 zwork(:) = 0.0 258 zsum(:,:,:) = 0.0 259 189 260 IF( lwp .AND. kt==nit000+nn_dct-1 ) THEN 190 261 WRITE(numout,*) " " … … 194 265 ENDIF 195 266 196 197 ! Compute transport and write only at nn_dctwri 198 IF( MOD(kt,nn_dct)==0 ) THEN 267 IF ( MOD(kt,nn_dct)==0 .or. & ! compute transport every nn_dct time steps 268 (ln_NOOS .and. kt==nn_it000 ) ) THEN ! also include first time step when calculating NOOS 25 hour averages 199 269 200 270 DO jsec=1,nb_sec 201 271 202 !debug this section computing ?203 272 lldebug=.FALSE. 204 273 IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. kt==nit000+nn_dct-1 .AND. lwp ) lldebug=.TRUE. 205 274 206 275 !Compute transport through section 207 CALL transport(secs(jsec),lldebug )276 CALL transport(secs(jsec),lldebug,jsec) 208 277 209 278 ENDDO … … 211 280 IF( MOD(kt,nn_dctwri)==0 )THEN 212 281 213 IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)" diadct: write at kt = ",kt282 IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)" diadct: average and write at kt = ",kt 214 283 284 !! divide arrays by nn_dctwri/nn_dct to obtain average 285 transports_3d(:,:,:,:)=transports_3d(:,:,:,:)/(nn_dctwri/nn_dct) 286 transports_2d(:,:,:) =transports_2d(:,:,:) /(nn_dctwri/nn_dct) 287 288 ! Sum over each class 289 DO jsec=1,nb_sec 290 CALL dia_dct_sum(secs(jsec),jsec) 291 ENDDO 292 215 293 !Sum on all procs 216 294 IF( lk_mpp )THEN 217 ish(1) = nb_sec_max*nb_type_class*nb_class_max 218 ish2 = (/nb_sec_max,nb_type_class,nb_class_max/) 295 zsum(:,:,:)=0.0_wp 296 ish(1) = nb_sec_max*nb_type*nb_class_max 297 ish2 = (/nb_sec_max,nb_type,nb_class_max/) 219 298 DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO 220 299 zwork(:)= RESHAPE(zsum(:,:,:), ish ) … … 227 306 DO jsec=1,nb_sec 228 307 229 IF( lwp )CALL dia_dct_wri(kt,jsec,secs(jsec)) 308 IF( lwp .and. .not. ln_NOOS )CALL dia_dct_wri(kt,jsec,secs(jsec)) 309 IF( lwp .and. ln_NOOS )CALL dia_dct_wri_NOOS(kt,jsec,secs(jsec)) ! use NOOS specific formatting 230 310 231 311 !nullify transports values after writing 312 transports_3d(:,jsec,:,:)=0.0 313 transports_2d(:,jsec,:)=0.0 232 314 secs(jsec)%transport(:,:)=0. 315 IF ( ln_NOOS ) CALL transport(secs(jsec),lldebug,jsec) ! reinitialise for next 25 hour instantaneous average (overlapping values) 233 316 234 317 ENDDO … … 238 321 ENDIF 239 322 323 IF ( MOD(kt,nn_dct_h)==0 ) THEN ! compute transport every nn_dct_h time steps 324 325 DO jsec=1,nb_sec 326 327 lldebug=.FALSE. 328 IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. kt==nit000+nn_dct_h-1 .AND. lwp ) lldebug=.TRUE. 329 330 !Compute transport through section 331 CALL transport_h(secs(jsec),lldebug,jsec) 332 333 ENDDO 334 335 IF( MOD(kt,nn_dctwri_h)==0 )THEN 336 337 IF( lwp .AND. kt==nit000+nn_dctwri_h-1 )WRITE(numout,*)" diadct: average and write hourly files at kt = ",kt 338 339 !! divide arrays by nn_dctwri/nn_dct to obtain average 340 transports_3d_h(:,:,:,:)=transports_3d_h(:,:,:,:)/(nn_dctwri_h/nn_dct_h) 341 transports_2d_h(:,:,:) =transports_2d_h(:,:,:) /(nn_dctwri_h/nn_dct_h) 342 343 ! Sum over each class 344 DO jsec=1,nb_sec 345 CALL dia_dct_sum_h(secs(jsec),jsec) 346 ENDDO 347 348 !Sum on all procs 349 IF( lk_mpp )THEN 350 ish(1) = nb_sec_max*nb_type*nb_class_max 351 ish2 = (/nb_sec_max,nb_type,nb_class_max/) 352 DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport_h(:,:) ; ENDDO 353 zwork(:)= RESHAPE(zsum(:,:,:), ish ) 354 CALL mpp_sum(zwork, ish(1)) 355 zsum(:,:,:)= RESHAPE(zwork,ish2) 356 DO jsec=1,nb_sec ; secs(jsec)%transport_h(:,:) = zsum(jsec,:,:) ; ENDDO 357 ENDIF 358 359 !Write the transport 360 DO jsec=1,nb_sec 361 362 IF( lwp .and. ln_NOOS )CALL dia_dct_wri_NOOS_h(kt/nn_dctwri_h,jsec,secs(jsec)) ! use NOOS specific formatting 363 364 !nullify transports values after writing 365 transports_3d_h(:,jsec,:,:)=0.0 366 transports_2d_h(:,jsec,:)=0.0 367 secs(jsec)%transport_h(:,:)=0. 368 IF ( ln_NOOS ) CALL transport_h(secs(jsec),lldebug,jsec) ! reinitialise for next 25 hour instantaneous average (overlapping values) 369 370 ENDDO 371 372 ENDIF 373 374 ENDIF 375 240 376 IF( lk_mpp )THEN 241 itotal = nb_sec_max*nb_type _class*nb_class_max242 CALL wrk_dealloc( itotal 243 CALL wrk_dealloc( nb_sec_max,nb_type _class,nb_class_max , zsum )377 itotal = nb_sec_max*nb_type*nb_class_max 378 CALL wrk_dealloc( itotal , zwork ) 379 CALL wrk_dealloc( nb_sec_max,nb_type,nb_class_max , zsum ) 244 380 ENDIF 245 381 … … 299 435 secs(jsec)%zlay = 99._wp 300 436 secs(jsec)%transport = 0._wp ; secs(jsec)%nb_point = 0 437 secs(jsec)%transport_h = 0._wp ; secs(jsec)%nb_point = 0 301 438 302 439 !read section's number / name / computing choices / classes / slopeSection / points number … … 331 468 332 469 WRITE(numout,*) " Section name : ",TRIM(secs(jsec)%name) 333 WRITE(numout,*) " Compute heat and salt transport? ",secs(jsec)%llstrpond470 WRITE(numout,*) " Compute temperature and salinity transports ? ",secs(jsec)%llstrpond 334 471 WRITE(numout,*) " Compute ice transport ? ",secs(jsec)%ll_ice_section 335 472 WRITE(numout,*) " Section crosses date-line ? ",secs(jsec)%ll_date_line … … 362 499 WRITE(numout,*)" List of points in global domain:" 363 500 DO jpt=1,iptglo 364 WRITE(numout,*)' # I J ',jpt,coordtemp(jpt) 501 WRITE(numout,*)' # I J ',jpt,coordtemp(jpt),directemp(jpt) 365 502 ENDDO 366 503 ENDIF … … 403 540 404 541 IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN 405 WRITE(narea+200,*)'avant secs(jsec)%nb_point iptloc ',secs(jsec)%nb_point,iptloc406 542 DO jpt = 1,iptloc 407 543 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1 408 544 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1 409 WRITE(narea+200,*)'avant # I J : ',iiglo,ijglo410 545 ENDDO 411 546 ENDIF … … 421 556 ENDIF 422 557 IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN 423 WRITE(narea+200,*)'apres secs(jsec)%nb_point iptloc ',secs(jsec)%nb_point,iptloc424 558 DO jpt = 1,secs(jsec)%nb_point 425 559 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1 426 560 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1 427 WRITE(narea+200,*)'apres # I J : ',iiglo,ijglo428 561 ENDDO 429 562 ENDIF … … 534 667 CALL wrk_dealloc( nb_point_max, idirec ) 535 668 CALL wrk_dealloc( 2, nb_point_max, icoord ) 669 536 670 END SUBROUTINE removepoints 537 671 538 SUBROUTINE transport(sec,ld_debug )672 SUBROUTINE transport(sec,ld_debug,jsec) 539 673 !!------------------------------------------------------------------------------------------- 540 674 !! *** ROUTINE transport *** 541 675 !! 542 !! ** Purpose : Compute the transport through a section 543 !! 544 !! ** Method :Transport through a given section is equal to the sum of transports 545 !! computed on each proc. 546 !! On each proc,transport is equal to the sum of transport computed through 547 !! segments linking each point of sec%listPoint with the next one. 548 !! 549 !! !BE carefull : 550 !! one section is a sum of segments 551 !! one segment is defined by 2 consectuve points in sec%listPoint 552 !! all points of sec%listPoint are positioned on the F-point of the cell. 676 !! Purpose :: Compute the transport for each point in a section 677 !! 678 !! Method :: Loop over each segment, and each vertical level and add the transport 679 !! Be aware : 680 !! One section is a sum of segments 681 !! One segment is defined by 2 consecutive points in sec%listPoint 682 !! All points of sec%listPoint are positioned on the F-point of the cell 553 683 !! 554 !! There are several loops: 555 !! loop on the density/temperature/salinity/level classes 684 !! There are two loops: 556 685 !! loop on the segment between 2 nodes 557 686 !! loop on the level jk 558 !! test on the density/temperature/salinity/level 559 !! 560 !! ** Output: sec%transport: volume/mass/ice/heat/salt transport in the 2 directions 561 !! 687 !! 688 !! ** Output: Arrays containing the volume,density,salinity,temperature etc 689 !! transports for each point in a section, summed over each nn_dct. 562 690 !! 563 691 !!------------------------------------------------------------------------------------------- … … 565 693 TYPE(SECTION),INTENT(INOUT) :: sec 566 694 LOGICAL ,INTENT(IN) :: ld_debug 695 INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section 567 696 568 697 !! * Local variables 569 INTEGER :: jk,jseg,jclass, &!loop on level/segment/classes 570 isgnu , isgnv ! 571 INTEGER :: ii, ij ! local integer 572 REAL(wp):: zumid , zvmid ,&!U/V velocity on a cell segment 573 zumid_ice , zvmid_ice ,&!U/V ice velocity 574 zTnorm ,&!transport of velocity through one cell's sides 575 ztransp1 , ztransp2 ,&!total transport in directions 1 and 2 576 ztemp1 , ztemp2 ,&!temperature transport " 577 zrhoi1 , zrhoi2 ,&!mass transport " 578 zrhop1 , zrhop2 ,&!mass transport " 579 zsal1 , zsal2 ,&!salinity transport " 580 zice_vol_pos , zice_vol_neg ,&!volume ice transport " 581 zice_surf_pos, zice_surf_neg !surface ice transport " 698 INTEGER :: jk,jseg,jclass, &!loop on level/segment/classes 699 isgnu , isgnv ! 700 REAL(wp):: zumid , zvmid , &!U/V velocity on a cell segment 701 zumid_ice , zvmid_ice , &!U/V ice velocity 702 zTnorm !transport of velocity through one cell's sides 582 703 REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 583 704 584 705 TYPE(POINT_SECTION) :: k 585 REAL(wp), POINTER, DIMENSION(:,:):: zsum ! 2D work array586 706 !!-------------------------------------------------------- 587 CALL wrk_alloc( nb_type_class , nb_class_max , zsum )588 707 589 708 IF( ld_debug )WRITE(numout,*)' Compute transport' 590 591 !----------------!592 ! INITIALIZATION !593 !----------------!594 zsum = 0._wp595 zice_surf_neg = 0._wp ; zice_surf_pos = 0._wp596 zice_vol_pos = 0._wp ; zice_vol_neg = 0._wp597 709 598 710 !---------------------------! … … 626 738 ELSE ; isgnv = 1 627 739 ENDIF 628 629 IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv 740 IF( sec%slopeSection .GE. 9999. ) isgnv = 1 741 742 IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv 630 743 631 744 !--------------------------------------! … … 670 783 END SELECT 671 784 672 !------------------------------- 673 ! LOOP ON THE DENSITY CLASSES | 674 !------------------------------- 675 !The computation is made for each density class 676 DO jclass=1,MAX(1,sec%nb_class-1) 677 678 ztransp1=0._wp ; zrhoi1=0._wp ; zrhop1=0._wp ; ztemp1=0._wp ;zsal1=0._wp 679 ztransp2=0._wp ; zrhoi2=0._wp ; zrhop2=0._wp ; ztemp2=0._wp ;zsal2=0._wp 680 681 !---------------------------| 682 ! LOOP ON THE LEVEL | 683 !---------------------------| 684 !Sum of the transport on the vertical 685 DO jk=1,jpk 686 687 688 ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point 689 SELECT CASE( sec%direction(jseg) ) 690 CASE(0,1) 691 ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 692 zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 693 zrhop = interp(k%I,k%J,jk,'V',rhop) 694 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 695 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) * vmask(k%I,k%J,1) 696 CASE(2,3) 697 ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 698 zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 699 zrhop = interp(k%I,k%J,jk,'U',rhop) 700 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 701 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1) 702 END SELECT 703 704 zfsdep= gdept(k%I,k%J,jk) 705 706 !----------------------------------------------! 707 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL! 708 !----------------------------------------------! 709 710 IF ( ( ((( zrhop .GE. (sec%zsigp(jclass)+1000. )) .AND. & 711 ( zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR. & 712 ( sec%zsigp(jclass) .EQ. 99.)) .AND. & 713 ((( zrhoi .GE. (sec%zsigi(jclass) + 1000. )) .AND. & 714 ( zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR. & 715 ( sec%zsigi(jclass) .EQ. 99.)) .AND. & 716 ((( zsn .GT. sec%zsal(jclass)) .AND. & 717 ( zsn .LE. sec%zsal(jclass+1))) .OR. & 718 ( sec%zsal(jclass) .EQ. 99.)) .AND. & 719 ((( ztn .GE. sec%ztem(jclass)) .AND. & 720 ( ztn .LE. sec%ztem(jclass+1))) .OR. & 721 ( sec%ztem(jclass) .EQ.99.)) .AND. & 722 ((( zfsdep .GE. sec%zlay(jclass)) .AND. & 723 ( zfsdep .LE. sec%zlay(jclass+1))) .OR. & 724 ( sec%zlay(jclass) .EQ. 99. )))) THEN 725 726 727 !compute velocity with the correct direction 728 SELECT CASE( sec%direction(jseg) ) 729 CASE(0,1) 730 zumid=0. 731 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 732 CASE(2,3) 733 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 734 zvmid=0. 735 END SELECT 736 737 !velocity* cell's length * cell's thickness 738 zTnorm=zumid*e2u(k%I,k%J)* fse3u(k%I,k%J,jk)+ & 739 zvmid*e1v(k%I,k%J)* fse3v(k%I,k%J,jk) 785 !---------------------------| 786 ! LOOP ON THE LEVEL | 787 !---------------------------| 788 !Sum of the transport on the vertical 789 DO jk=1,mbathy(k%I,k%J) 790 791 ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point 792 SELECT CASE( sec%direction(jseg) ) 793 CASE(0,1) 794 ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 795 zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 796 zrhop = interp(k%I,k%J,jk,'V',rhop) 797 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 798 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) * vmask(k%I,k%J,1) 799 CASE(2,3) 800 ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 801 zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 802 zrhop = interp(k%I,k%J,jk,'U',rhop) 803 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 804 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1) 805 END SELECT 806 807 zfsdep= gdept(k%I,k%J,jk) 808 809 !compute velocity with the correct direction 810 SELECT CASE( sec%direction(jseg) ) 811 CASE(0,1) 812 zumid=0. 813 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 814 CASE(2,3) 815 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 816 zvmid=0. 817 END SELECT 818 819 !zTnorm=transport through one cell; 820 !velocity* cell's length * cell's thickness 821 zTnorm=zumid*e2u(k%I,k%J)* fse3u(k%I,k%J,jk)+ & 822 zvmid*e1v(k%I,k%J)* fse3v(k%I,k%J,jk) 740 823 741 824 #if ! defined key_vvl 742 743 744 745 746 825 !add transport due to free surface 826 IF( jk==1 )THEN 827 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + & 828 zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 829 ENDIF 747 830 #endif 748 !COMPUTE TRANSPORT 749 !zTnorm=transport through one cell for one class 750 !ztransp1 or ztransp2=transport through one cell i 751 ! for one class for one direction 752 IF( zTnorm .GE. 0 )THEN 753 754 ztransp1=zTnorm+ztransp1 755 756 IF ( sec%llstrpond ) THEN 757 ztemp1 = ztemp1 + zTnorm * ztn 758 zsal1 = zsal1 + zTnorm * zsn 759 zrhoi1 = zrhoi1 + zTnorm * zrhoi 760 zrhop1 = zrhop1 + zTnorm * zrhop 761 ENDIF 762 763 ELSE 764 765 ztransp2=(zTnorm)+ztransp2 766 767 IF ( sec%llstrpond ) THEN 768 ztemp2 = ztemp2 + zTnorm * ztn 769 zsal2 = zsal2 + zTnorm * zsn 770 zrhoi2 = zrhoi2 + zTnorm * zrhoi 771 zrhop2 = zrhop2 + zTnorm * zrhop 772 ENDIF 773 ENDIF 774 775 776 ENDIF ! end of density test 777 ENDDO!end of loop on the level 778 779 !ZSUM=TRANSPORT FOR EACH CLASSES FOR THE DIRECTIONS 780 !--------------------------------------------------- 781 zsum(1,jclass) = zsum(1,jclass)+ztransp1 782 zsum(2,jclass) = zsum(2,jclass)+ztransp2 783 IF( sec%llstrpond )THEN 784 zsum(3 ,jclass) = zsum( 3,jclass)+zrhoi1 785 zsum(4 ,jclass) = zsum( 4,jclass)+zrhoi2 786 zsum(5 ,jclass) = zsum( 5,jclass)+zrhop1 787 zsum(6 ,jclass) = zsum( 6,jclass)+zrhop2 788 zsum(7 ,jclass) = zsum( 7,jclass)+ztemp1 789 zsum(8 ,jclass) = zsum( 8,jclass)+ztemp2 790 zsum(9 ,jclass) = zsum( 9,jclass)+zsal1 791 zsum(10,jclass) = zsum(10,jclass)+zsal2 831 !COMPUTE TRANSPORT 832 833 transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm 834 835 IF ( sec%llstrpond ) THEN 836 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk) + zTnorm * zrhoi 837 transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk) + zTnorm * zrhop 838 transports_3d(4,jsec,jseg,jk) = transports_3d(4,jsec,jseg,jk) + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp 839 transports_3d(5,jsec,jseg,jk) = transports_3d(5,jsec,jseg,jk) + zTnorm * 0.001 * zsn * 1026._wp 792 840 ENDIF 793 794 ENDDO !end of loop on the density classes841 842 ENDDO !end of loop on the level 795 843 796 844 #if defined key_lim2 || defined key_lim3 … … 816 864 zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J) 817 865 818 IF( zTnorm .GE. 0)THEN 819 zice_vol_pos = (zTnorm)* & 820 (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 821 *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) + & 822 hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 823 +zice_vol_pos 824 zice_surf_pos = (zTnorm)* & 825 (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 826 +zice_surf_pos 827 ELSE 828 zice_vol_neg=(zTnorm)* & 866 transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)* & 829 867 (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 830 868 *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) + & 831 869 hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 832 +zice_vol_neg 833 zice_surf_neg=(zTnorm)* & 834 (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 835 +zice_surf_neg 836 ENDIF 837 838 zsum(11,1) = zsum(11,1)+zice_vol_pos 839 zsum(12,1) = zsum(12,1)+zice_vol_neg 840 zsum(13,1) = zsum(13,1)+zice_surf_pos 841 zsum(14,1) = zsum(14,1)+zice_surf_neg 870 +zice_vol_pos 871 transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)* & 872 (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 873 +zice_surf_pos 842 874 843 875 ENDIF !end of ice case … … 846 878 ENDDO !end of loop on the segment 847 879 848 849 ELSE !if sec%nb_point =0850 zsum(1:2,:)=0.851 IF (sec%llstrpond) zsum(3:10,:)=0.852 zsum( 11:14,:)=0.853 880 ENDIF !end of sec%nb_point =0 case 854 855 !-------------------------------|856 !FINISH COMPUTING TRANSPORTS |857 !-------------------------------|858 DO jclass=1,MAX(1,sec%nb_class-1)859 sec%transport(1,jclass)=sec%transport(1,jclass)+zsum(1,jclass)*1.E-6860 sec%transport(2,jclass)=sec%transport(2,jclass)+zsum(2,jclass)*1.E-6861 IF( sec%llstrpond ) THEN862 IF( zsum(1,jclass) .NE. 0._wp ) THEN863 sec%transport( 3,jclass) = sec%transport( 3,jclass) + zsum( 3,jclass)/zsum(1,jclass)864 sec%transport( 5,jclass) = sec%transport( 5,jclass) + zsum( 5,jclass)/zsum(1,jclass)865 sec%transport( 7,jclass) = sec%transport( 7,jclass) + zsum( 7,jclass)866 sec%transport( 9,jclass) = sec%transport( 9,jclass) + zsum( 9,jclass)867 ENDIF868 IF( zsum(2,jclass) .NE. 0._wp )THEN869 sec%transport( 4,jclass) = sec%transport( 4,jclass) + zsum( 4,jclass)/zsum(2,jclass)870 sec%transport( 6,jclass) = sec%transport( 6,jclass) + zsum( 6,jclass)/zsum(2,jclass)871 sec%transport( 8,jclass) = sec%transport( 8,jclass) + zsum( 8,jclass)872 sec%transport(10,jclass) = sec%transport(10,jclass) + zsum(10,jclass)873 ENDIF874 ELSE875 sec%transport( 3,jclass) = 0._wp876 sec%transport( 4,jclass) = 0._wp877 sec%transport( 5,jclass) = 0._wp878 sec%transport( 6,jclass) = 0._wp879 sec%transport( 7,jclass) = 0._wp880 sec%transport( 8,jclass) = 0._wp881 sec%transport(10,jclass) = 0._wp882 ENDIF883 ENDDO884 885 IF( sec%ll_ice_section ) THEN886 sec%transport( 9,1)=sec%transport( 9,1)+zsum( 9,1)*1.E-6887 sec%transport(10,1)=sec%transport(10,1)+zsum(10,1)*1.E-6888 sec%transport(11,1)=sec%transport(11,1)+zsum(11,1)*1.E-6889 sec%transport(12,1)=sec%transport(12,1)+zsum(12,1)*1.E-6890 ENDIF891 892 CALL wrk_dealloc( nb_type_class , nb_class_max , zsum )893 881 ! 894 882 END SUBROUTINE transport 895 883 884 SUBROUTINE transport_h(sec,ld_debug,jsec) 885 !!------------------------------------------------------------------------------------------- 886 !! *** ROUTINE hourly_transport *** 887 !! 888 !! Exactly as routine transport but for data summed at 889 !! each time step and averaged each hour 890 !! 891 !! Purpose :: Compute the transport for each point in a section 892 !! 893 !! Method :: Loop over each segment, and each vertical level and add the transport 894 !! Be aware : 895 !! One section is a sum of segments 896 !! One segment is defined by 2 consecutive points in sec%listPoint 897 !! All points of sec%listPoint are positioned on the F-point of the cell 898 !! 899 !! There are two loops: 900 !! loop on the segment between 2 nodes 901 !! loop on the level jk 902 !! 903 !! ** Output: Arrays containing the volume,density,salinity,temperature etc 904 !! transports for each point in a section, summed over each nn_dct. 905 !! 906 !!------------------------------------------------------------------------------------------- 907 !! * Arguments 908 TYPE(SECTION),INTENT(INOUT) :: sec 909 LOGICAL ,INTENT(IN) :: ld_debug 910 INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section 911 912 !! * Local variables 913 INTEGER :: jk,jseg,jclass, &!loop on level/segment/classes 914 isgnu , isgnv ! 915 REAL(wp):: zumid , zvmid , &!U/V velocity on a cell segment 916 zumid_ice , zvmid_ice , &!U/V ice velocity 917 zTnorm !transport of velocity through one cell's sides 918 REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 919 920 TYPE(POINT_SECTION) :: k 921 !!-------------------------------------------------------- 922 923 IF( ld_debug )WRITE(numout,*)' Compute transport' 924 925 !---------------------------! 926 ! COMPUTE TRANSPORT ! 927 !---------------------------! 928 IF(sec%nb_point .NE. 0)THEN 929 930 !---------------------------------------------------------------------------------------------------- 931 !Compute sign for velocities: 932 ! 933 !convention: 934 ! non horizontal section: direction + is toward left hand of section 935 ! horizontal section: direction + is toward north of section 936 ! 937 ! 938 ! slopeSection < 0 slopeSection > 0 slopeSection=inf slopeSection=0 939 ! ---------------- ----------------- --------------- -------------- 940 ! 941 ! isgnv=1 direction + 942 ! ______ _____ ______ 943 ! | //| | | direction + 944 ! | isgnu=1 // | |isgnu=1 |isgnu=1 /|\ 945 ! |_______ // ______| \\ | ---\ | 946 ! | | isgnv=-1 \\ | | ---/ direction + ____________ 947 ! | | __\\| | 948 ! | | direction + | isgnv=1 949 ! 950 !---------------------------------------------------------------------------------------------------- 951 isgnu = 1 952 IF( sec%slopeSection .GT. 0 ) THEN ; isgnv = -1 953 ELSE ; isgnv = 1 954 ENDIF 955 956 IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv 957 958 !--------------------------------------! 959 ! LOOP ON THE SEGMENT BETWEEN 2 NODES ! 960 !--------------------------------------! 961 DO jseg=1,MAX(sec%nb_point-1,0) 962 963 !------------------------------------------------------------------------------------------- 964 ! Select the appropriate coordinate for computing the velocity of the segment 965 ! 966 ! CASE(0) Case (2) 967 ! ------- -------- 968 ! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j) 969 ! F(i,j)----------V(i+1,j)-------F(i+1,j) | 970 ! | 971 ! | 972 ! | 973 ! Case (3) U(i,j) 974 ! -------- | 975 ! | 976 ! listPoint(jseg+1) F(i,j+1) | 977 ! | | 978 ! | | 979 ! | listPoint(jseg+1) F(i,j-1) 980 ! | 981 ! | 982 ! U(i,j+1) 983 ! | Case(1) 984 ! | ------ 985 ! | 986 ! | listPoint(jseg+1) listPoint(jseg) 987 ! | F(i-1,j)-----------V(i,j) -------f(jseg) 988 ! listPoint(jseg) F(i,j) 989 ! 990 !------------------------------------------------------------------------------------------- 991 992 SELECT CASE( sec%direction(jseg) ) 993 CASE(0) ; k = sec%listPoint(jseg) 994 CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 995 CASE(2) ; k = sec%listPoint(jseg) 996 CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 997 END SELECT 998 999 !---------------------------| 1000 ! LOOP ON THE LEVEL | 1001 !---------------------------| 1002 !Sum of the transport on the vertical 1003 DO jk=1,mbathy(k%I,k%J) 1004 1005 ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point 1006 SELECT CASE( sec%direction(jseg) ) 1007 CASE(0,1) 1008 ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 1009 zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 1010 zrhop = interp(k%I,k%J,jk,'V',rhop) 1011 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 1012 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) * vmask(k%I,k%J,1) 1013 CASE(2,3) 1014 ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 1015 zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 1016 zrhop = interp(k%I,k%J,jk,'U',rhop) 1017 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 1018 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1) 1019 END SELECT 1020 1021 zfsdep= gdept(k%I,k%J,jk) 1022 1023 !compute velocity with the correct direction 1024 SELECT CASE( sec%direction(jseg) ) 1025 CASE(0,1) 1026 zumid=0. 1027 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 1028 CASE(2,3) 1029 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 1030 zvmid=0. 1031 END SELECT 1032 1033 !zTnorm=transport through one cell; 1034 !velocity* cell's length * cell's thickness 1035 zTnorm=zumid*e2u(k%I,k%J)* fse3u(k%I,k%J,jk)+ & 1036 zvmid*e1v(k%I,k%J)* fse3v(k%I,k%J,jk) 1037 1038 #if ! defined key_vvl 1039 !add transport due to free surface 1040 IF( jk==1 )THEN 1041 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + & 1042 zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk) 1043 ENDIF 1044 #endif 1045 !COMPUTE TRANSPORT 1046 1047 transports_3d_h(1,jsec,jseg,jk) = transports_3d_h(1,jsec,jseg,jk) + zTnorm 1048 1049 IF ( sec%llstrpond ) THEN 1050 transports_3d_h(2,jsec,jseg,jk) = transports_3d_h(2,jsec,jseg,jk) + zTnorm * zrhoi 1051 transports_3d_h(3,jsec,jseg,jk) = transports_3d_h(3,jsec,jseg,jk) + zTnorm * zrhop 1052 transports_3d_h(4,jsec,jseg,jk) = transports_3d_h(4,jsec,jseg,jk) + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp 1053 transports_3d_h(5,jsec,jseg,jk) = transports_3d_h(5,jsec,jseg,jk) + zTnorm * 0.001 * zsn * 1026._wp 1054 ENDIF 1055 1056 ENDDO !end of loop on the level 1057 1058 #if defined key_lim2 || defined key_lim3 1059 1060 !ICE CASE 1061 !------------ 1062 IF( sec%ll_ice_section )THEN 1063 SELECT CASE (sec%direction(jseg)) 1064 CASE(0) 1065 zumid_ice = 0 1066 zvmid_ice = isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1)) 1067 CASE(1) 1068 zumid_ice = 0 1069 zvmid_ice = isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1)) 1070 CASE(2) 1071 zvmid_ice = 0 1072 zumid_ice = isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1)) 1073 CASE(3) 1074 zvmid_ice = 0 1075 zumid_ice = isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1)) 1076 END SELECT 1077 1078 zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J) 1079 1080 transports_2d_h(1,jsec,jseg) = transports_2d_h(1,jsec,jseg) + (zTnorm)* & 1081 (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 1082 *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) + & 1083 hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 1084 +zice_vol_pos 1085 transports_2d_h(2,jsec,jseg) = transports_2d_h(2,jsec,jseg) + (zTnorm)* & 1086 (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) & 1087 +zice_surf_pos 1088 1089 ENDIF !end of ice case 1090 #endif 1091 1092 ENDDO !end of loop on the segment 1093 1094 ENDIF !end of sec%nb_point =0 case 1095 ! 1096 END SUBROUTINE transport_h 1097 1098 SUBROUTINE dia_dct_sum(sec,jsec) 1099 !!------------------------------------------------------------- 1100 !! Purpose: Average the transport over nn_dctwri time steps 1101 !! and sum over the density/salinity/temperature/depth classes 1102 !! 1103 !! Method: 1104 !! Sum over relevant grid cells to obtain values 1105 !! for each 1106 !! There are several loops: 1107 !! loop on the segment between 2 nodes 1108 !! loop on the level jk 1109 !! loop on the density/temperature/salinity/level classes 1110 !! test on the density/temperature/salinity/level 1111 !! 1112 !! ** Method :Transport through a given section is equal to the sum of transports 1113 !! computed on each proc. 1114 !! On each proc,transport is equal to the sum of transport computed through 1115 !! segments linking each point of sec%listPoint with the next one. 1116 !! 1117 !!------------------------------------------------------------- 1118 !! * arguments 1119 TYPE(SECTION),INTENT(INOUT) :: sec 1120 INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section 1121 1122 TYPE(POINT_SECTION) :: k 1123 INTEGER :: jk,jseg,jclass !loop on level/segment/classes 1124 REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 1125 !!------------------------------------------------------------- 1126 1127 !! Sum the relevant segments to obtain values for each class 1128 IF(sec%nb_point .NE. 0)THEN 1129 1130 !--------------------------------------! 1131 ! LOOP ON THE SEGMENT BETWEEN 2 NODES ! 1132 !--------------------------------------! 1133 DO jseg=1,MAX(sec%nb_point-1,0) 1134 1135 !------------------------------------------------------------------------------------------- 1136 ! Select the appropriate coordinate for computing the velocity of the segment 1137 ! 1138 ! CASE(0) Case (2) 1139 ! ------- -------- 1140 ! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j) 1141 ! F(i,j)----------V(i+1,j)-------F(i+1,j) | 1142 ! | 1143 ! | 1144 ! | 1145 ! Case (3) U(i,j) 1146 ! -------- | 1147 ! | 1148 ! listPoint(jseg+1) F(i,j+1) | 1149 ! | | 1150 ! | | 1151 ! | listPoint(jseg+1) F(i,j-1) 1152 ! | 1153 ! | 1154 ! U(i,j+1) 1155 ! | Case(1) 1156 ! | ------ 1157 ! | 1158 ! | listPoint(jseg+1) listPoint(jseg) 1159 ! | F(i-1,j)-----------V(i,j) -------f(jseg) 1160 ! listPoint(jseg) F(i,j) 1161 ! 1162 !------------------------------------------------------------------------------------------- 1163 1164 SELECT CASE( sec%direction(jseg) ) 1165 CASE(0) ; k = sec%listPoint(jseg) 1166 CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 1167 CASE(2) ; k = sec%listPoint(jseg) 1168 CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 1169 END SELECT 1170 1171 !---------------------------| 1172 ! LOOP ON THE LEVEL | 1173 !---------------------------| 1174 !Sum of the transport on the vertical 1175 DO jk=1,mbathy(k%I,k%J) 1176 1177 ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point 1178 SELECT CASE( sec%direction(jseg) ) 1179 CASE(0,1) 1180 ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 1181 zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 1182 zrhop = interp(k%I,k%J,jk,'V',rhop) 1183 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 1184 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) * vmask(k%I,k%J,1) 1185 CASE(2,3) 1186 ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 1187 zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 1188 zrhop = interp(k%I,k%J,jk,'U',rhop) 1189 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 1190 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1) 1191 END SELECT 1192 1193 zfsdep= gdept(k%I,k%J,jk) 1194 1195 !------------------------------- 1196 ! LOOP ON THE DENSITY CLASSES | 1197 !------------------------------- 1198 !The computation is made for each density/heat/salt/... class 1199 DO jclass=1,MAX(1,sec%nb_class-1) 1200 1201 !----------------------------------------------! 1202 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL! 1203 !----------------------------------------------! 1204 1205 IF ( ( & 1206 ((( zrhop .GE. (sec%zsigp(jclass)+1000. )) .AND. & 1207 ( zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR. & 1208 ( sec%zsigp(jclass) .EQ. 99.)) .AND. & 1209 1210 ((( zrhoi .GE. (sec%zsigi(jclass) + 1000. )) .AND. & 1211 ( zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR. & 1212 ( sec%zsigi(jclass) .EQ. 99.)) .AND. & 1213 1214 ((( zsn .GT. sec%zsal(jclass)) .AND. & 1215 ( zsn .LE. sec%zsal(jclass+1))) .OR. & 1216 ( sec%zsal(jclass) .EQ. 99.)) .AND. & 1217 1218 ((( ztn .GE. sec%ztem(jclass)) .AND. & 1219 ( ztn .LE. sec%ztem(jclass+1))) .OR. & 1220 ( sec%ztem(jclass) .EQ.99.)) .AND. & 1221 1222 ((( zfsdep .GE. sec%zlay(jclass)) .AND. & 1223 ( zfsdep .LE. sec%zlay(jclass+1))) .OR. & 1224 ( sec%zlay(jclass) .EQ. 99. )) & 1225 )) THEN 1226 1227 !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS 1228 !---------------------------------------------------------------------------- 1229 IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN 1230 sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk) 1231 ELSE 1232 sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk) 1233 ENDIF 1234 IF( sec%llstrpond )THEN 1235 1236 IF( transports_3d(1,jsec,jseg,jk) .NE. 0._wp ) THEN 1237 1238 IF (transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN 1239 sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) 1240 ELSE 1241 sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) 1242 ENDIF 1243 1244 IF ( transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN 1245 sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) 1246 ELSE 1247 sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) 1248 ENDIF 1249 1250 ENDIF 1251 1252 IF ( transports_3d(4,jsec,jseg,jk) .GE. 0.0 ) THEN 1253 sec%transport(7,jclass) = sec%transport(7,jclass)+transports_3d(4,jsec,jseg,jk) 1254 ELSE 1255 sec%transport(8,jclass) = sec%transport(8,jclass)+transports_3d(4,jsec,jseg,jk) 1256 ENDIF 1257 1258 IF ( transports_3d(5,jsec,jseg,jk) .GE. 0.0 ) THEN 1259 sec%transport( 9,jclass) = sec%transport( 9,jclass)+transports_3d(5,jsec,jseg,jk) 1260 ELSE 1261 sec%transport(10,jclass) = sec%transport(10,jclass)+transports_3d(5,jsec,jseg,jk) 1262 ENDIF 1263 1264 ELSE 1265 sec%transport( 3,jclass) = 0._wp 1266 sec%transport( 4,jclass) = 0._wp 1267 sec%transport( 5,jclass) = 0._wp 1268 sec%transport( 6,jclass) = 0._wp 1269 sec%transport( 7,jclass) = 0._wp 1270 sec%transport( 8,jclass) = 0._wp 1271 sec%transport( 9,jclass) = 0._wp 1272 sec%transport(10,jclass) = 0._wp 1273 ENDIF 1274 1275 ENDIF ! end of test if point is in class 1276 1277 ENDDO ! end of loop on the classes 1278 1279 ENDDO ! loop over jk 1280 1281 #if defined key_lim2 || defined key_lim3 1282 1283 !ICE CASE 1284 IF( sec%ll_ice_section )THEN 1285 1286 IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN 1287 sec%transport(11,1) = sec%transport(11,1)+transports_2d(1,jsec,jseg) 1288 ELSE 1289 sec%transport(12,1) = sec%transport(12,1)+transports_2d(1,jsec,jseg) 1290 ENDIF 1291 1292 IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN 1293 sec%transport(13,1) = sec%transport(13,1)+transports_2d(2,jsec,jseg) 1294 ELSE 1295 sec%transport(14,1) = sec%transport(14,1)+transports_2d(2,jsec,jseg) 1296 ENDIF 1297 1298 ENDIF !end of ice case 1299 #endif 1300 1301 ENDDO !end of loop on the segment 1302 1303 ELSE !if sec%nb_point =0 1304 sec%transport(1:2,:)=0. 1305 IF (sec%llstrpond) sec%transport(3:10,:)=0. 1306 IF (sec%ll_ice_section) sec%transport( 11:14,:)=0. 1307 ENDIF !end of sec%nb_point =0 case 1308 1309 END SUBROUTINE dia_dct_sum 1310 1311 SUBROUTINE dia_dct_sum_h(sec,jsec) 1312 !!------------------------------------------------------------- 1313 !! Exactly as dia_dct_sum but for hourly files containing data summed at each time step 1314 !! 1315 !! Purpose: Average the transport over nn_dctwri time steps 1316 !! and sum over the density/salinity/temperature/depth classes 1317 !! 1318 !! Method: 1319 !! Sum over relevant grid cells to obtain values 1320 !! for each 1321 !! There are several loops: 1322 !! loop on the segment between 2 nodes 1323 !! loop on the level jk 1324 !! loop on the density/temperature/salinity/level classes 1325 !! test on the density/temperature/salinity/level 1326 !! 1327 !! ** Method :Transport through a given section is equal to the sum of transports 1328 !! computed on each proc. 1329 !! On each proc,transport is equal to the sum of transport computed through 1330 !! segments linking each point of sec%listPoint with the next one. 1331 !! 1332 !!------------------------------------------------------------- 1333 !! * arguments 1334 TYPE(SECTION),INTENT(INOUT) :: sec 1335 INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section 1336 1337 TYPE(POINT_SECTION) :: k 1338 INTEGER :: jk,jseg,jclass !loop on level/segment/classes 1339 REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 1340 !!------------------------------------------------------------- 1341 1342 !! Sum the relevant segments to obtain values for each class 1343 IF(sec%nb_point .NE. 0)THEN 1344 1345 !--------------------------------------! 1346 ! LOOP ON THE SEGMENT BETWEEN 2 NODES ! 1347 !--------------------------------------! 1348 DO jseg=1,MAX(sec%nb_point-1,0) 1349 1350 !------------------------------------------------------------------------------------------- 1351 ! Select the appropriate coordinate for computing the velocity of the segment 1352 ! 1353 ! CASE(0) Case (2) 1354 ! ------- -------- 1355 ! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j) 1356 ! F(i,j)----------V(i+1,j)-------F(i+1,j) | 1357 ! | 1358 ! | 1359 ! | 1360 ! Case (3) U(i,j) 1361 ! -------- | 1362 ! | 1363 ! listPoint(jseg+1) F(i,j+1) | 1364 ! | | 1365 ! | | 1366 ! | listPoint(jseg+1) F(i,j-1) 1367 ! | 1368 ! | 1369 ! U(i,j+1) 1370 ! | Case(1) 1371 ! | ------ 1372 ! | 1373 ! | listPoint(jseg+1) listPoint(jseg) 1374 ! | F(i-1,j)-----------V(i,j) -------f(jseg) 1375 ! listPoint(jseg) F(i,j) 1376 ! 1377 !------------------------------------------------------------------------------------------- 1378 1379 SELECT CASE( sec%direction(jseg) ) 1380 CASE(0) ; k = sec%listPoint(jseg) 1381 CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J) 1382 CASE(2) ; k = sec%listPoint(jseg) 1383 CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1) 1384 END SELECT 1385 1386 !---------------------------| 1387 ! LOOP ON THE LEVEL | 1388 !---------------------------| 1389 !Sum of the transport on the vertical 1390 DO jk=1,mbathy(k%I,k%J) 1391 1392 ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point 1393 SELECT CASE( sec%direction(jseg) ) 1394 CASE(0,1) 1395 ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 1396 zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 1397 zrhop = interp(k%I,k%J,jk,'V',rhop) 1398 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 1399 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) * vmask(k%I,k%J,1) 1400 CASE(2,3) 1401 ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 1402 zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 1403 zrhop = interp(k%I,k%J,jk,'U',rhop) 1404 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 1405 zsshn = 0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) * umask(k%I,k%J,1) 1406 END SELECT 1407 1408 zfsdep= gdept(k%I,k%J,jk) 1409 1410 !------------------------------- 1411 ! LOOP ON THE DENSITY CLASSES | 1412 !------------------------------- 1413 !The computation is made for each density/heat/salt/... class 1414 DO jclass=1,MAX(1,sec%nb_class-1) 1415 1416 !----------------------------------------------! 1417 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL! 1418 !----------------------------------------------! 1419 1420 IF ( ( & 1421 ((( zrhop .GE. (sec%zsigp(jclass)+1000. )) .AND. & 1422 ( zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR. & 1423 ( sec%zsigp(jclass) .EQ. 99.)) .AND. & 1424 1425 ((( zrhoi .GE. (sec%zsigi(jclass) + 1000. )) .AND. & 1426 ( zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR. & 1427 ( sec%zsigi(jclass) .EQ. 99.)) .AND. & 1428 1429 ((( zsn .GT. sec%zsal(jclass)) .AND. & 1430 ( zsn .LE. sec%zsal(jclass+1))) .OR. & 1431 ( sec%zsal(jclass) .EQ. 99.)) .AND. & 1432 1433 ((( ztn .GE. sec%ztem(jclass)) .AND. & 1434 ( ztn .LE. sec%ztem(jclass+1))) .OR. & 1435 ( sec%ztem(jclass) .EQ.99.)) .AND. & 1436 1437 ((( zfsdep .GE. sec%zlay(jclass)) .AND. & 1438 ( zfsdep .LE. sec%zlay(jclass+1))) .OR. & 1439 ( sec%zlay(jclass) .EQ. 99. )) & 1440 )) THEN 1441 1442 !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS 1443 !---------------------------------------------------------------------------- 1444 IF (transports_3d_h(1,jsec,jseg,jk) .GE. 0.0) THEN 1445 sec%transport_h(1,jclass) = sec%transport_h(1,jclass)+transports_3d_h(1,jsec,jseg,jk) 1446 ELSE 1447 sec%transport_h(2,jclass) = sec%transport_h(2,jclass)+transports_3d_h(1,jsec,jseg,jk) 1448 ENDIF 1449 IF( sec%llstrpond )THEN 1450 1451 IF( transports_3d_h(1,jsec,jseg,jk) .NE. 0._wp ) THEN 1452 1453 IF (transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) .GE. 0.0 ) THEN 1454 sec%transport_h(3,jclass) = sec%transport_h(3,jclass)+transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) 1455 ELSE 1456 sec%transport_h(4,jclass) = sec%transport_h(4,jclass)+transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) 1457 ENDIF 1458 1459 IF ( transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) .GE. 0.0 ) THEN 1460 sec%transport_h(5,jclass) = sec%transport_h(5,jclass)+transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) 1461 ELSE 1462 sec%transport_h(6,jclass) = sec%transport_h(6,jclass)+transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) 1463 ENDIF 1464 1465 ENDIF 1466 1467 IF ( transports_3d_h(4,jsec,jseg,jk) .GE. 0.0 ) THEN 1468 sec%transport_h(7,jclass) = sec%transport_h(7,jclass)+transports_3d_h(4,jsec,jseg,jk) 1469 ELSE 1470 sec%transport_h(8,jclass) = sec%transport_h(8,jclass)+transports_3d_h(4,jsec,jseg,jk) 1471 ENDIF 1472 1473 IF ( transports_3d_h(5,jsec,jseg,jk) .GE. 0.0 ) THEN 1474 sec%transport_h( 9,jclass) = sec%transport_h( 9,jclass)+transports_3d_h(5,jsec,jseg,jk) 1475 ELSE 1476 sec%transport_h(10,jclass) = sec%transport_h(10,jclass)+transports_3d_h(5,jsec,jseg,jk) 1477 ENDIF 1478 1479 ELSE 1480 sec%transport_h( 3,jclass) = 0._wp 1481 sec%transport_h( 4,jclass) = 0._wp 1482 sec%transport_h( 5,jclass) = 0._wp 1483 sec%transport_h( 6,jclass) = 0._wp 1484 sec%transport_h( 7,jclass) = 0._wp 1485 sec%transport_h( 8,jclass) = 0._wp 1486 sec%transport_h( 9,jclass) = 0._wp 1487 sec%transport_h(10,jclass) = 0._wp 1488 ENDIF 1489 1490 ENDIF ! end of test if point is in class 1491 1492 ENDDO ! end of loop on the classes 1493 1494 ENDDO ! loop over jk 1495 1496 #if defined key_lim2 || defined key_lim3 1497 1498 !ICE CASE 1499 IF( sec%ll_ice_section )THEN 1500 1501 IF ( transports_2d_h(1,jsec,jseg) .GE. 0.0 ) THEN 1502 sec%transport_h(11,1) = sec%transport_h(11,1)+transports_2d_h(1,jsec,jseg) 1503 ELSE 1504 sec%transport_h(12,1) = sec%transport_h(12,1)+transports_2d_h(1,jsec,jseg) 1505 ENDIF 1506 1507 IF ( transports_2d_h(3,jsec,jseg) .GE. 0.0 ) THEN 1508 sec%transport_h(13,1) = sec%transport_h(13,1)+transports_2d_h(2,jsec,jseg) 1509 ELSE 1510 sec%transport_h(14,1) = sec%transport_h(14,1)+transports_2d_h(2,jsec,jseg) 1511 ENDIF 1512 1513 ENDIF !end of ice case 1514 #endif 1515 1516 ENDDO !end of loop on the segment 1517 1518 ELSE !if sec%nb_point =0 1519 sec%transport_h(1:2,:)=0. 1520 IF (sec%llstrpond) sec%transport_h(3:10,:)=0. 1521 IF (sec%ll_ice_section) sec%transport_h( 11:14,:)=0. 1522 ENDIF !end of sec%nb_point =0 case 1523 1524 END SUBROUTINE dia_dct_sum_h 1525 1526 SUBROUTINE dia_dct_wri_NOOS(kt,ksec,sec) 1527 !!------------------------------------------------------------- 1528 !! Write transport output in numdct using NOOS formatting 1529 !! 1530 !! Purpose: Write transports in ascii files 1531 !! 1532 !! Method: 1533 !! 1. Write volume transports in "volume_transport" 1534 !! Unit: Sv : area * Velocity / 1.e6 1535 !! 1536 !! 2. Write heat transports in "heat_transport" 1537 !! Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15 1538 !! 1539 !! 3. Write salt transports in "salt_transport" 1540 !! Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6 1541 !! 1542 !!------------------------------------------------------------- 1543 !!arguments 1544 INTEGER, INTENT(IN) :: kt ! time-step 1545 TYPE(SECTION), INTENT(INOUT) :: sec ! section to write 1546 INTEGER ,INTENT(IN) :: ksec ! section number 1547 1548 !!local declarations 1549 INTEGER :: jclass,ji ! Dummy loop 1550 CHARACTER(len=2) :: classe ! Classname 1551 REAL(wp) :: zbnd1,zbnd2 ! Class bounds 1552 REAL(wp) :: zslope ! section's slope coeff 1553 ! 1554 REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace 1555 !!------------------------------------------------------------- 1556 CALL wrk_alloc(nb_type , zsumclasses ) 1557 1558 zsumclasses(:)=0._wp 1559 zslope = sec%slopeSection 1560 1561 WRITE(numdct_NOOS,'(I4,a1,I2,a1,I2,a12,i3,a17,i3,a10,a25)') nyear,'.',nmonth,'.',nday,' Transect:',ksec-1,' No. of layers:',sec%nb_class-1,' Name: ',sec%name 1562 1563 DO jclass=1,MAX(1,sec%nb_class-1) 1564 zsumclasses(1:nb_type)=zsumclasses(1:nb_type)+sec%transport(1:nb_type,jclass) 1565 ENDDO 1566 1567 classe = 'total ' 1568 zbnd1 = 0._wp 1569 zbnd2 = 0._wp 1570 1571 IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 1572 WRITE(numdct_NOOS,'(9e12.4E2)') -(zsumclasses( 1)+zsumclasses( 2)), -zsumclasses( 2),-zsumclasses( 1), & 1573 -(zsumclasses( 7)+zsumclasses( 8)), -zsumclasses( 8),-zsumclasses( 7), & 1574 -(zsumclasses( 9)+zsumclasses(10)), -zsumclasses(10),-zsumclasses( 9) 1575 ELSE 1576 WRITE(numdct_NOOS,'(9e12.4E2)') zsumclasses( 1)+zsumclasses( 2) , zsumclasses( 1), zsumclasses( 2), & 1577 zsumclasses( 7)+zsumclasses( 8) , zsumclasses( 7), zsumclasses( 8), & 1578 zsumclasses( 9)+zsumclasses(10) , zsumclasses( 9), zsumclasses(10) 1579 ENDIF 1580 1581 DO jclass=1,MAX(1,sec%nb_class-1) 1582 1583 classe = 'N ' 1584 zbnd1 = 0._wp 1585 zbnd2 = 0._wp 1586 1587 !insitu density classes transports 1588 IF( ( sec%zsigi(jclass) .NE. 99._wp ) .AND. & 1589 ( sec%zsigi(jclass+1) .NE. 99._wp ) )THEN 1590 classe = 'DI ' 1591 zbnd1 = sec%zsigi(jclass) 1592 zbnd2 = sec%zsigi(jclass+1) 1593 ENDIF 1594 !potential density classes transports 1595 IF( ( sec%zsigp(jclass) .NE. 99._wp ) .AND. & 1596 ( sec%zsigp(jclass+1) .NE. 99._wp ) )THEN 1597 classe = 'DP ' 1598 zbnd1 = sec%zsigp(jclass) 1599 zbnd2 = sec%zsigp(jclass+1) 1600 ENDIF 1601 !depth classes transports 1602 IF( ( sec%zlay(jclass) .NE. 99._wp ) .AND. & 1603 ( sec%zlay(jclass+1) .NE. 99._wp ) )THEN 1604 classe = 'Z ' 1605 zbnd1 = sec%zlay(jclass) 1606 zbnd2 = sec%zlay(jclass+1) 1607 ENDIF 1608 !salinity classes transports 1609 IF( ( sec%zsal(jclass) .NE. 99._wp ) .AND. & 1610 ( sec%zsal(jclass+1) .NE. 99._wp ) )THEN 1611 classe = 'S ' 1612 zbnd1 = sec%zsal(jclass) 1613 zbnd2 = sec%zsal(jclass+1) 1614 ENDIF 1615 !temperature classes transports 1616 IF( ( sec%ztem(jclass) .NE. 99._wp ) .AND. & 1617 ( sec%ztem(jclass+1) .NE. 99._wp ) ) THEN 1618 classe = 'T ' 1619 zbnd1 = sec%ztem(jclass) 1620 zbnd2 = sec%ztem(jclass+1) 1621 ENDIF 1622 1623 !write volume transport per class 1624 IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 1625 WRITE(numdct_NOOS,'(9e12.4E2)') -(sec%transport( 1,jclass)+sec%transport( 2,jclass)),-sec%transport( 2,jclass),-sec%transport( 1,jclass), & 1626 -(sec%transport( 7,jclass)+sec%transport( 8,jclass)),-sec%transport( 8,jclass),-sec%transport( 7,jclass), & 1627 -(sec%transport( 9,jclass)+sec%transport(10,jclass)),-sec%transport(10,jclass),-sec%transport( 9,jclass) 1628 ELSE 1629 WRITE(numdct_NOOS,'(9e12.4E2)') sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), & 1630 sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), & 1631 sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass) 1632 ENDIF 1633 1634 ENDDO 1635 1636 CALL wrk_dealloc(nb_type , zsumclasses ) 1637 1638 END SUBROUTINE dia_dct_wri_NOOS 1639 1640 SUBROUTINE dia_dct_wri_NOOS_h(hr,ksec,sec) 1641 !!------------------------------------------------------------- 1642 !! As routine dia_dct_wri_NOOS but for hourly output files 1643 !! 1644 !! Write transport output in numdct using NOOS formatting 1645 !! 1646 !! Purpose: Write transports in ascii files 1647 !! 1648 !! Method: 1649 !! 1. Write volume transports in "volume_transport" 1650 !! Unit: Sv : area * Velocity / 1.e6 1651 !! 1652 !!------------------------------------------------------------- 1653 !!arguments 1654 INTEGER, INTENT(IN) :: hr ! hour 1655 TYPE(SECTION), INTENT(INOUT) :: sec ! section to write 1656 INTEGER ,INTENT(IN) :: ksec ! section number 1657 1658 !!local declarations 1659 INTEGER :: jclass,jhr ! Dummy loop 1660 CHARACTER(len=2) :: classe ! Classname 1661 REAL(wp) :: zbnd1,zbnd2 ! Class bounds 1662 REAL(wp) :: zslope ! section's slope coeff 1663 ! 1664 REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace 1665 !!------------------------------------------------------------- 1666 1667 CALL wrk_alloc(nb_type , zsumclasses ) 1668 1669 zsumclasses(:)=0._wp 1670 zslope = sec%slopeSection 1671 1672 DO jclass=1,MAX(1,sec%nb_class-1) 1673 zsumclasses(1:nb_type)=zsumclasses(1:nb_type)+sec%transport_h(1:nb_type,jclass) 1674 ENDDO 1675 1676 !write volume transport per class 1677 IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 1678 z_hr_output(ksec,hr,1)=-(zsumclasses(1)+zsumclasses(2)) 1679 ELSE 1680 z_hr_output(ksec,hr,1)= (zsumclasses(1)+zsumclasses(2)) 1681 ENDIF 1682 1683 DO jclass=1,MAX(1,sec%nb_class-1) 1684 IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN 1685 z_hr_output(ksec,hr,jclass+1)=-(sec%transport_h(1,jclass)+sec%transport_h(2,jclass)) 1686 ELSE 1687 z_hr_output(ksec,hr,jclass+1)= (sec%transport_h(1,jclass)+sec%transport_h(2,jclass)) 1688 ENDIF 1689 ENDDO 1690 1691 IF ( hr .eq. 48._wp ) THEN 1692 WRITE(numdct_NOOS_h,'(I4,a1,I2,a1,I2,a12,i3,a17,i3)') nyear,'.',nmonth,'.',nday,' Transect:',ksec-1,' No. of layers:',sec%nb_class-1 1693 DO jhr=25,48 1694 WRITE(numdct_NOOS_h,'(11F12.1)') z_hr_output(ksec,jhr,1), (z_hr_output(ksec,jhr,jclass+1), jclass=1,MAX(1,10) ) 1695 ENDDO 1696 ENDIF 1697 1698 CALL wrk_dealloc(nb_type , zsumclasses ) 1699 1700 END SUBROUTINE dia_dct_wri_NOOS_h 1701 896 1702 SUBROUTINE dia_dct_wri(kt,ksec,sec) 897 1703 !!------------------------------------------------------------- … … 917 1723 918 1724 !!local declarations 919 INTEGER :: jcl ,ji! Dummy loop1725 INTEGER :: jclass ! Dummy loop 920 1726 CHARACTER(len=2) :: classe ! Classname 921 1727 REAL(wp) :: zbnd1,zbnd2 ! Class bounds 922 1728 REAL(wp) :: zslope ! section's slope coeff 923 1729 ! 924 REAL(wp), POINTER, DIMENSION(:):: zsumclass ! 1D workspace1730 REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace 925 1731 !!------------------------------------------------------------- 926 CALL wrk_alloc(nb_type _class , zsumclass )927 928 zsumclass (:)=0._wp1732 CALL wrk_alloc(nb_type , zsumclasses ) 1733 1734 zsumclasses(:)=0._wp 929 1735 zslope = sec%slopeSection 930 931 932 DO jcl=1,MAX(1,sec%nb_class-1) 933 934 ! Mean computation 935 sec%transport(:,jcl)=sec%transport(:,jcl)/(nn_dctwri/nn_dct) 1736 1737 DO jclass=1,MAX(1,sec%nb_class-1) 1738 936 1739 classe = 'N ' 937 1740 zbnd1 = 0._wp 938 1741 zbnd2 = 0._wp 939 zsumclass (1:nb_type_class)=zsumclass(1:nb_type_class)+sec%transport(1:nb_type_class,jcl)1742 zsumclasses(1:nb_type)=zsumclasses(1:nb_type)+sec%transport(1:nb_type,jclass) 940 1743 941 1744 942 1745 !insitu density classes transports 943 IF( ( sec%zsigi(jcl ) .NE. 99._wp ) .AND. &944 ( sec%zsigi(jcl +1) .NE. 99._wp ) )THEN1746 IF( ( sec%zsigi(jclass) .NE. 99._wp ) .AND. & 1747 ( sec%zsigi(jclass+1) .NE. 99._wp ) )THEN 945 1748 classe = 'DI ' 946 zbnd1 = sec%zsigi(jcl )947 zbnd2 = sec%zsigi(jcl +1)1749 zbnd1 = sec%zsigi(jclass) 1750 zbnd2 = sec%zsigi(jclass+1) 948 1751 ENDIF 949 1752 !potential density classes transports 950 IF( ( sec%zsigp(jcl ) .NE. 99._wp ) .AND. &951 ( sec%zsigp(jcl +1) .NE. 99._wp ) )THEN1753 IF( ( sec%zsigp(jclass) .NE. 99._wp ) .AND. & 1754 ( sec%zsigp(jclass+1) .NE. 99._wp ) )THEN 952 1755 classe = 'DP ' 953 zbnd1 = sec%zsigp(jcl )954 zbnd2 = sec%zsigp(jcl +1)1756 zbnd1 = sec%zsigp(jclass) 1757 zbnd2 = sec%zsigp(jclass+1) 955 1758 ENDIF 956 1759 !depth classes transports 957 IF( ( sec%zlay(jcl ) .NE. 99._wp ) .AND. &958 ( sec%zlay(jcl +1) .NE. 99._wp ) )THEN1760 IF( ( sec%zlay(jclass) .NE. 99._wp ) .AND. & 1761 ( sec%zlay(jclass+1) .NE. 99._wp ) )THEN 959 1762 classe = 'Z ' 960 zbnd1 = sec%zlay(jcl )961 zbnd2 = sec%zlay(jcl +1)1763 zbnd1 = sec%zlay(jclass) 1764 zbnd2 = sec%zlay(jclass+1) 962 1765 ENDIF 963 1766 !salinity classes transports 964 IF( ( sec%zsal(jcl ) .NE. 99._wp ) .AND. &965 ( sec%zsal(jcl +1) .NE. 99._wp ) )THEN1767 IF( ( sec%zsal(jclass) .NE. 99._wp ) .AND. & 1768 ( sec%zsal(jclass+1) .NE. 99._wp ) )THEN 966 1769 classe = 'S ' 967 zbnd1 = sec%zsal(jcl )968 zbnd2 = sec%zsal(jcl +1)1770 zbnd1 = sec%zsal(jclass) 1771 zbnd2 = sec%zsal(jclass+1) 969 1772 ENDIF 970 1773 !temperature classes transports 971 IF( ( sec%ztem(jcl ) .NE. 99._wp ) .AND. &972 ( sec%ztem(jcl +1) .NE. 99._wp ) ) THEN1774 IF( ( sec%ztem(jclass) .NE. 99._wp ) .AND. & 1775 ( sec%ztem(jclass+1) .NE. 99._wp ) ) THEN 973 1776 classe = 'T ' 974 zbnd1 = sec%ztem(jcl )975 zbnd2 = sec%ztem(jcl +1)1777 zbnd1 = sec%ztem(jclass) 1778 zbnd2 = sec%ztem(jclass+1) 976 1779 ENDIF 977 1780 978 1781 !write volume transport per class 979 1782 WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, & 980 jcl ,classe,zbnd1,zbnd2,&981 sec%transport(1,jcl ),sec%transport(2,jcl), &982 sec%transport(1,jcl )+sec%transport(2,jcl)1783 jclass,classe,zbnd1,zbnd2,& 1784 sec%transport(1,jclass),sec%transport(2,jclass), & 1785 sec%transport(1,jclass)+sec%transport(2,jclass) 983 1786 984 1787 IF( sec%llstrpond )THEN 985 1788 986 1789 !write heat transport per class: 987 WRITE(numdct_ heat,119) ndastp,kt,ksec,sec%name,zslope, &988 jcl ,classe,zbnd1,zbnd2,&989 sec%transport(7,jcl )*1000._wp*rcp/1.e15,sec%transport(8,jcl)*1000._wp*rcp/1.e15, &990 ( sec%transport(7,jcl )+sec%transport(8,jcl) )*1000._wp*rcp/1.e151790 WRITE(numdct_temp,119) ndastp,kt,ksec,sec%name,zslope, & 1791 jclass,classe,zbnd1,zbnd2,& 1792 sec%transport(7,jclass)*1000._wp*rcp/1.e15,sec%transport(8,jclass)*1000._wp*rcp/1.e15, & 1793 ( sec%transport(7,jclass)+sec%transport(8,jclass) )*1000._wp*rcp/1.e15 991 1794 !write salt transport per class 992 WRITE(numdct_sal t,119) ndastp,kt,ksec,sec%name,zslope, &993 jcl ,classe,zbnd1,zbnd2,&994 sec%transport(9,jcl )*1000._wp/1.e9,sec%transport(10,jcl)*1000._wp/1.e9,&995 (sec%transport(9,jcl )+sec%transport(10,jcl))*1000._wp/1.e91795 WRITE(numdct_sal ,119) ndastp,kt,ksec,sec%name,zslope, & 1796 jclass,classe,zbnd1,zbnd2,& 1797 sec%transport(9,jclass)*1000._wp/1.e9,sec%transport(10,jclass)*1000._wp/1.e9,& 1798 (sec%transport(9,jclass)+sec%transport(10,jclass))*1000._wp/1.e9 996 1799 ENDIF 997 1800 … … 1000 1803 zbnd1 = 0._wp 1001 1804 zbnd2 = 0._wp 1002 jcl =01805 jclass=0 1003 1806 1004 1807 !write total volume transport 1005 1808 WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, & 1006 jcl ,"total",zbnd1,zbnd2,&1007 zsumclass (1),zsumclass(2),zsumclass(1)+zsumclass(2)1809 jclass,"total",zbnd1,zbnd2,& 1810 zsumclasses(1),zsumclasses(2),zsumclasses(1)+zsumclasses(2) 1008 1811 1009 1812 IF( sec%llstrpond )THEN 1010 1813 1011 1814 !write total heat transport 1012 WRITE(numdct_ heat,119) ndastp,kt,ksec,sec%name,zslope, &1013 jcl ,"total",zbnd1,zbnd2,&1014 zsumclass (7)* 1000._wp*rcp/1.e15,zsumclass(8)* 1000._wp*rcp/1.e15,&1015 (zsumclass (7)+zsumclass(8) )* 1000._wp*rcp/1.e151815 WRITE(numdct_temp,119) ndastp,kt,ksec,sec%name,zslope, & 1816 jclass,"total",zbnd1,zbnd2,& 1817 zsumclasses(7)* 1000._wp*rcp/1.e15,zsumclasses(8)* 1000._wp*rcp/1.e15,& 1818 (zsumclasses(7)+zsumclasses(8) )* 1000._wp*rcp/1.e15 1016 1819 !write total salt transport 1017 WRITE(numdct_sal t,119) ndastp,kt,ksec,sec%name,zslope, &1018 jcl ,"total",zbnd1,zbnd2,&1019 zsumclass (9)*1000._wp/1.e9,zsumclass(10)*1000._wp/1.e9,&1020 (zsumclass (9)+zsumclass(10))*1000._wp/1.e91820 WRITE(numdct_sal ,119) ndastp,kt,ksec,sec%name,zslope, & 1821 jclass,"total",zbnd1,zbnd2,& 1822 zsumclasses(9)*1000._wp/1.e9,zsumclasses(10)*1000._wp/1.e9,& 1823 (zsumclasses(9)+zsumclasses(10))*1000._wp/1.e9 1021 1824 ENDIF 1022 1825 … … 1025 1828 !write total ice volume transport 1026 1829 WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,& 1027 jcl ,"ice_vol",zbnd1,zbnd2,&1028 sec%transport( 9,1),sec%transport(10,1),&1029 sec%transport( 9,1)+sec%transport(10,1)1830 jclass,"ice_vol",zbnd1,zbnd2,& 1831 sec%transport(11,1),sec%transport(12,1),& 1832 sec%transport(11,1)+sec%transport(12,1) 1030 1833 !write total ice surface transport 1031 1834 WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,& 1032 jcl ,"ice_surf",zbnd1,zbnd2,&1033 sec%transport(1 1,1),sec%transport(12,1), &1034 sec%transport(1 1,1)+sec%transport(12,1)1835 jclass,"ice_surf",zbnd1,zbnd2,& 1836 sec%transport(13,1),sec%transport(14,1), & 1837 sec%transport(13,1)+sec%transport(14,1) 1035 1838 ENDIF 1036 1839 … … 1038 1841 119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6) 1039 1842 1040 CALL wrk_dealloc(nb_type _class , zsumclass )1843 CALL wrk_dealloc(nb_type , zsumclasses ) 1041 1844 END SUBROUTINE dia_dct_wri 1042 1845 -
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DIA/diaptr.F90
r7363 r7367 332 332 !!---------------------------------------------------------------------- 333 333 USE oce, vt => ua ! use ua as workspace 334 USE oce, vs => ua ! use ua as workspace334 USE oce, vs => va ! use va as workspace 335 335 IMPLICIT none 336 336 !! … … 378 378 zv = ( vn(ji,jj,jk) + vn(ji,jj-1,jk) ) * 0.5_wp 379 379 #endif 380 vt( :,jj,jk) = zv * tsn(:,jj,jk,jp_tem)381 vs( :,jj,jk) = zv * tsn(:,jj,jk,jp_sal)380 vt(ji,jj,jk) = zv * tsn(ji,jj,jk,jp_tem) 381 vs(ji,jj,jk) = zv * tsn(ji,jj,jk,jp_sal) 382 382 END DO 383 383 END DO -
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
r7363 r7367 44 44 USE iom 45 45 USE ioipsl 46 USE diafoam, ONLY: dia_wri_foam 47 !CEOD USE insitu_tem, ONLY: insitu_t, theta2t 48 USE bartrop_uv, ONLY: un_dm, vn_dm, bartrop_vel 46 49 #if defined key_lim2 47 50 USE limwri_2 51 USE ice_2 ! LIM_2 ice model variables 52 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 53 #endif 54 #if defined key_lim3 55 USE ice_3 ! LIM_3 ice model variables 56 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 57 #endif 58 USE daymod ! calendar 59 USE insitu_tem, ONLY: insitu_t, theta2t 60 #if defined key_top 61 USE par_trc ! biogeochemical variables 62 USE trc 63 #endif 64 #if defined key_spm 65 USE spm_con, ONLY: Eps0XS 66 #endif 67 #if defined key_zdftke 68 USE zdftke, ONLY: en 69 #endif 70 USE zdf_oce, ONLY: avt, avm 71 #if defined key_zdfgls 72 USE zdfgls, ONLY: mxln, en 48 73 #endif 49 74 USE lib_mpp ! MPP library … … 54 79 PRIVATE 55 80 81 PUBLIC dia_wri_tmb_init ! Called by nemogcm module 56 82 PUBLIC dia_wri ! routines called by step.F90 57 83 PUBLIC dia_wri_state 58 84 PUBLIC dia_wri_alloc ! Called by nemogcm module 85 PUBLIC dia_wri_tide_init ! Called by nemogcm module 59 86 60 87 INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file … … 65 92 INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV 66 93 INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V 94 95 !! * variables for calculating 25-hourly means 96 REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: tn_25h , sn_25h , insitu_t_25h 97 REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:) :: sshn_25h, hmld_kara_25h 98 REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: un_25h , vn_25h , wn_25h 99 REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avt_25h , avm_25h 100 #if defined key_zdfgls || key_zdftke 101 REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: en_25h 102 #endif 103 #if defined key_zdfgls 104 REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: mxln_25h 105 #endif 106 INTEGER, SAVE :: cnt_25h ! Counter for 25 hour means 107 108 67 109 68 110 !! * Substitutions … … 125 167 REAL(wp), POINTER, DIMENSION(:,:) :: z2d ! 2D workspace 126 168 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d ! 3D workspace 169 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmb ! temporary workspace for tmb 127 170 !!---------------------------------------------------------------------- 128 171 ! … … 131 174 CALL wrk_alloc( jpi , jpj , z2d ) 132 175 CALL wrk_alloc( jpi , jpj, jpk , z3d ) 176 CALL wrk_alloc( jpi , jpj, 3 , zwtmb ) 133 177 ! 134 178 ! Output the initial state and forcings … … 138 182 ENDIF 139 183 184 IF (ln_diatide) THEN 185 CALL dia_wri_tide(kt) 186 ENDIF 187 140 188 CALL iom_put( "toce" , tsn(:,:,:,jp_tem) ) ! temperature 189 CALL theta2t ! in-situ temperature conversion 190 !CEOD CALL iom_put( "tinsitu", insitu_t(:,:,:) ) ! in-situ temperature 141 191 CALL iom_put( "soce" , tsn(:,:,:,jp_sal) ) ! salinity 142 192 CALL iom_put( "sst" , tsn(:,:,1,jp_tem) ) ! sea surface temperature … … 146 196 CALL iom_put( "uoce" , un ) ! i-current 147 197 CALL iom_put( "voce" , vn ) ! j-current 148 198 CALL iom_put( "ssu" , un(:,:,1) ) ! sea surface U velocity 199 CALL iom_put( "ssv" , vn(:,:,1) ) ! sea surface V velocity 200 IF( cp_cfg == "natl" .OR. cp_cfg == "ind12" ) CALL bartrop_vel ! barotropic velocity conversion 201 !These don't exist independently in this branch so we remove them to get a CO5 202 !that works on the Cray 203 !CEOD CALL iom_put( "uocebt" , un_dm ) ! barotropic i-current 204 !CEOD CALL iom_put( "vocebt" , vn_dm ) ! barotropic j-current 149 205 CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef. 150 206 CALL iom_put( "avm" , avmu ) ! T vert. eddy visc. coef. 151 207 IF( lk_zdfddm ) THEN 152 208 CALL iom_put( "avs" , fsavs(:,:,:) ) ! S vert. eddy diff. coef. 209 ENDIF 210 ! 211 ! If we want tmb values 212 213 IF (ln_diatmb) THEN 214 CALL dia_wri_calctmb( tsn(:,:,:,jp_tem),zwtmb ) 215 !ssh already output but here we output it masked 216 CALL iom_put( "sshnmasked" , sshn(:,:)*tmask(:,:,1)+missing_val*(1-tmask(:,:,1 ) ) ) ! tmb Temperature 217 CALL iom_put( "top_temp" , zwtmb(:,:,1) ) ! tmb Temperature 218 CALL iom_put( "mid_temp" , zwtmb(:,:,2) ) ! tmb Temperature 219 CALL iom_put( "bot_temp" , zwtmb(:,:,3) ) ! tmb Temperature 220 ! CALL iom_put( "sotrefml" , hmld_tref(:,:) ) ! "T criterion Mixed Layer Depth 221 222 CALL dia_wri_calctmb( tsn(:,:,:,jp_sal),zwtmb ) 223 CALL iom_put( "top_sal" , zwtmb(:,:,1) ) ! tmb Salinity 224 CALL iom_put( "mid_sal" , zwtmb(:,:,2) ) ! tmb Salinity 225 CALL iom_put( "bot_sal" , zwtmb(:,:,3) ) ! tmb Salinity 226 227 CALL dia_wri_calctmb( un(:,:,:),zwtmb ) 228 CALL iom_put( "top_u" , zwtmb(:,:,1) ) ! tmb U Velocity 229 CALL iom_put( "mid_u" , zwtmb(:,:,2) ) ! tmb U Velocity 230 CALL iom_put( "bot_u" , zwtmb(:,:,3) ) ! tmb U Velocity 231 !Called in dynspg_ts.F90 CALL iom_put( "baro_u" , un_b ) ! Barotropic U Velocity 232 233 CALL dia_wri_calctmb( vn(:,:,:),zwtmb ) 234 CALL iom_put( "top_v" , zwtmb(:,:,1) ) ! tmb V Velocity 235 CALL iom_put( "mid_v" , zwtmb(:,:,2) ) ! tmb V Velocity 236 CALL iom_put( "bot_v" , zwtmb(:,:,3) ) ! tmb V Velocity 237 !Called in dynspg_ts.F90 CALL iom_put( "baro_v" , vn_b ) ! Barotropic V Velocity 153 238 ENDIF 154 239 … … 171 256 z3d(:,:,jpk) = 0.e0 172 257 DO jk = 1, jpkm1 173 z3d(:,:,jk) = rau0 * un(:,:,jk) * e 1u(:,:) * fse3u(:,:,jk)258 z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) 174 259 END DO 175 260 CALL iom_put( "u_masstr", z3d ) ! mass transport in i-direction … … 186 271 CALL iom_put( "u_heattr", z2d ) ! heat transport in i-direction 187 272 DO jk = 1, jpkm1 188 z3d(:,:,jk) = rau0 * vn(:,:,jk) * e 2v(:,:) * fse3v(:,:,jk)273 z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) 189 274 END DO 190 275 CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction … … 251 336 ENDIF 252 337 ! 338 ! -1. Alternative routine 339 !------------------------ 340 IF (ln_diafoam) THEN 341 CALL dia_wri_foam( kt ) 342 RETURN 343 END IF 344 ! 253 345 ! 0. Initialisation 254 346 ! ----------------- … … 673 765 #endif 674 766 767 SUBROUTINE dia_wri_calctmb( infield,outtmb ) 768 !!--------------------------------------------------------------------- 769 !! *** ROUTINE dia_tmb *** 770 !! 771 !! ** Purpose : Write diagnostics for Top,Mid, and Bottom of water Column 772 !! 773 !! ** Method : 774 !! use mbathy to find surface, mid and bottom of model levels 775 !! 776 !! History : 777 !! 3.4 ! 04-13 (E. O'Dea) Routine taken from old dia_wri_foam 778 !!---------------------------------------------------------------------- 779 !! * Modules used 780 781 ! Routine to map 3d field to top, middle, bottom 782 IMPLICIT NONE 783 784 ! Routine arguments 785 REAL(wp), DIMENSION(jpi, jpj, jpk), INTENT(IN ) :: infield ! Input 3d field and mask 786 REAL(wp), DIMENSION(jpi, jpj, 3 ), INTENT( OUT) :: outtmb ! Output top, middle, bottom 787 788 ! Local variables 789 INTEGER :: ji,jj,jk ! Dummy loop indices 790 791 ! Calculate top 792 outtmb(:,:,1) = infield(:,:,1)*tmask(:,:,1) + missing_val*(1-tmask(:,:,1)) 793 794 ! Calculate middle 795 DO ji = 1,jpi 796 DO jj = 1,jpj 797 jk = max(1,mbathy(ji,jj)/2) 798 outtmb(ji,jj,2) = infield(ji,jj,jk)*tmask(ji,jj,jk) + missing_val*(1-tmask(ji,jj,jk)) 799 END DO 800 END DO 801 802 ! Calculate bottom 803 DO ji = 1,jpi 804 DO jj = 1,jpj 805 jk = max(1,mbathy(ji,jj) - 1) 806 outtmb(ji,jj,3) = infield(ji,jj,jk)*tmask(ji,jj,jk) + missing_val*(1-tmask(ji,jj,jk)) 807 END DO 808 END DO 809 810 END SUBROUTINE dia_wri_calctmb 811 812 SUBROUTINE dia_wri_tmb_init 813 !!--------------------------------------------------------------------------- 814 !! *** ROUTINE dia_wri_tmb_init *** 815 !! 816 !! ** Purpose: Initialization of tmb namelist 817 !! 818 !! ** Method : Read namelist 819 !! History 820 !! 3.4 ! 04-13 (E. O'Dea) Routine to initialize dia_wri_tmb 821 !!--------------------------------------------------------------------------- 822 !! 823 INTEGER :: ierror ! local integer 824 !! 825 NAMELIST/nam_diatmb/ ln_diatmb 826 !! 827 !!---------------------------------------------------------------------- 828 ! 829 REWIND ( numnam ) ! Read Namelist nam_diatmb 830 READ ( numnam, nam_diatmb ) 831 ! 832 IF(lwp) THEN ! Control print 833 WRITE(numout,*) 834 WRITE(numout,*) 'dia_wri_tmb_init : Output Top, Middle, Bottom Diagnostics' 835 WRITE(numout,*) '~~~~~~~~~~~~' 836 WRITE(numout,*) ' Namelist nam_diatmb : set tmb outputs ' 837 WRITE(numout,*) ' Switch for TMB diagnostics (T) or not (F) ln_diatmb = ', ln_diatmb 838 ENDIF 839 840 END SUBROUTINE dia_wri_tmb_init 841 842 675 843 SUBROUTINE dia_wri_state( cdfile_name, kt ) 676 844 !!--------------------------------------------------------------------- … … 798 966 END SUBROUTINE dia_wri_state 799 967 !!====================================================================== 968 !!====================================================================== 969 970 SUBROUTINE dia_wri_tide( kt ) 971 !!--------------------------------------------------------------------- 972 !! *** ROUTINE dia_tide *** 973 !! 974 !! ** Purpose : Write diagnostics with M2/S2 tide removed 975 !! 976 !! ** Method : 977 !! 25hr mean outputs for shelf seas 978 !! 979 !! History : 980 !! ?.0 ! 07-04 (A. Hines) New routine, developed from dia_wri_foam 981 !! 3.4 ! 02-13 (J. Siddorn) Routine taken from old dia_wri_foam 982 !!---------------------------------------------------------------------- 983 !! * Modules used 984 985 IMPLICIT NONE 986 987 !! * Arguments 988 INTEGER, INTENT( in ) :: kt ! ocean time-step index 989 990 991 !! * Local declarations 992 INTEGER :: ji, jj, jk 993 994 LOGICAL :: ll_print = .FALSE. ! =T print and flush numout 995 REAL(wp) :: zsto, zout, zmax, zjulian, zdt, zmdi ! temporary integers 996 INTEGER :: i_steps ! no of timesteps per hour 997 REAL(wp), DIMENSION(jpi,jpj ) :: zw2d, un_dm, vn_dm ! temporary workspace 998 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! temporary workspace 999 REAL(wp), DIMENSION(jpi,jpj,3) :: zwtmb ! temporary workspace 1000 INTEGER :: nyear0, nmonth0,nday0 ! start year,month,day 1001 !#if defined key_top 1002 ! CHARACTER (len=20) :: cltra, cltrau 1003 ! CHARACTER (len=80) :: cltral 1004 ! INTEGER :: jn, jl 1005 !#endif 1006 !#if defined key_spm 1007 ! ! variables needed to calculate visibility field from sediment fields 1008 ! REAL(wp), DIMENSION(jpi,jpj,jpk) :: vis3d ! derived 3D visibility field 1009 ! REAL(wp) :: epsessX = 0.07d-03 ! attenuation coefficient applied to the sediment (as used in ERSEM) 1010 ! REAL(wp) :: tiny = 1.0d-15 ! to prevent division by zero in visibility calculation 1011 !#endif 1012 1013 !!---------------------------------------------------------------------- 1014 1015 ! 0. Initialisation 1016 ! ----------------- 1017 ! Define frequency of summing to create 25 h mean 1018 zdt = rdt 1019 IF( nacc == 1 ) zdt = rdtmin 1020 1021 IF( MOD( 3600,INT(zdt) ) == 0 ) THEN 1022 i_steps = 3600/INT(zdt) 1023 ELSE 1024 CALL ctl_stop('STOP', 'dia_wri_tide: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible') 1025 ENDIF 1026 1027 #if defined key_lim3 || defined key_lim2 1028 CALL ctl_stop('STOP', 'dia_wri_tide not setup yet to do tidemean ice') 1029 #endif 1030 #if defined key_spm || defined key_MOersem 1031 CALL ctl_stop('STOP', 'dia_wri_tide not setup yet to do tidemean ERSEM') 1032 #endif 1033 1034 ! local variable for debugging 1035 ll_print = ll_print .AND. lwp 1036 1037 ! Sum of 25 hourly instantaneous values to give a 25h mean from 24hours 1038 ! every day 1039 IF( MOD( kt, i_steps ) == 0 .and. kt .ne. nn_it000 ) THEN 1040 1041 IF (lwp) THEN 1042 WRITE(numout,*) 'dia_wri_tide : Summing instantaneous hourly diagnostics at timestep ',kt 1043 WRITE(numout,*) '~~~~~~~~~~~~ ' 1044 ENDIF 1045 1046 tn_25h(:,:,:) = tn_25h(:,:,:) + tsn(:,:,:,jp_tem) 1047 sn_25h(:,:,:) = sn_25h(:,:,:) + tsn(:,:,:,jp_sal) 1048 CALL theta2t 1049 insitu_t_25h(:,:,:) = insitu_t_25h(:,:,:) + insitu_t(:,:,:) 1050 sshn_25h(:,:) = sshn_25h(:,:) + sshn (:,:) 1051 ! hmld_kara_25h(:,:) = hmld_kara_25h(:,:) + hmld_kara(:,:) 1052 #if defined key_lim3 || defined key_lim2 1053 hsnif_25h(:,:) = hsnif_25h(:,:) + hsnif(:,:) 1054 hicif_25h(:,:) = hicif_25h(:,:) + hicif(:,:) 1055 frld_25h(:,:) = frld_25h(:,:) + frld(:,:) 1056 #endif 1057 #if defined key_spm || defined key_MOersem 1058 trn_25h(:,:,:,:) = trn_25h(:,:,:,:) + trn (:,:,:,:) 1059 trc3d_25h(:,:,:,:) = trc3d_25h(:,:,:,:) + trc3d(:,:,:,:) 1060 trc2d_25h(:,:,:) = trc2d_25h(:,:,:) + trc2d(:,:,:) 1061 #endif 1062 un_25h(:,:,:) = un_25h(:,:,:) + un(:,:,:) 1063 vn_25h(:,:,:) = vn_25h(:,:,:) + vn(:,:,:) 1064 wn_25h(:,:,:) = wn_25h(:,:,:) + wn(:,:,:) 1065 avt_25h(:,:,:) = avt_25h(:,:,:) + avt(:,:,:) 1066 avm_25h(:,:,:) = avm_25h(:,:,:) + avm(:,:,:) 1067 # if defined key_zdfgls || defined key_zdftke 1068 en_25h(:,:,:) = en_25h(:,:,:) + en(:,:,:) 1069 #endif 1070 # if defined key_zdfgls 1071 mxln_25h(:,:,:) = mxln_25h(:,:,:) + mxln(:,:,:) 1072 #endif 1073 cnt_25h = cnt_25h + 1 1074 1075 IF (lwp) THEN 1076 WRITE(numout,*) 'dia_tide : Summed the following number of hourly values so far',cnt_25h 1077 ENDIF 1078 1079 ENDIF ! MOD( kt, i_steps ) == 0 1080 1081 ! Write data for 25 hour mean output streams 1082 IF( cnt_25h .EQ. 25 .AND. MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN 1083 1084 IF(lwp) THEN 1085 WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt 1086 WRITE(numout,*) '~~~~~~~~~~~~ ' 1087 ENDIF 1088 1089 tn_25h(:,:,:) = tn_25h(:,:,:) / 25.0_wp 1090 sn_25h(:,:,:) = sn_25h(:,:,:) / 25.0_wp 1091 insitu_t_25h(:,:,:) = insitu_t_25h(:,:,:) / 25.0_wp 1092 sshn_25h(:,:) = sshn_25h(:,:) / 25.0_wp 1093 ! hmld_kara_25h(:,:) = hmld_kara_25h(:,:) / 25.0_wp 1094 #if defined key_lim3 || defined key_lim2 1095 hsnif_25h(:,:) = hsnif_25h(:,:) / 25.0_wp 1096 hicif_25h(:,:) = hicif_25h(:,:) / 25.0_wp 1097 frld_25h(:,:) = frld_25h(:,:) / 25.0_wp 1098 #endif 1099 #if defined key_spm || defined key_MOersem 1100 trn_25h(:,:,:,:) = trn_25h(:,:,:,:) / 25.0_wp 1101 trc3d_25h(:,:,:,:) = trc3d_25h(:,:,:,:) / 25.0_wp 1102 trc2d_25h(:,:,:) = trc2d_25h(:,:,:) / 25.0_wp 1103 #endif 1104 un_25h(:,:,:) = un_25h(:,:,:) / 25.0_wp 1105 vn_25h(:,:,:) = vn_25h(:,:,:) / 25.0_wp 1106 wn_25h(:,:,:) = wn_25h(:,:,:) / 25.0_wp 1107 avt_25h(:,:,:) = avt_25h(:,:,:) / 25.0_wp 1108 avm_25h(:,:,:) = avm_25h(:,:,:) / 25.0_wp 1109 # if defined key_zdfgls || defined key_zdftke 1110 en_25h(:,:,:) = en_25h(:,:,:) / 25.0_wp 1111 #endif 1112 # if defined key_zdfgls 1113 mxln_25h(:,:,:) = mxln_25h(:,:,:) / 25.0_wp 1114 #endif 1115 1116 IF (lwp) WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output' 1117 1118 zmdi=missing_val ! for masking 1119 ! write tracers (instantaneous) 1120 zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 1121 CALL iom_put("temper25h", zw3d) ! potential temperature 1122 CALL theta2t ! calculate insitu temp 1123 zw3d(:,:,:) = insitu_t_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 1124 CALL iom_put("tempis25h", zw3d) ! in-situ temperature 1125 zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 1126 CALL iom_put( "salin25h", zw3d ) ! salinity 1127 zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 1128 CALL iom_put( "ssh25h", zw2d ) ! sea surface 1129 ! zw2d(:,:) = hmld_kara_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 1130 ! CALL iom_put( "kara25h", zw2d ) ! mixed layer 1131 1132 #if defined key_lim3 || defined key_lim2 1133 ! Write ice model variables (instantaneous) 1134 zw2d(:,:) = hsnif_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 1135 CALL iom_put("isnowthi", zw2d ) ! ice thickness 1136 zw2d(:,:) = hicif_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 1137 CALL iom_put("iicethic", zw2d ) ! ice thickness 1138 zw2d(:,:) = (1.0-frld_25h(:,:))*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 1139 CALL iom_put("iiceconc", zw2d ) ! ice concetration 1140 #endif 1141 #if defined key_spm || key_MOersem 1142 ! output biogeochemical variables: 1143 ! output main tracers 1144 DO jn = 1, jptra 1145 cltra = ctrcnm(jn) ! short title for tracer 1146 zw3d(:,:,:) = trn_25h(:,:,:,jn)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 1147 IF( lutsav(jn) ) CALL iom_put( cltra, zw3d ) ! temperature 1148 END DO 1149 ! more 3D horizontal arrays from diagnostics 1150 DO jl = 1, jpdia3d 1151 cltra = ctrc3d(jl) ! short title for 3D diagnostic 1152 zw3d(:,:,:) = trc3d_25h(:,:,:,jl)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 1153 CALL iom_put( cltra, zw3d ) 1154 END DO 1155 ! more 2D horizontal arrays from diagnostics 1156 DO jl = 1, jpdia2d 1157 cltra = ctrc2d(jl) ! short title for 2D diagnostic 1158 zw2d(:,:) = trc2d_25h(:,:,jl)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 1159 CALL iom_put(cltra, zw2d ) 1160 END DO 1161 #endif 1162 1163 ! Write velocities (instantaneous) 1164 zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:)) 1165 CALL iom_put("vozocrtx25h", zw3d) ! i-current 1166 zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:)) 1167 CALL iom_put("vomecrty25h", zw3d ) ! j-current 1168 1169 zw3d(:,:,:) = wn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 1170 CALL iom_put("vomecrtz25h", zw3d ) ! k-current 1171 zw3d(:,:,:) = avt_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 1172 CALL iom_put("avt25h", zw3d ) ! diffusivity 1173 zw3d(:,:,:) = avm_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 1174 CALL iom_put("avm25h", zw3d) ! viscosity 1175 #if defined key_zdftke || defined key_zdfgls 1176 zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 1177 CALL iom_put("tke25h", zw3d) ! tke 1178 #endif 1179 #if defined key_zdfgls 1180 zw3d(:,:,:) = mxln_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 1181 CALL iom_put( "mxln25h",zw3d) 1182 #endif 1183 1184 ! After the write reset the values to cnt=1 and sum values equal current value 1185 tn_25h(:,:,:) = tsn(:,:,:,jp_tem) 1186 sn_25h(:,:,:) = tsn(:,:,:,jp_sal) 1187 CALL theta2t 1188 insitu_t_25h(:,:,:) = insitu_t(:,:,:) 1189 sshn_25h(:,:) = sshn (:,:) 1190 ! hmld_kara_25h(:,:) = hmld_kara(:,:) 1191 #if defined key_lim3 || defined key_lim2 1192 hsnif_25h(:,:) = hsnif(:,:) 1193 hicif_25h(:,:) = hicif(:,:) 1194 frld_25h(:,:) = frld(:,:) 1195 #endif 1196 #if defined key_spm || defined key_MOersem 1197 trn_25h(:,:,:,:) = trn (:,:,:,:) 1198 trc3d_25h(:,:,:,:) = trc3d(:,:,:,:) 1199 trc2d_25h(:,:,:) = trc2d(:,:,:) 1200 #endif 1201 un_25h(:,:,:) = un(:,:,:) 1202 vn_25h(:,:,:) = vn(:,:,:) 1203 wn_25h(:,:,:) = wn(:,:,:) 1204 avt_25h(:,:,:) = avt(:,:,:) 1205 avm_25h(:,:,:) = avm(:,:,:) 1206 # if defined key_zdfgls || defined key_zdftke 1207 en_25h(:,:,:) = en(:,:,:) 1208 #endif 1209 # if defined key_zdfgls 1210 mxln_25h(:,:,:) = mxln(:,:,:) 1211 #endif 1212 cnt_25h = 1 1213 IF (lwp) WRITE(numout,*) 'dia_wri_tide : After 25hr mean write, reset sum to current value and cnt_25h to one for overlapping average',cnt_25h 1214 1215 ENDIF ! cnt_25h .EQ. 25 .AND. MOD( kt, i_steps * 24) == 0 .AND. kt .NE. nn_it000 1216 1217 END SUBROUTINE dia_wri_tide 1218 !!====================================================================== 1219 1220 SUBROUTINE dia_wri_tide_init 1221 !!--------------------------------------------------------------------------- 1222 !! *** ROUTINE dia_wri_tide_init *** 1223 !! 1224 !! ** Purpose: Initialization of 25hour mean variables for detided output 1225 !! 1226 !! ** Method : Read namelist, allocate and assign initial values 1227 !! History 1228 !! 3.4 ! 03-13 (E. O'Dea) Routine to initialize dia_wri_tide 1229 !!--------------------------------------------------------------------------- 1230 !! 1231 INTEGER :: ierror ! local integer 1232 !! 1233 NAMELIST/nam_diatide/ ln_diatide 1234 !! 1235 !!---------------------------------------------------------------------- 1236 ! 1237 REWIND ( numnam ) ! Read Namelist nam_tiatide 1238 READ ( numnam, nam_diatide ) 1239 ! 1240 IF(lwp) THEN ! Control print 1241 WRITE(numout,*) 1242 WRITE(numout,*) 'dia_wri_tide_init : Output 25 hour Mean Diagnostics' 1243 WRITE(numout,*) '~~~~~~~~~~~~' 1244 WRITE(numout,*) ' Namelist nam_diatide : set 25hour mean outputs ' 1245 WRITE(numout,*) ' Switch for 25 hour mean diagnostics (T) or not (F) ln_diatide = ', ln_diatide 1246 ENDIF 1247 IF( .NOT. ln_diatide ) RETURN 1248 1249 ! ------------------- ! 1250 ! 1 - Allocate memory ! 1251 ! ------------------- ! 1252 ALLOCATE( tn_25h(jpi,jpj,jpk), STAT=ierror ) 1253 IF( ierror > 0 ) THEN 1254 CALL ctl_stop( 'dia_tide: unable to allocate tn_25h' ) ; RETURN 1255 ENDIF 1256 ALLOCATE( sn_25h(jpi,jpj,jpk), STAT=ierror ) 1257 IF( ierror > 0 ) THEN 1258 CALL ctl_stop( 'dia_tide: unable to allocate sn_25h' ) ; RETURN 1259 ENDIF 1260 ALLOCATE( insitu_t_25h(jpi,jpj,jpk), STAT=ierror ) 1261 IF( ierror > 0 ) THEN 1262 CALL ctl_stop( 'dia_tide: unable to allocate insitu_t_25h' ) ; RETURN 1263 ENDIF 1264 ALLOCATE( un_25h(jpi,jpj,jpk), STAT=ierror ) 1265 IF( ierror > 0 ) THEN 1266 CALL ctl_stop( 'dia_tide: unable to allocate un_25h' ) ; RETURN 1267 ENDIF 1268 ALLOCATE( vn_25h(jpi,jpj,jpk), STAT=ierror ) 1269 IF( ierror > 0 ) THEN 1270 CALL ctl_stop( 'dia_tide: unable to allocate vn_25h' ) ; RETURN 1271 ENDIF 1272 ALLOCATE( wn_25h(jpi,jpj,jpk), STAT=ierror ) 1273 IF( ierror > 0 ) THEN 1274 CALL ctl_stop( 'dia_tide: unable to allocate wn_25h' ) ; RETURN 1275 ENDIF 1276 ALLOCATE( avt_25h(jpi,jpj,jpk), STAT=ierror ) 1277 IF( ierror > 0 ) THEN 1278 CALL ctl_stop( 'dia_tide: unable to allocate avt_25h' ) ; RETURN 1279 ENDIF 1280 ALLOCATE( avm_25h(jpi,jpj,jpk), STAT=ierror ) 1281 IF( ierror > 0 ) THEN 1282 CALL ctl_stop( 'dia_tide: unable to allocate avm_25h' ) ; RETURN 1283 ENDIF 1284 # if defined key_zdfgls || defined key_zdftke 1285 ALLOCATE( en_25h(jpi,jpj,jpk), STAT=ierror ) 1286 IF( ierror > 0 ) THEN 1287 CALL ctl_stop( 'dia_tide: unable to allocate en_25h' ) ; RETURN 1288 ENDIF 1289 #endif 1290 # if defined key_zdfgls 1291 ALLOCATE( mxln_25h(jpi,jpj,jpk), STAT=ierror ) 1292 IF( ierror > 0 ) THEN 1293 CALL ctl_stop( 'dia_tide: unable to allocate mxln_25h' ) ; RETURN 1294 ENDIF 1295 #endif 1296 ALLOCATE( sshn_25h(jpi,jpj), STAT=ierror ) 1297 IF( ierror > 0 ) THEN 1298 CALL ctl_stop( 'dia_tide: unable to allocate sshn_25h' ) ; RETURN 1299 ENDIF 1300 ALLOCATE( hmld_kara_25h(jpi,jpj), STAT=ierror ) 1301 IF( ierror > 0 ) THEN 1302 CALL ctl_stop( 'dia_tide: unable to allocate hmld_kara_25h' ) ; RETURN 1303 ENDIF 1304 ! ------------------------- ! 1305 ! 2 - Assign Initial Values ! 1306 ! ------------------------- ! 1307 cnt_25h = 1 ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible) 1308 tn_25h(:,:,:) = tsb(:,:,:,jp_tem) 1309 sn_25h(:,:,:) = tsb(:,:,:,jp_sal) 1310 CALL theta2t 1311 insitu_t_25h(:,:,:) = insitu_t(:,:,:) 1312 sshn_25h(:,:) = sshb(:,:) 1313 ! hmld_kara_25h(:,:) = hmld_kara(:,:) 1314 un_25h(:,:,:) = ub(:,:,:) 1315 vn_25h(:,:,:) = vb(:,:,:) 1316 wn_25h(:,:,:) = wn(:,:,:) 1317 avt_25h(:,:,:) = avt(:,:,:) 1318 avm_25h(:,:,:) = avm(:,:,:) 1319 # if defined key_zdfgls || defined key_zdftke 1320 en_25h(:,:,:) = en(:,:,:) 1321 #endif 1322 # if defined key_zdfgls 1323 mxln_25h(:,:,:) = mxln(:,:,:) 1324 #endif 1325 #if defined key_lim3 || defined key_lim2 1326 CALL ctl_stop('STOP', 'dia_wri_tide not setup yet to do tidemean ice') 1327 #endif 1328 1329 ! -------------------------- ! 1330 ! 3 - Return to dia_wri_tide ! 1331 ! -------------------------- ! 1332 1333 1334 END SUBROUTINE dia_wri_tide_init 1335 1336 800 1337 END MODULE diawri -
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DOM/closea.F90
r7363 r7367 7 7 !! 8.5 ! 02-06 (E. Durand, G. Madec) F90 8 8 !! 9.0 ! 06-07 (G. Madec) add clo_rnf, clo_ups, clo_bat 9 !! NEMO 3.4 ! 03-12 (P.G. Fogli) sbc_clo bug fix & mpp reproducibility 9 10 !!---------------------------------------------------------------------- 10 11 … … 20 21 USE in_out_manager ! I/O manager 21 22 USE sbc_oce ! ocean surface boundary conditions 22 USE lib_mpp ! distributed memory computing library 23 USE lbclnk ! ??? 23 USE lib_fortran, ONLY: glob_sum, DDPDD 24 USE lbclnk ! lateral boundary condition - MPP exchanges 25 USE lib_mpp ! MPP library 26 USE timing 24 27 25 28 IMPLICIT NONE … … 85 88 SELECT CASE ( jp_cfg ) 86 89 ! ! ======================= 90 CASE ( 1 ) ! ORCA_R1 configuration 91 ! ! ======================= 92 ncsnr(1) = 1 ; ncstt(1) = 0 ! Caspian Sea 93 ncsi1(1) = 332 ; ncsj1(1) = 203 94 ncsi2(1) = 344 ; ncsj2(1) = 235 95 ncsir(1,1) = 1 ; ncsjr(1,1) = 1 96 ! 97 ! ! ======================= 87 98 CASE ( 2 ) ! ORCA_R2 configuration 88 99 ! ! ======================= … … 177 188 INTEGER, INTENT(in) :: kt ! ocean model time step 178 189 ! 179 INTEGER :: ji, jj, jc, jn ! dummy loop indices 180 REAL(wp) :: zze2 181 REAL(wp), DIMENSION (jpncs) :: zfwf 182 !!---------------------------------------------------------------------- 183 ! 190 INTEGER :: ji, jj, jc, jn ! dummy loop indices 191 REAL(wp), PARAMETER :: rsmall = 1.e-20_wp ! Closed sea correction epsilon 192 REAL(wp) :: zze2, ztmp, zcorr ! 193 COMPLEX(wp) :: ctmp 194 REAL(wp), DIMENSION(jpncs) :: zfwf ! 1D workspace 195 !!---------------------------------------------------------------------- 196 ! 197 IF( nn_timing == 1 ) CALL timing_start('sbc_clo') 184 198 ! !------------------! 185 199 IF( kt == nit000 ) THEN ! Initialisation ! … … 189 203 IF(lwp) WRITE(numout,*)'~~~~~~~' 190 204 191 ! Total surface of ocean 192 surf(jpncs+1) = SUM( e1t(:,:) * e2t(:,:) * tmask_i(:,:) ) 193 194 DO jc = 1, jpncs 195 surf(jc) =0.e0 196 DO jj = ncsj1(jc), ncsj2(jc) 197 DO ji = ncsi1(jc), ncsi2(jc) 198 surf(jc) = surf(jc) + e1t(ji,jj) * e2t(ji,jj) * tmask_i(ji,jj) ! surface of closed seas 205 surf(:) = 0.e0_wp 206 ! 207 surf(jpncs+1) = glob_sum( e1e2t(:,:) ) ! surface of the global ocean 208 ! 209 ! ! surface of closed seas 210 IF( lk_mpp_rep ) THEN ! MPP reproductible calculation 211 DO jc = 1, jpncs 212 ctmp = CMPLX( 0.e0, 0.e0, wp ) 213 DO jj = ncsj1(jc), ncsj2(jc) 214 DO ji = ncsi1(jc), ncsi2(jc) 215 ztmp = e1e2t(ji,jj) * tmask_i(ji,jj) 216 CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 217 END DO 199 218 END DO 200 END DO 201 END DO 202 IF( lk_mpp ) CALL mpp_sum ( surf, jpncs+1 ) ! mpp: sum over all the global domain 219 IF( lk_mpp ) CALL mpp_sum( ctmp ) 220 surf(jc) = REAL(ctmp,wp) 221 END DO 222 ELSE ! Standard calculation 223 DO jc = 1, jpncs 224 DO jj = ncsj1(jc), ncsj2(jc) 225 DO ji = ncsi1(jc), ncsi2(jc) 226 surf(jc) = surf(jc) + e1e2t(ji,jj) * tmask_i(ji,jj) ! surface of closed seas 227 END DO 228 END DO 229 END DO 230 IF( lk_mpp ) CALL mpp_sum ( surf, jpncs ) ! mpp: sum over all the global domain 231 ENDIF 203 232 204 233 IF(lwp) WRITE(numout,*)' Closed sea surfaces' … … 215 244 ! !--------------------! 216 245 ! ! update emp, emps ! 217 zfwf = 0.e0 !--------------------! 218 DO jc = 1, jpncs 219 DO jj = ncsj1(jc), ncsj2(jc) 220 DO ji = ncsi1(jc), ncsi2(jc) 221 zfwf(jc) = zfwf(jc) + e1t(ji,jj) * e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj) 222 END DO 223 END DO 224 END DO 225 IF( lk_mpp ) CALL mpp_sum ( zfwf(:) , jpncs ) ! mpp: sum over all the global domain 246 zfwf = 0.e0_wp !--------------------! 247 IF( lk_mpp_rep ) THEN ! MPP reproductible calculation 248 DO jc = 1, jpncs 249 ctmp = CMPLX( 0.e0, 0.e0, wp ) 250 DO jj = ncsj1(jc), ncsj2(jc) 251 DO ji = ncsi1(jc), ncsi2(jc) 252 ztmp = e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj) 253 CALL DDPDD( CMPLX( ztmp, 0.e0, wp ), ctmp ) 254 END DO 255 END DO 256 IF( lk_mpp ) CALL mpp_sum( ctmp ) 257 zfwf(jc) = REAL(ctmp,wp) 258 END DO 259 ELSE ! Standard calculation 260 DO jc = 1, jpncs 261 DO jj = ncsj1(jc), ncsj2(jc) 262 DO ji = ncsi1(jc), ncsi2(jc) 263 zfwf(jc) = zfwf(jc) + e1e2t(ji,jj) * ( emp(ji,jj)-rnf(ji,jj) ) * tmask_i(ji,jj) 264 END DO 265 END DO 266 END DO 267 IF( lk_mpp ) CALL mpp_sum ( zfwf(:) , jpncs ) ! mpp: sum over all the global domain 268 ENDIF 226 269 227 270 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! Black Sea case for ORCA_R2 configuration 228 zze2 = ( zfwf(3) + zfwf(4) ) / 2.271 zze2 = ( zfwf(3) + zfwf(4) ) * 0.5_wp 229 272 zfwf(3) = zze2 230 273 zfwf(4) = zze2 231 274 ENDIF 232 275 276 zcorr = 0._wp 277 233 278 DO jc = 1, jpncs 234 279 ! 235 IF( ncstt(jc) == 0 ) THEN 236 ! water/evap excess is shared by all open ocean 237 emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1) 238 emps(:,:) = emps(:,:) + zfwf(jc) / surf(jpncs+1) 239 ELSEIF( ncstt(jc) == 1 ) THEN 240 ! Excess water in open sea, at outflow location, excess evap shared 241 IF ( zfwf(jc) <= 0.e0 ) THEN 242 DO jn = 1, ncsnr(jc) 280 ! The following if avoids the redistribution of the round off 281 IF ( ABS(zfwf(jc) / surf(jpncs+1) ) > rsmall) THEN 282 ! 283 IF( ncstt(jc) == 0 ) THEN ! water/evap excess is shared by all open ocean 284 emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1) 285 emps(:,:) = emps(:,:) + zfwf(jc) / surf(jpncs+1) 286 ! accumulate closed seas correction 287 zcorr = zcorr + zfwf(jc) / surf(jpncs+1) 288 ! 289 ELSEIF( ncstt(jc) == 1 ) THEN ! Excess water in open sea, at outflow location, excess evap shared 290 IF ( zfwf(jc) <= 0.e0_wp ) THEN 291 DO jn = 1, ncsnr(jc) 292 ji = mi0(ncsir(jc,jn)) 293 jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean 294 IF ( ji > 1 .AND. ji < jpi & 295 .AND. jj > 1 .AND. jj < jpj ) THEN 296 emp (ji,jj) = emp (ji,jj) + zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) ) 297 emps(ji,jj) = emps(ji,jj) + zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) ) 298 ENDIF 299 END DO 300 ELSE 301 emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1) 302 emps(:,:) = emps(:,:) + zfwf(jc) / surf(jpncs+1) 303 ! accumulate closed seas correction 304 zcorr = zcorr + zfwf(jc) / surf(jpncs+1) 305 ENDIF 306 ELSEIF( ncstt(jc) == 2 ) THEN ! Excess e-p-r (either sign) goes to open ocean, at outflow location 307 DO jn = 1, ncsnr(jc) 243 308 ji = mi0(ncsir(jc,jn)) 244 309 jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean 245 IF ( ji > 1 .AND. ji < jpi & 246 .AND. jj > 1 .AND. jj < jpj ) THEN 247 emp (ji,jj) = emp (ji,jj) + zfwf(jc) / & 248 (FLOAT(ncsnr(jc)) * e1t(ji,jj) * e2t(ji,jj)) 249 emps(ji,jj) = emps(ji,jj) + zfwf(jc) / & 250 (FLOAT(ncsnr(jc)) * e1t(ji,jj) * e2t(ji,jj)) 251 END IF 252 END DO 253 ELSE 254 emp (:,:) = emp (:,:) + zfwf(jc) / surf(jpncs+1) 255 emps(:,:) = emps(:,:) + zfwf(jc) / surf(jpncs+1) 256 ENDIF 257 ELSEIF( ncstt(jc) == 2 ) THEN 258 ! Excess e-p+r (either sign) goes to open ocean, at outflow location 259 IF( ji > 1 .AND. ji < jpi & 260 .AND. jj > 1 .AND. jj < jpj ) THEN 261 DO jn = 1, ncsnr(jc) 262 ji = mi0(ncsir(jc,jn)) 263 jj = mj0(ncsjr(jc,jn)) ! Location of outflow in open ocean 264 emp (ji,jj) = emp (ji,jj) + zfwf(jc) & 265 / (FLOAT(ncsnr(jc)) * e1t(ji,jj) * e2t(ji,jj) ) 266 emps(ji,jj) = emps(ji,jj) + zfwf(jc) & 267 / (FLOAT(ncsnr(jc)) * e1t(ji,jj) * e2t(ji,jj) ) 268 END DO 310 IF( ji > 1 .AND. ji < jpi & 311 .AND. jj > 1 .AND. jj < jpj ) THEN 312 emp (ji,jj) = emp (ji,jj) + zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) ) 313 emps(ji,jj) = emps(ji,jj) + zfwf(jc) / ( REAL(ncsnr(jc)) * e1e2t(ji,jj) ) 314 ENDIF 315 END DO 269 316 ENDIF 270 ENDIF271 !272 DO jj = ncsj1(jc), ncsj2(jc)273 DO ji = ncsi1(jc), ncsi2(jc)274 emp (ji,jj) = emp(ji,jj) - zfwf(jc) / surf(jc)275 emps(ji,jj) = emps(ji,jj) - zfwf(jc) / surf(jc)276 END DO 277 END DO278 !317 ! 318 DO jj = ncsj1(jc), ncsj2(jc) 319 DO ji = ncsi1(jc), ncsi2(jc) 320 emp (ji,jj) = emp (ji,jj) - zfwf(jc) / surf(jc) 321 emps(ji,jj) = emps(ji,jj) - zfwf(jc) / surf(jc) 322 END DO 323 END DO 324 ! 325 END IF 279 326 END DO 280 ! 281 CALL lbc_lnk( emp , 'T', 1. ) 282 CALL lbc_lnk( emps, 'T', 1. ) 327 328 IF ( ABS(zcorr) > rsmall ) THEN ! remove the global correction from the closed seas 329 DO jc = 1, jpncs ! only if it is large enough 330 DO jj = ncsj1(jc), ncsj2(jc) 331 DO ji = ncsi1(jc), ncsi2(jc) 332 emp (ji,jj) = emp (ji,jj) - zcorr 333 emps(ji,jj) = emps(ji,jj) - zcorr 334 END DO 335 END DO 336 END DO 337 ENDIF 338 ! 339 emp (:,:) = emp (:,:) * tmask(:,:,1) 340 emps(:,:) = emps(:,:) * tmask(:,:,1) 341 ! 342 CALL lbc_lnk( emp , 'T', 1._wp ) 343 CALL lbc_lnk( emps, 'T', 1._wp ) 344 ! 345 IF( nn_timing == 1 ) CALL timing_stop('sbc_clo') 283 346 ! 284 347 END SUBROUTINE sbc_clo 285 286 348 349 287 350 SUBROUTINE clo_rnf( p_rnfmsk ) 288 351 !!--------------------------------------------------------------------- … … 308 371 ii = mi0( ncsir(jc,jn) ) 309 372 ij = mj0( ncsjr(jc,jn) ) 310 p_rnfmsk(ii,ij) = MAX( p_rnfmsk(ii,ij), 1.0 )373 p_rnfmsk(ii,ij) = MAX( p_rnfmsk(ii,ij), 1.0_wp ) 311 374 END DO 312 375 ENDIF … … 336 399 DO jj = ncsj1(jc), ncsj2(jc) 337 400 DO ji = ncsi1(jc), ncsi2(jc) 338 p_upsmsk(ji,jj) = 0.5 401 p_upsmsk(ji,jj) = 0.5_wp ! mixed upstream/centered scheme over closed seas 339 402 END DO 340 403 END DO … … 374 437 !!====================================================================== 375 438 END MODULE closea 439 -
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90
r7363 r7367 116 116 117 117 ! number of seconds since the beginning of current year/month/week/day at the middle of the time-step 118 nsec_year = nday_year * nsecd - ndt 05! 1 time step before the middle of the first time step119 nsec_month = nday * nsecd - ndt 05! because day will be called at the beginning of step120 nsec_week = idweek * nsecd - ndt 05121 nsec_day = nsecd - ndt 05118 nsec_year = nday_year * nsecd - ndt ! 1 time step before the middle of the first time step 119 nsec_month = nday * nsecd - ndt ! because day will be called at the beginning of step 120 nsec_week = idweek * nsecd - ndt 121 nsec_day = nsecd - ndt 122 122 123 123 ! control print … … 219 219 IF( ABS(adatrj - REAL(NINT(adatrj ),wp)) < zprec ) adatrj = REAL(NINT(adatrj ),wp) ! avoid truncation error 220 220 221 IF( nsec_day > nsecd ) THEN ! New day221 IF( nsec_day >= nsecd ) THEN ! New day 222 222 ! 223 223 nday = nday + 1 -
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r7363 r7367 52 52 REAL(wp), PUBLIC :: rdtmax !: maximum time step on tracers 53 53 REAL(wp), PUBLIC :: rdth !: depth variation of tracer step 54 INTEGER , PUBLIC :: nclosea !: =0 suppress closed sea/lake from the ORCA domain or not (=1)55 54 56 55 ! !!! associated variables … … 125 124 LOGICAL, PUBLIC :: ln_zps = .FALSE. !: z-coordinate - partial step 126 125 LOGICAL, PUBLIC :: ln_sco = .FALSE. !: s-coordinate or hybrid z-s coordinate 126 LOGICAL, PUBLIC :: ln_read_zenv = .FALSE. !: Whether to read zenv or calculate it 127 127 128 128 !! All coordinates … … 173 173 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: scosrf, scobot !: ocean surface and bottom topographies 174 174 ! ! (if deviating from coordinate surfaces in HYBRID) 175 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rx1 !: maximum grid stiffness ratio 175 176 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hifv , hiff !: interface depth between stretching at V--F 176 177 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hift , hifu !: and quasi-uniform spacing T--U points (m) 178 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: zenv !: Envelope Batymetry, calcualted or read in 177 179 178 180 !!---------------------------------------------------------------------- … … 295 297 & scosrf(jpi,jpj) , scobot(jpi,jpj) , & 296 298 & hifv (jpi,jpj) , hiff (jpi,jpj) , & 297 & hift (jpi,jpj) , hifu (jpi,jpj) , STAT=ierr(8) )299 & hift (jpi,jpj) , hifu (jpi,jpj) , rx1(jpi,jpj), STAT=ierr(8) ) 298 300 299 301 ALLOCATE( mbathy(jpi,jpj) , bathy(jpi,jpj) , & 300 302 & tmask_i(jpi,jpj) , bmask(jpi,jpj) , & 301 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , STAT=ierr(9) )303 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv(jpi,jpj) , zenv(jpi,jpj), STAT=ierr(9) ) 302 304 303 305 ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk), & -
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r7363 r7367 36 36 USE dyncor_c1d ! Coriolis term (c1d case) (cor_c1d routine) 37 37 USE timing ! Timing 38 USE lbclnk 38 39 39 40 IMPLICIT NONE … … 84 85 CALL dom_zgr ! Vertical mesh and bathymetry 85 86 CALL dom_msk ! Masks 87 IF( ln_sco ) CALL dom_stiff ! Maximum stiffness ratio/hydrostatic consistency 86 88 IF( lk_vvl ) CALL dom_vvl ! Vertical variable mesh 87 89 ! … … 123 125 !!---------------------------------------------------------------------- 124 126 USE ioipsl 127 NAMELIST/namrun/ ln_NOOS 125 128 NAMELIST/namrun/ nn_no , cn_exp , cn_ocerst_in, cn_ocerst_out, ln_rstart , nn_rstctl, & 129 & nn_stocklist, & 126 130 & nn_it000, nn_itend , nn_date0 , nn_leapy , nn_istate , nn_stock , & 127 & nn_write, ln_dimgnnn, ln_mskland , ln_clobber , nn_chunksz 131 & nn_write, ln_dimgnnn, ln_mskland , ln_clobber , nn_chunksz, & 132 & ln_diafoam, nn_diafoam, ln_depwri 128 133 NAMELIST/namdom/ nn_bathy , rn_e3zps_min, rn_e3zps_rat, nn_msh , rn_hmin, & 129 134 & nn_acc , rn_atfp , rn_rdt , rn_rdtmin , & 130 135 & rn_rdtmax, rn_rdth , nn_baro , nn_closea 131 136 NAMELIST/namcla/ nn_cla 137 NAMELIST/namrun/ ln_rstdate 132 138 #if defined key_netcdf4 133 139 NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip 134 140 #endif 135 !!---------------------------------------------------------------------- 136 141 NAMELIST/namrun/ ln_diatide 142 !!---------------------------------------------------------------------- 143 NAMELIST/namrun/ ln_diatmb 144 145 146 NAMELIST/namrun/ cn_rst_dir ! moved here to allow merge with CO5 branches (ln_NOOS) 137 147 REWIND( numnam ) ! Namelist namrun : parameters of the run 138 148 READ ( numnam, namrun ) … … 152 162 WRITE(numout,*) ' leap year calendar (0/1) nn_leapy = ', nn_leapy 153 163 WRITE(numout,*) ' initial state output nn_istate = ', nn_istate 154 WRITE(numout,*) ' frequency of restart file nn_stock = ', nn_stock 164 IF ( ALL( nn_stocklist == 0 ) ) THEN 165 WRITE(numout,*) ' frequency of restart file nn_stock = ', nn_stock 166 ELSE 167 WRITE(numout,*) ' list of restart times nn_stocklist = ', nn_stocklist(1:10) 168 ENDIF 155 169 WRITE(numout,*) ' frequency of output file nn_write = ', nn_write 156 170 WRITE(numout,*) ' multi file dimgout ln_dimgnnn = ', ln_dimgnnn 171 WRITE(numout,*) ' use date in restart name ln_rstdate = ', ln_rstdate 157 172 WRITE(numout,*) ' mask land points ln_mskland = ', ln_mskland 173 WRITE(numout,*) ' NOOS transect diagnostics ln_NOOS = ', ln_NOOS 158 174 WRITE(numout,*) ' overwrite an existing file ln_clobber = ', ln_clobber 175 WRITE(numout,*) ' restart directory cn_rst_dir = ', cn_rst_dir 159 176 WRITE(numout,*) ' NetCDF chunksize (bytes) nn_chunksz = ', nn_chunksz 177 WRITE(numout,*) ' Met Office FOAM diagnostics ln_diafoam = ', ln_diafoam 178 WRITE(numout,*) ' FOAM diagnostic choices nn_diafoam = ', nn_diafoam 179 WRITE(numout,*) ' depths file output logical ln_depwri = ', ln_depwri 160 180 ENDIF 161 181 … … 169 189 ninist = nn_istate 170 190 nstock = nn_stock 191 nstock_list = nn_stocklist 171 192 nwrite = nn_write 172 193 … … 238 259 rdtmax = rn_rdtmin 239 260 rdth = rn_rdth 240 nclosea = nn_closea241 261 242 262 REWIND( numnam ) ! Namelist cross land advection … … 274 294 275 295 296 !!====================================================================== 276 297 SUBROUTINE dom_ctl 277 298 !!---------------------------------------------------------------------- … … 323 344 END SUBROUTINE dom_ctl 324 345 346 SUBROUTINE dom_stiff 347 !!---------------------------------------------------------------------- 348 !! *** ROUTINE dom_stiff *** 349 !! 350 !! ** Purpose : Diagnose maximum grid stiffness/hydrostatic consistency 351 !! 352 !! ** Method : Compute Haney (1991) hydrostatic condition ratio 353 !! Save the maximum in the vertical direction 354 !! (this number is only relevant in s-coordinates) 355 !! 356 !! Haney, R. L., 1991: On the pressure gradient force 357 !! over steep topography in sigma coordinate ocean models. 358 !! J. Phys. Oceanogr., 21, 610???619. 359 !!---------------------------------------------------------------------- 360 INTEGER :: ji, jj, jk 361 REAL(wp) :: zrxmax 362 REAL(wp), DIMENSION(4) :: zr1 363 !!---------------------------------------------------------------------- 364 rx1(:,:) = 0.e0 365 zrxmax = 0.e0 366 zr1(:) = 0.e0 367 368 DO ji = 2, jpim1 369 DO jj = 2, jpjm1 370 DO jk = 1, jpkm1 371 zr1(1) = umask(ji-1,jj ,jk) *abs( (gdepw(ji ,jj ,jk )-gdepw(ji-1,jj ,jk ) & 372 & +gdepw(ji ,jj ,jk+1)-gdepw(ji-1,jj ,jk+1)) & 373 & /(gdepw(ji ,jj ,jk )+gdepw(ji-1,jj ,jk ) & 374 & -gdepw(ji ,jj ,jk+1)-gdepw(ji-1,jj ,jk+1) + rsmall) ) 375 zr1(2) = umask(ji ,jj ,jk) *abs( (gdepw(ji+1,jj ,jk )-gdepw(ji ,jj ,jk ) & 376 & +gdepw(ji+1,jj ,jk+1)-gdepw(ji ,jj ,jk+1)) & 377 & /(gdepw(ji+1,jj ,jk )+gdepw(ji ,jj ,jk ) & 378 & -gdepw(ji+1,jj ,jk+1)-gdepw(ji ,jj ,jk+1) + rsmall) ) 379 zr1(3) = vmask(ji ,jj ,jk) *abs( (gdepw(ji ,jj+1,jk )-gdepw(ji ,jj ,jk ) & 380 & +gdepw(ji ,jj+1,jk+1)-gdepw(ji ,jj ,jk+1)) & 381 & /(gdepw(ji ,jj+1,jk )+gdepw(ji ,jj ,jk ) & 382 & -gdepw(ji ,jj+1,jk+1)-gdepw(ji ,jj ,jk+1) + rsmall) ) 383 zr1(4) = vmask(ji ,jj-1,jk) *abs( (gdepw(ji ,jj ,jk )-gdepw(ji ,jj-1,jk ) & 384 & +gdepw(ji ,jj ,jk+1)-gdepw(ji ,jj-1,jk+1)) & 385 & /(gdepw(ji ,jj ,jk )+gdepw(ji ,jj-1,jk ) & 386 & -gdepw(ji, jj ,jk+1)-gdepw(ji ,jj-1,jk+1) + rsmall) ) 387 zrxmax = MAXVAL(zr1(1:4)) 388 rx1(ji,jj) = MAX(rx1(ji,jj), zrxmax) 389 END DO 390 END DO 391 END DO 392 393 CALL lbc_lnk( rx1, 'T', 1. ) 394 395 zrxmax = MAXVAL(rx1) 396 397 IF( lk_mpp ) CALL mpp_max( zrxmax ) ! max over the global domain 398 399 IF(lwp) THEN 400 WRITE(numout,*) 401 WRITE(numout,*) 'dom_stiff : maximum grid stiffness ratio: ', zrxmax 402 WRITE(numout,*) '~~~~~~~~~' 403 ENDIF 404 405 END SUBROUTINE dom_stiff 406 325 407 !!====================================================================== 326 408 END MODULE domain -
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
r7363 r7367 172 172 173 173 IF( ln_sco ) THEN ! s-coordinate 174 CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt ) ! ! depth175 CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu ) 174 CALL iom_rstput( 0, 0, inum4, 'hbatt', hbatt ) 175 CALL iom_rstput( 0, 0, inum4, 'hbatu', hbatu ) 176 176 CALL iom_rstput( 0, 0, inum4, 'hbatv', hbatv ) 177 177 CALL iom_rstput( 0, 0, inum4, 'hbatf', hbatf ) … … 187 187 CALL iom_rstput( 0, 0, inum4, 'e3v', e3v ) 188 188 CALL iom_rstput( 0, 0, inum4, 'e3w', e3w ) 189 ! 190 CALL iom_rstput( 0, 0, inum4, 'gdept_0' , gdept_0 ) ! ! stretched system 191 CALL iom_rstput( 0, 0, inum4, 'gdepw_0' , gdepw_0 ) 189 CALL iom_rstput( 0, 0, inum4, 'rx1', rx1 ) ! ! Max. grid stiffness ratio 190 ! 191 CALL iom_rstput( 0, 0, inum4, 'gdept' , gdept ) ! ! stretched system 192 CALL iom_rstput( 0, 0, inum4, 'gdepw' , gdepw ) 192 193 ENDIF 193 194 -
branches/UKMO/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r7363 r7367 15 15 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 16 16 !! 3.3 ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level 17 !! 3.4 ! 2012-08 (J. Siddorn) added Siddorn adn Furner stretching functio 17 18 !!---------------------------------------------------------------------- 18 19 … … 28 29 !! zgr_sco : s-coordinate 29 30 !! fssig : sigma coordinate non-dimensional function 31 !! fgamma : Siddorn and Furner stretching function 30 32 !! dfssig : derivative of the sigma coordinate function !!gm (currently missing!) 31 33 !!--------------------------------------------------------------------- … … 47 49 48 50 ! !!* Namelist namzgr_sco * 51 LOGICAL :: ln_s_sh94 = .false. ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T) 52 LOGICAL :: ln_s_sf12 = .true. ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T) 53 LOGICAL :: ln_sigcrit = .false. ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch 54 ! 49 55 REAL(wp) :: rn_sbot_min = 300._wp ! minimum depth of s-bottom surface (>0) (m) 50 56 REAL(wp) :: rn_sbot_max = 5250._wp ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 57 REAL(wp) :: rn_rmax = 0.15_wp ! maximum cut-off r-value allowed (0<rn_rmax<1) 58 REAL(wp) :: rn_hc = 150._wp ! Critical depth for transition from sigma to stretched coordinates 59 ! Song and Haidvogel 1994 stretching parameters 51 60 REAL(wp) :: rn_theta = 6.00_wp ! surface control parameter (0<=rn_theta<=20) 52 61 REAL(wp) :: rn_thetb = 0.75_wp ! bottom control parameter (0<=rn_thetb<= 1) 53 REAL(wp) :: rn_rmax = 0.15_wp ! maximum cut-off r-value allowed (0<rn_rmax<1) 54 LOGICAL :: ln_s_sigma = .false. ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T) 55 REAL(wp) :: rn_bb = 0.80_wp ! stretching parameter for song and haidvogel stretching 62 REAL(wp) :: rn_bb = 0.80_wp ! stretching parameter 56 63 ! ! ( rn_bb=0; top only, rn_bb =1; top and bottom) 57 REAL(wp) :: rn_hc = 150._wp ! Critical depth for s-sigma coordinates 64 ! Siddorn and Furner stretching parameters 65 REAL(wp) :: rn_alpha = 4.4_wp ! control parameter ( > 1 stretch towards surface, < 1 towards seabed) 66 REAL(wp) :: rn_efold = 0.0_wp ! efold length scale for transition to stretched coord 67 REAL(wp) :: rn_zs = 1.0_wp ! depth of surface grid box 68 ! bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b 69 REAL(wp) :: rn_zb_a = 0.024_wp ! bathymetry scaling factor for calculating Zb 70 REAL(wp) :: rn_zb_b = -0.2_wp ! offset for calculating Zb 58 71 59 72 !! * Substitutions … … 86 99 INTEGER :: ioptio = 0 ! temporary integer 87 100 ! 88 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco 101 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_read_zenv 89 102 !!---------------------------------------------------------------------- 90 103 ! … … 102 115 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 103 116 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 117 WRITE(numout,*) ' Read Zenv from Bathy T/F ln_read_zenv = ', ln_read_zenv 104 118 ENDIF 105 119 … … 243 257 END DO 244 258 ELSE ! Madec & Imbard 1996 function 259 # if key_levels==1 260 !Hard wire a deep and shallow level 261 !NOTE this configuration is for use with NEMOVAR, 262 !it is not set-up for NEMO 263 CALL ctl_warn("Single level model, depth of first layer set to 1cm."//& 264 & "\nThis configuration is designed to be used with NEMOVAR only") 265 gdepw_0(1)=0 266 gdept_0(1)=0.01 267 gdepw_0(2)=7000 268 gdept_0(2)=14000 269 e3w_0(:)=7000 270 e3t_0(1)=6999.99 271 e3t_0(2)=7000 272 # else 245 273 IF( .NOT. ldbletanh ) THEN 246 274 DO jk = 1, jpk … … 267 295 END DO 268 296 ENDIF 297 # endif 269 298 gdepw_0(1) = 0._wp ! force first w-level to be exactly at zero 270 299 ENDIF … … 422 451 CALL iom_close( inum ) 423 452 mbathy(:,:) = INT( bathy(:,:) ) 424 ! ! =====================453 ! 425 454 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 426 ! ! =====================455 ! 427 456 IF( nn_cla == 0 ) THEN 428 457 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open … … 453 482 CALL iom_open ( 'bathy_meter.nc', inum ) 454 483 CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy ) 484 IF ( ln_read_zenv ) THEN ! Whether we should read zenv or not 485 CALL iom_get ( inum, jpdom_data, 'zenv', zenv ) 486 ENDIF 455 487 CALL iom_close( inum ) 456 ! ! =====================488 ! 457 489 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 458 ! ! =====================490 ! 459 491 IF( nn_cla == 0 ) THEN 460 492 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open … … 489 521 ENDIF 490 522 ! 491 ! ! =========================== ! 492 IF( nclosea == 0 ) THEN ! NO closed seas or lakes ! 493 DO jl = 1, jpncs ! =========================== ! 494 DO jj = ncsj1(jl), ncsj2(jl) 495 DO ji = ncsi1(jl), ncsi2(jl) 496 mbathy(ji,jj) = 0 ! suppress closed seas and lakes from bathymetry 497 bathy (ji,jj) = 0._wp 498 END DO 499 END DO 500 END DO 501 ENDIF 502 ! 503 ! ! =========================== ! 504 ! ! set a minimum depth ! 505 ! ! =========================== ! 506 IF ( .not. ln_sco ) THEN 523 IF( nn_closea == 0 ) CALL clo_bat( bathy, mbathy ) !== NO closed seas or lakes ==! 524 ! 525 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 507 526 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level 508 527 ELSE ; ik = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth … … 1047 1066 END SUBROUTINE zgr_zps 1048 1067 1049 1050 FUNCTION fssig( pk ) RESULT( pf )1051 !!----------------------------------------------------------------------1052 !! *** ROUTINE eos_init ***1053 !!1054 !! ** Purpose : provide the analytical function in s-coordinate1055 !!1056 !! ** Method : the function provide the non-dimensional position of1057 !! T and W (i.e. between 0 and 1)1058 !! T-points at integer values (between 1 and jpk)1059 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)1060 !!----------------------------------------------------------------------1061 REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate1062 REAL(wp) :: pf ! sigma value1063 !!----------------------------------------------------------------------1064 !1065 pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) &1066 & - TANH( rn_thetb * rn_theta ) ) &1067 & * ( COSH( rn_theta ) &1068 & + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &1069 & / ( 2._wp * SINH( rn_theta ) )1070 !1071 END FUNCTION fssig1072 1073 1074 FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )1075 !!----------------------------------------------------------------------1076 !! *** ROUTINE eos_init ***1077 !!1078 !! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate1079 !!1080 !! ** Method : the function provides the non-dimensional position of1081 !! T and W (i.e. between 0 and 1)1082 !! T-points at integer values (between 1 and jpk)1083 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)1084 !!----------------------------------------------------------------------1085 REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate1086 REAL(wp), INTENT(in) :: pbb ! Stretching coefficient1087 REAL(wp) :: pf1 ! sigma value1088 !!----------------------------------------------------------------------1089 !1090 IF ( rn_theta == 0 ) then ! uniform sigma1091 pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )1092 ELSE ! stretched sigma1093 pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) &1094 & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta ) ) &1095 & / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )1096 ENDIF1097 !1098 END FUNCTION fssig11099 1100 1101 1068 SUBROUTINE zgr_sco 1102 1069 !!---------------------------------------------------------------------- … … 1123 1090 !! esigt(k) = fsdsig(k ) 1124 1091 !! esigw(k) = fsdsig(k-0.5) 1125 !! Th is routine is given as an example, it mustbe modified1126 !! following the user s desiderata. nevertheless, the output as1092 !! Three options for stretching are give, and they can be modified 1093 !! following the users requirements. Nevertheless, the output as 1127 1094 !! well as the way to compute the model levels and scale factors 1128 !! must be respected in order to insure second order a !!uracy1095 !! must be respected in order to insure second order accuracy 1129 1096 !! schemes. 1130 1097 !! 1131 !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 1098 !! The three methods for stretching available are: 1099 !! 1100 !! s_sh94 (Song and Haidvogel 1994) 1101 !! a sinh/tanh function that allows sigma and stretched sigma 1102 !! 1103 !! s_sf12 (Siddorn and Furner 2012?) 1104 !! allows the maintenance of fixed surface and or 1105 !! bottom cell resolutions (cf. geopotential coordinates) 1106 !! within an analytically derived stretched S-coordinate framework. 1107 !! 1108 !! s_tanh (Madec et al 1996) 1109 !! a cosh/tanh function that gives stretched coordinates 1110 !! 1132 1111 !!---------------------------------------------------------------------- 1133 1112 ! 1134 1113 INTEGER :: ji, jj, jk, jl ! dummy loop argument 1135 1114 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers 1136 REAL(wp) :: z coeft, zcoefw, zrmax, ztaper ! temporary scalars1137 ! 1138 REAL(wp), POINTER, DIMENSION(:,: ) :: z env, ztmp, zmsk, zri, zrj, zhbat1115 REAL(wp) :: zrmax, ztaper ! temporary scalars 1116 ! 1117 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmp, zmsk, zri, zrj, zhbat 1139 1118 REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 1140 1119 REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 1141 1120 1142 NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc 1143 !!---------------------------------------------------------------------- 1121 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 1122 rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 1123 !!---------------------------------------------------------------------- 1144 1124 ! 1145 1125 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1146 1126 ! 1147 CALL wrk_alloc( jpi, jpj, z env, ztmp, zmsk, zri, zrj, zhbat )1127 CALL wrk_alloc( jpi, jpj, ztmp, zmsk, zri, zrj, zhbat ) 1148 1128 CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3 ) 1149 1129 CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) … … 1157 1137 WRITE(numout,*) '~~~~~~~~~~~' 1158 1138 WRITE(numout,*) ' Namelist namzgr_sco' 1159 WRITE(numout,*) ' sigma-stretching coeffs ' 1160 WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ' ,rn_sbot_max 1161 WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ' ,rn_sbot_min 1162 WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ', rn_theta 1163 WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ', rn_thetb 1164 WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ', rn_rmax 1165 WRITE(numout,*) ' Hybrid s-sigma-coordinate ln_s_sigma = ', ln_s_sigma 1166 WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ', rn_bb 1167 WRITE(numout,*) ' Critical depth rn_hc = ', rn_hc 1168 ENDIF 1169 1170 gsigw3 = 0._wp ; gsigt3 = 0._wp ; gsi3w3 = 0._wp 1171 esigt3 = 0._wp ; esigw3 = 0._wp 1172 esigtu3 = 0._wp ; esigtv3 = 0._wp ; esigtf3 = 0._wp 1173 esigwu3 = 0._wp ; esigwv3 = 0._wp 1139 WRITE(numout,*) ' stretching coeffs ' 1140 WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max 1141 WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min 1142 WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc 1143 WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax 1144 WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94 1145 WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients' 1146 WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta 1147 WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb 1148 WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb 1149 WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12 1150 WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients' 1151 WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha 1152 WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold 1153 WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs 1154 WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a 1155 WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b 1156 WRITE(numout,*) ' Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b' 1157 ENDIF 1174 1158 1175 1159 hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate … … 1190 1174 ! ! ============================= 1191 1175 ! use r-value to create hybrid coordinates 1176 1177 ! Smooth the bathymetry (if required) 1178 scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea) 1179 scobot(:,:) = bathy(:,:) ! ocean bottom depth 1180 IF( ln_read_zenv) THEN 1181 WRITE(numout,*) ' Zenv is not calculated but read from Bathy File ln_read_zenv = ', ln_read_zenv 1182 ELSE 1183 IF ( jpnij .gt.1) CALL ctl_stop( ' zgr_zps : ln_read_zenv=false and jpnij > 1, Calculating zenv on more than one Proc is not safe, calculate on one proc only ' ) 1184 1192 1185 DO jj = 1, jpj 1193 1186 DO ji = 1, jpi … … 1196 1189 END DO 1197 1190 ! 1198 ! Smooth the bathymetry (if required)1199 scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)1200 scobot(:,:) = bathy(:,:) ! ocean bottom depth1201 1191 ! 1202 1192 jl = 0 … … 1270 1260 ! 1271 1261 ! ! envelop bathymetry saved in hbatt 1262 ENDIF ! End of IF block for reading from a file or calculating zenv 1272 1263 hbatt(:,:) = zenv(:,:) 1273 1264 IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN … … 1365 1356 ! non-dimensional "sigma" for model level depth at w- and t-levels 1366 1357 1367 IF( ln_s_sigma ) THEN ! Song and Haidvogel style stretched sigma for depths 1368 ! ! below rn_hc, with uniform sigma in shallower waters 1369 DO ji = 1, jpi 1370 DO jj = 1, jpj 1371 1372 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 1373 DO jk = 1, jpk 1374 gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 1375 gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) 1376 END DO 1377 ELSE ! shallow water, uniform sigma 1378 DO jk = 1, jpk 1379 gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp) 1380 gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 1381 END DO 1382 ENDIF 1383 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk) 1384 ! 1385 DO jk = 1, jpkm1 1386 esigt3(ji,jj,jk ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 1387 esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 1388 END DO 1389 esigw3(ji,jj,1 ) = 2._wp * ( gsigt3(ji,jj,1 ) - gsigw3(ji,jj,1 ) ) 1390 esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) ) 1391 ! 1392 ! Coefficients for vertical depth as the sum of e3w scale factors 1393 gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1) 1394 DO jk = 2, jpk 1395 gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 1396 END DO 1397 ! 1398 DO jk = 1, jpk 1399 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1400 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1401 gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 1402 gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 1403 gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 1404 END DO 1405 ! 1406 END DO ! for all jj's 1407 END DO ! for all ji's 1408 1409 DO ji = 1, jpim1 1410 DO jj = 1, jpjm1 1411 DO jk = 1, jpk 1412 esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) & 1413 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1414 esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) & 1415 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1416 esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) & 1417 & + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) & 1418 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1419 esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) & 1420 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1421 esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) & 1422 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1423 ! 1424 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1425 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1426 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1427 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1428 ! 1429 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1430 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1431 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1432 END DO 1433 END DO 1434 END DO 1435 1436 CALL lbc_lnk( e3t , 'T', 1._wp ) 1437 CALL lbc_lnk( e3u , 'U', 1._wp ) 1438 CALL lbc_lnk( e3v , 'V', 1._wp ) 1439 CALL lbc_lnk( e3f , 'F', 1._wp ) 1440 CALL lbc_lnk( e3w , 'W', 1._wp ) 1441 CALL lbc_lnk( e3uw, 'U', 1._wp ) 1442 CALL lbc_lnk( e3vw, 'V', 1._wp ) 1443 1444 ! 1445 ELSE ! not ln_s_sigma 1446 ! 1447 DO jk = 1, jpk 1448 gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 1449 gsigt(jk) = -fssig( REAL(jk,wp) ) 1450 END DO 1451 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk) 1452 ! 1453 ! Coefficients for vertical scale factors at w-, t- levels 1454 !!gm bug : define it from analytical function, not like juste bellow.... 1455 !!gm or betteroffer the 2 possibilities.... 1456 DO jk = 1, jpkm1 1457 esigt(jk ) = gsigw(jk+1) - gsigw(jk) 1458 esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 1459 END DO 1460 esigw( 1 ) = 2._wp * ( gsigt(1 ) - gsigw(1 ) ) 1461 esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) ) 1462 1463 !!gm original form 1464 !!org DO jk = 1, jpk 1465 !!org esigt(jk)=fsdsig( FLOAT(jk) ) 1466 !!org esigw(jk)=fsdsig( FLOAT(jk)-0.5 ) 1467 !!org END DO 1468 !!gm 1469 ! 1470 ! Coefficients for vertical depth as the sum of e3w scale factors 1471 gsi3w(1) = 0.5_wp * esigw(1) 1472 DO jk = 2, jpk 1473 gsi3w(jk) = gsi3w(jk-1) + esigw(jk) 1474 END DO 1475 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 1476 DO jk = 1, jpk 1477 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1478 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1479 gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft ) 1480 gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw ) 1481 gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft ) 1482 END DO 1483 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 1484 DO jj = 1, jpj 1485 DO ji = 1, jpi 1486 DO jk = 1, jpk 1487 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1488 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1489 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1490 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 1491 ! 1492 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1493 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1494 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1495 END DO 1496 END DO 1497 END DO 1498 ! 1499 ENDIF ! ln_s_sigma 1500 1501 1358 1359 !======================================================================== 1360 ! Song and Haidvogel 1994 (ln_s_sh94=T) 1361 ! Siddorn and Furner 2012 (ln_sf12=T) 1362 ! or tanh function (both false) 1363 !======================================================================== 1364 IF ( ln_s_sh94 ) THEN 1365 CALL s_sh94() 1366 ELSE IF ( ln_s_sf12 ) THEN 1367 CALL s_sf12() 1368 ELSE 1369 CALL s_tanh() 1370 ENDIF 1371 1372 CALL lbc_lnk( e3t , 'T', 1._wp ) 1373 CALL lbc_lnk( e3u , 'U', 1._wp ) 1374 CALL lbc_lnk( e3v , 'V', 1._wp ) 1375 CALL lbc_lnk( e3f , 'F', 1._wp ) 1376 CALL lbc_lnk( e3w , 'W', 1._wp ) 1377 CALL lbc_lnk( e3uw, 'U', 1._wp ) 1378 CALL lbc_lnk( e3vw, 'V', 1._wp ) 1379 1380 fsdepw(:,:,:) = gdepw (:,:,:) 1381 fsde3w(:,:,:) = gdep3w(:,:,:) 1502 1382 ! 1503 1383 where (e3t (:,:,:).eq.0.0) e3t(:,:,:) = 1.0 … … 1557 1437 & ' w ', MAXVAL( fse3w (:,:,:) ) 1558 1438 ENDIF 1559 ! 1439 ! END DO 1560 1440 IF(lwp) THEN ! selected vertical profiles 1561 1441 WRITE(numout,*) … … 1587 1467 ENDIF 1588 1468 1589 !!gm bug? no more necessary? if ! defined key_helsinki 1590 DO jk = 1, jpk 1469 !================================================================================ 1470 ! check the coordinate makes sense 1471 !================================================================================ 1472 DO ji = 1, jpi 1591 1473 DO jj = 1, jpj 1592 DO ji = 1, jpi 1593 IF( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 1594 WRITE(ctmp1,*) 'zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1595 CALL ctl_stop( ctmp1 ) 1596 ENDIF 1597 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 1598 WRITE(ctmp1,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1599 CALL ctl_stop( ctmp1 ) 1600 ENDIF 1601 END DO 1602 END DO 1603 END DO 1604 !!gm bug #endif 1605 ! 1606 CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat ) 1474 1475 IF( hbatt(ji,jj) > 0._wp) THEN 1476 DO jk = 1, mbathy(ji,jj) 1477 ! check coordinate is monotonically increasing 1478 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 1479 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1480 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1481 WRITE(numout,*) 'e3w',fse3w(ji,jj,:) 1482 WRITE(numout,*) 'e3t',fse3t(ji,jj,:) 1483 CALL ctl_stop( ctmp1 ) 1484 ENDIF 1485 ! and check it has never gone negative 1486 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 1487 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1488 WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1489 WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 1490 WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 1491 CALL ctl_stop( ctmp1 ) 1492 ENDIF 1493 ! and check it never exceeds the total depth 1494 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN 1495 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 1496 WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 1497 WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 1498 CALL ctl_stop( ctmp1 ) 1499 ENDIF 1500 END DO 1501 1502 DO jk = 1, mbathy(ji,jj)-1 1503 ! and check it never exceeds the total depth 1504 IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN 1505 WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 1506 WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 1507 WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 1508 CALL ctl_stop( ctmp1 ) 1509 ENDIF 1510 END DO 1511 1512 ENDIF 1513 1514 END DO 1515 END DO 1516 ! 1517 CALL wrk_dealloc( jpi, jpj, ztmp, zmsk, zri, zrj, zhbat ) 1607 1518 CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3 ) 1608 1519 CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) … … 1612 1523 END SUBROUTINE zgr_sco 1613 1524 1525 !!====================================================================== 1526 SUBROUTINE s_sh94() 1527 1528 !!---------------------------------------------------------------------- 1529 !! *** ROUTINE s_sh94 *** 1530 !! 1531 !! ** Purpose : stretch the s-coordinate system 1532 !! 1533 !! ** Method : s-coordinate stretch using the Song and Haidvogel 1994 1534 !! mixed S/sigma coordinate 1535 !! 1536 !! Reference : Song and Haidvogel 1994. 1537 !!---------------------------------------------------------------------- 1538 ! 1539 INTEGER :: ji, jj, jk ! dummy loop argument 1540 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 1541 ! 1542 REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 1543 REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 1544 1545 CALL wrk_alloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3 ) 1546 CALL wrk_alloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 1547 1548 gsigw3 = 0._wp ; gsigt3 = 0._wp ; gsi3w3 = 0._wp 1549 esigt3 = 0._wp ; esigw3 = 0._wp 1550 esigtu3 = 0._wp ; esigtv3 = 0._wp ; esigtf3 = 0._wp 1551 esigwu3 = 0._wp ; esigwv3 = 0._wp 1552 1553 DO ji = 1, jpi 1554 DO jj = 1, jpj 1555 1556 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 1557 DO jk = 1, jpk 1558 gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 1559 gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) 1560 END DO 1561 ELSE ! shallow water, uniform sigma 1562 DO jk = 1, jpk 1563 gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp) 1564 gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 1565 END DO 1566 ENDIF 1567 ! 1568 DO jk = 1, jpkm1 1569 esigt3(ji,jj,jk ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 1570 esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 1571 END DO 1572 esigw3(ji,jj,1 ) = 2._wp * ( gsigt3(ji,jj,1 ) - gsigw3(ji,jj,1 ) ) 1573 esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) ) 1574 ! 1575 ! Coefficients for vertical depth as the sum of e3w scale factors 1576 gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1) 1577 DO jk = 2, jpk 1578 gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 1579 END DO 1580 ! 1581 DO jk = 1, jpk 1582 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1583 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1584 gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 1585 gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 1586 gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 1587 END DO 1588 ! 1589 END DO ! for all jj's 1590 END DO ! for all ji's 1591 1592 DO ji = 1, jpim1 1593 DO jj = 1, jpjm1 1594 DO jk = 1, jpk 1595 esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) & 1596 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1597 esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) & 1598 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1599 esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) & 1600 & + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) & 1601 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1602 esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) & 1603 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1604 esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) & 1605 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1606 ! 1607 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1608 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1609 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1610 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1611 ! 1612 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1613 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1614 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 1615 END DO 1616 END DO 1617 END DO 1618 1619 CALL wrk_dealloc( jpi, jpj, jpk, gsigw3, gsigt3, gsi3w3 ) 1620 CALL wrk_dealloc( jpi, jpj, jpk, esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 ) 1621 1622 END SUBROUTINE s_sh94 1623 1624 SUBROUTINE s_sf12 1625 1626 !!---------------------------------------------------------------------- 1627 !! *** ROUTINE s_sf12 *** 1628 !! 1629 !! ** Purpose : stretch the s-coordinate system 1630 !! 1631 !! ** Method : s-coordinate stretch using the Siddorn and Furner 2012? 1632 !! mixed S/sigma/Z coordinate 1633 !! 1634 !! This method allows the maintenance of fixed surface and or 1635 !! bottom cell resolutions (cf. geopotential coordinates) 1636 !! within an analytically derived stretched S-coordinate framework. 1637 !! 1638 !! 1639 !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 1640 !!---------------------------------------------------------------------- 1641 ! 1642 INTEGER :: ji, jj, jk ! dummy loop argument 1643 REAL(wp) :: fsmth ! smoothing around critical depth 1644 REAL(wp) :: zss, zbb ! Surface and bottom cell thickness in sigma space 1645 ! 1646 REAL(wp), POINTER, DIMENSION(:,:,:) :: gsigw3, gsigt3, gsi3w3 1647 REAL(wp), POINTER, DIMENSION(:,:,:) :: esigt3, esigw3, esigtu3, esigtv3, esigtf3, esigwu3, esigwv3 1648 1649 !<