- Timestamp:
- 2020-06-24T14:38:26+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE
- Files:
-
- 3 added
- 15 edited
-
asminc.F90 (modified) (4 diffs)
-
diawri.F90 (modified) (12 diffs)
-
dom_oce.F90 (modified) (5 diffs)
-
domain.F90 (modified) (6 diffs)
-
domvvl.F90 (modified) (11 diffs)
-
domzgr_substitute.h90 (added)
-
dynadv.F90 (modified) (1 diff)
-
dynatf.F90 (modified) (2 diffs)
-
dynkeg.F90 (modified) (2 diffs)
-
dynldf_lap_blp.F90 (modified) (7 diffs)
-
dynvor.F90 (modified) (21 diffs)
-
ldfdyn.F90 (modified) (3 diffs)
-
nemogcm.F90 (modified) (2 diffs)
-
sbcice_cice.F90 (modified) (3 diffs)
-
step.F90 (modified) (15 diffs)
-
stepLF.F90 (added)
-
stpRK3.F90 (added)
-
stpctl.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/asminc.F90
r12614 r13151 95 95 !! * Substitutions 96 96 # include "do_loop_substitute.h90" 97 !!st10 98 # include "domzgr_substitute.h90" 99 !!st10 97 100 !!---------------------------------------------------------------------- 98 101 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 417 420 & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * u_bkginc(ji-1,jj,jk) & 418 421 & + e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * v_bkginc(ji,jj ,jk) & 419 & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk) ) / e3t(ji,jj,jk,Kmm) 422 & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk) ) & 423 & / e3t(ji,jj,jk,Kmm) 420 424 END_2D 421 425 CALL lbc_lnk( 'asminc', zhdiv, 'T', 1. ) ! lateral boundary cond. (no sign change) … … 758 762 ! 759 763 ssh(:,:,Kbb) = ssh(:,:,Kmm) ! Update before fields 764 !!st11 765 #if ! defined key_qco 760 766 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 761 !!gm why not e3u(:,:,:,Kbb), e3v(:,:,:,Kbb), gdept(:,:,:,Kbb) ???? 767 #endif 768 !!st11 769 !!gm why not e3u(:,:,:,Kbb), e3v(:,:,:,Kbb), gdept(:,:,jk,Kbb) ???? 762 770 ! 763 771 DEALLOCATE( ssh_bkg ) … … 970 978 ! ! set to bottom of a level 971 979 ! DO jk = jpk-1, 2, -1 972 ! IF ((mld > gdepw(ji,jj,jk )) .and. (mld < gdepw(ji,jj,jk+1))) THEN973 ! mld=gdepw(ji,jj,jk+1 )980 ! IF ((mld > gdepw(ji,jj,jk,Kmm)) .and. (mld < gdepw(ji,jj,jk+1,Kmm))) THEN 981 ! mld=gdepw(ji,jj,jk+1,Kmm) 974 982 ! jkmax=jk 975 983 ! ENDIF -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/diawri.F90
r12667 r13151 85 85 !! * Substitutions 86 86 # include "do_loop_substitute.h90" 87 !!st12 88 # include "domzgr_substitute.h90" 87 89 !!---------------------------------------------------------------------- 88 90 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 137 139 CALL iom_put("e3v_0", e3v_0(:,:,:) ) 138 140 ! 139 CALL iom_put( "e3t" , e3t(:,:,:,Kmm) ) 140 CALL iom_put( "e3u" , e3u(:,:,:,Kmm) ) 141 CALL iom_put( "e3v" , e3v(:,:,:,Kmm) ) 142 CALL iom_put( "e3w" , e3w(:,:,:,Kmm) ) 143 IF( iom_use("e3tdef") ) & 144 CALL iom_put( "e3tdef" , ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 145 141 !!st13 142 #if ! defined key_qco 143 IF ( iom_use("e3t") .OR. iom_use("e3tdef") ) THEN ! time-varying e3t 144 DO jk = 1, jpk 145 z3d(:,:,jk) = e3t(:,:,jk,Kmm) 146 END DO 147 CALL iom_put( "e3t" , z3d(:,:,:) ) 148 CALL iom_put( "e3tdef" , ( ( z3d(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 ) 149 ENDIF 150 IF ( iom_use("e3u") ) THEN ! time-varying e3u 151 DO jk = 1, jpk 152 z3d(:,:,jk) = e3u(:,:,jk,Kmm) 153 END DO 154 CALL iom_put( "e3u" , z3d(:,:,:) ) 155 ENDIF 156 IF ( iom_use("e3v") ) THEN ! time-varying e3v 157 DO jk = 1, jpk 158 z3d(:,:,jk) = e3v(:,:,jk,Kmm) 159 END DO 160 CALL iom_put( "e3v" , z3d(:,:,:) ) 161 ENDIF 162 IF ( iom_use("e3w") ) THEN ! time-varying e3w 163 DO jk = 1, jpk 164 z3d(:,:,jk) = e3w(:,:,jk,Kmm) 165 END DO 166 CALL iom_put( "e3w" , z3d(:,:,:) ) 167 ENDIF 168 #endif 169 !!st13 146 170 IF( ll_wd ) THEN 147 171 CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) ) ! sea surface height (brought back to the reference used for wetting and drying) … … 351 375 ! 352 376 REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace 353 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace 377 !!st14 378 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d, ze3t, zgdept ! 3D workspace 354 379 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl ! ABL 3D workspace 355 380 !!---------------------------------------------------------------------- … … 391 416 it = kt 392 417 itmod = kt - nit000 + 1 393 418 !!st15 419 ! store e3t for subsitute 420 DO jk = 1, jpk 421 ze3t (:,:,jk) = e3t (:,:,jk,Kmm) 422 zgdept(:,:,jk) = gdept(:,:,jk,Kmm) 423 END DO 424 !!st15 394 425 395 426 ! 1. Define NETCDF files and fields at beginning of first time step … … 514 545 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 515 546 IF( .NOT.ln_linssh ) THEN 516 CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t(:,:,:,Kmm)547 CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! ze3t(:,:,:,Kmm) 517 548 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 518 CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t(:,:,:,Kmm)549 CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! ze3t(:,:,:,Kmm) 519 550 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 520 CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t(:,:,:,Kmm)551 CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! ze3t(:,:,:,Kmm) 521 552 & jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout ) 522 553 ENDIF … … 700 731 WRITE(numout,*) '~~~~~~ ' 701 732 ENDIF 702 733 !!st16 703 734 IF( .NOT.ln_linssh ) THEN 704 CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T ) ! heat content 705 CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T ) ! salt content 706 CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT ) ! sea surface heat content 707 CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT ) ! sea surface salinity content 735 CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T ) ! heat content 736 CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T ) ! salt content 737 CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content 738 CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content 739 !!st16 708 740 ELSE 709 741 CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T ) ! temperature … … 713 745 ENDIF 714 746 IF( .NOT.ln_linssh ) THEN 715 zw3d(:,:,:) = ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 716 CALL histwrite( nid_T, "vovvle3t", it, e3t (:,:,:,Kmm) , ndim_T , ndex_T ) ! level thickness 717 CALL histwrite( nid_T, "vovvldep", it, gdept(:,:,:,Kmm) , ndim_T , ndex_T ) ! t-point depth 747 !!st17 if ! defined key_qco 748 zw3d(:,:,:) = ( ( ze3t(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 749 CALL histwrite( nid_T, "vovvle3t", it, ze3t (:,:,:) , ndim_T , ndex_T ) ! level thickness 750 CALL histwrite( nid_T, "vovvldep", it, zgdept , ndim_T , ndex_T ) ! t-point depth 751 !!st17 718 752 CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation 719 753 ENDIF … … 854 888 !! 855 889 INTEGER :: inum, jk 890 !!st18 TBR 891 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, zgdept ! 3D workspace !!st patch to use substitution 856 892 !!---------------------------------------------------------------------- 857 893 ! … … 860 896 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created ' 861 897 IF(lwp) WRITE(numout,*) ' and named :', cdfile_name, '...nc' 862 898 !!st19 TBR 899 DO jk = 1, jpk 900 ze3t(:,:,jk) = e3t(:,:,jk,Kmm) 901 zgdept(:,:,jk) = gdept(:,:,jk,Kmm) 902 END DO 903 !!st19 863 904 #if defined key_si3 864 905 CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl ) … … 878 919 ENDIF 879 920 CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep ) ! now k-velocity 880 CALL iom_rstput( 0, 0, inum, 'ht' , ht 921 CALL iom_rstput( 0, 0, inum, 'ht' , ht(:,:) ) ! now water column height 881 922 882 923 IF ( ln_isf ) THEN … … 885 926 CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav ) ! now k-velocity 886 927 CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav ) ! now k-velocity 887 CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav, 8)) ! now k-velocity888 CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav, 8)) ! now k-velocity889 CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav, 8), ktype = jp_i1 )928 CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,wp) ) ! now k-velocity 929 CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,wp) ) ! now k-velocity 930 CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,wp), ktype = jp_i1 ) 890 931 END IF 891 932 IF (ln_isfpar_mlt) THEN 892 CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par, 8)) ! now k-velocity933 CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,wp) ) ! now k-velocity 893 934 CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par ) ! now k-velocity 894 935 CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par ) ! now k-velocity 895 936 CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par ) ! now k-velocity 896 CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par, 8)) ! now k-velocity897 CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par, 8)) ! now k-velocity898 CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par, 8), ktype = jp_i1 )937 CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,wp) ) ! now k-velocity 938 CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,wp) ) ! now k-velocity 939 CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,wp), ktype = jp_i1 ) 899 940 END IF 900 941 END IF 901 942 ! 902 943 IF( ALLOCATED(ahtu) ) THEN 903 944 CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point … … 914 955 CALL iom_rstput( 0, 0, inum, 'sozotaux', utau ) ! i-wind stress 915 956 CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress 916 IF( .NOT.ln_linssh ) THEN 917 CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept(:,:,:,Kmm) ) ! T-cell depth 918 CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t(:,:,:,Kmm) ) ! T-cell thickness 957 !!st20 TBR 958 IF( .NOT.ln_linssh ) THEN 959 CALL iom_rstput( 0, 0, inum, 'vovvldep', zgdept ) ! T-cell depth 960 CALL iom_rstput( 0, 0, inum, 'vovvle3t', ze3t ) ! T-cell thickness 919 961 END IF 920 962 IF( ln_wave .AND. ln_sdw ) THEN -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dom_oce.F90
r12614 r13151 71 71 ! ! = 6 cyclic East-West AND North fold F-point pivot 72 72 ! ! = 7 bi-cyclic East-West AND North-South 73 LOGICAL, PUBLIC :: l_Iperio, l_Jperio ! should we explicitely take care I/J periodicity 74 75 ! !domain MPP decomposition parameters73 LOGICAL, PUBLIC :: l_Iperio, l_Jperio ! should we explicitely take care I/J periodicity 74 75 ! !: domain MPP decomposition parameters 76 76 INTEGER , PUBLIC :: nimpp, njmpp !: i- & j-indexes for mpp-subdomain left bottom 77 77 INTEGER , PUBLIC :: nreci, nrecj !: overlap region in i and j … … 136 136 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3vw_0 !: vw-vert. scale factor [m] 137 137 ! ! time-dependent scale factors 138 !!st1 139 #if ! defined key_qco 138 140 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e3t, e3u, e3v, e3w, e3uw, e3vw !: vert. scale factor [m] 139 141 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3f !: F-point vert. scale factor [m] 142 #endif 143 ! ! time-dependent ratio ssh / h_0 144 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: r3t, r3u, r3v !: time-dependent ratio at t-, u- and v-point [-] 145 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r3f !: mid-time-level ratio at f-point [-] 146 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r3t_f, r3u_f, r3v_f !: now time-filtered ratio at t-, u- and v-point [-] 140 147 141 148 ! ! reference depths of cells 142 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_0 !: t- depth [m]143 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdepw_0 !: w- depth [m]144 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gde3w_0 !: w- depth (sum of e3w) [m]149 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_0 !: t- depth [m] 150 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdepw_0 !: w- depth [m] 151 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gde3w_0 !: w- depth (sum of e3w) [m] 145 152 ! ! time-dependent depths of cells 146 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: gdept, gdepw 147 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gde3w 153 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: gdept, gdepw 154 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gde3w 155 !!st2 156 ! ! reference heights of ocean water column and its inverse 157 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0, r1_ht_0 !: t-depth [m] and [1/m] 158 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0, r1_hu_0 !: u-depth [m] and [1/m] 159 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_0, r1_hv_0 !: v-depth [m] and [1/m] 160 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hf_0, r1_hf_0 !: f-depth [m] and [1/m] 161 ! ! time-dependent heights of ocean water column 162 #if ! defined key_qco 163 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht !: t-points [m] 164 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hu, r1_hu !: u-depth [m] and [1/m] 165 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hv, r1_hv !: v-depth [m] and [1/m] 166 #endif 167 !!st2 148 168 149 ! ! reference heights of water column150 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0 !: t-depth [m]151 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0 !: u-depth [m]152 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_0 !: v-depth [m]153 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hf_0 !: f-depth [m]154 ! ! inverse of reference heights of water column155 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ht_0 !: t-depth [1/m]156 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hu_0 !: u-depth [1/m]157 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hv_0 !: v-depth [1/m]158 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hf_0 !: f-depth [1/m]159 160 ! time-dependent heights of water column161 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht !: height of water column at T-points [m]162 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hu, hv, r1_hu, r1_hv !: height of water column [m] and reciprocal [1/m]163 164 169 INTEGER, PUBLIC :: nla10 !: deepest W level Above ~10m (nlb10 - 1) 165 170 INTEGER, PUBLIC :: nlb10 !: shallowest W level Bellow ~10m (nla10 + 1) … … 176 181 !! --------------------------------------------------------------------- 177 182 !!gm Proposition of new name for top/bottom vertical indices 178 ! INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mtk_t, mtk_u, mtk_v !: top first wet T-, U-, V-, F-level (ISF)179 ! INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbk_t, mbk_u, mbk_v !: bottom last wet T-, U-and V-level183 ! INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mtk_t, mtk_u, mtk_v !: top first wet T-, U-, and V-level (ISF) 184 ! INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbk_t, mbk_u, mbk_v !: bottom last wet T-, U-, and V-level 180 185 !!gm 181 186 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mbkt, mbku, mbkv !: bottom last wet T-, U- and V-level … … 244 249 END FUNCTION Agrif_CFixed 245 250 #endif 246 251 !!st3: dom_oce_alloc modified to ease the ifdef if necessary (gm stuff) 247 252 INTEGER FUNCTION dom_oce_alloc() 248 253 !!---------------------------------------------------------------------- 249 INTEGER, DIMENSION(12) :: ierr 254 INTEGER :: ii 255 INTEGER, DIMENSION(30) :: ierr 250 256 !!---------------------------------------------------------------------- 251 i err(:) = 0257 ii = 0 ; ierr(:) = 0 252 258 ! 253 ALLOCATE( mig(jpi), mjg(jpj), STAT=ierr(1) ) 254 ! 255 ALLOCATE( mi0(jpiglo) , mi1 (jpiglo), mj0(jpjglo) , mj1 (jpjglo) , & 256 & tpol(jpiglo) , fpol(jpiglo) , STAT=ierr(2) ) 257 ! 259 ii = ii+1 260 ALLOCATE( mig(jpi), mjg(jpj), STAT=ierr(ii) ) 261 ! 262 ii = ii+1 263 ALLOCATE( mi0 (jpiglo) , mi1 (jpiglo), mj0(jpjglo) , mj1 (jpjglo) , & 264 & tpol(jpiglo) , fpol(jpiglo) , STAT=ierr(ii) ) 265 ! 266 ii = ii+1 258 267 ALLOCATE( glamt(jpi,jpj) , glamu(jpi,jpj) , glamv(jpi,jpj) , glamf(jpi,jpj) , & 259 268 & gphit(jpi,jpj) , gphiu(jpi,jpj) , gphiv(jpi,jpj) , gphif(jpi,jpj) , & … … 266 275 & e1e2v(jpi,jpj) , r1_e1e2v(jpi,jpj) , e1_e2v(jpi,jpj) , & 267 276 & e1e2f(jpi,jpj) , r1_e1e2f(jpi,jpj) , & 268 & ff_f (jpi,jpj) , ff_t (jpi,jpj) , STAT=ierr(3) ) 269 ! 277 & ff_f (jpi,jpj) , ff_t (jpi,jpj) , STAT=ierr(ii) ) 278 ! 279 ii = ii+1 280 ALLOCATE( e3t_0(jpi,jpj,jpk) , e3u_0 (jpi,jpj,jpk) , e3v_0 (jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) , & 281 & e3w_0(jpi,jpj,jpk) , e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , STAT=ierr(ii) ) 282 ! 283 ii = ii+1 270 284 ALLOCATE( gdept_0(jpi,jpj,jpk) , gdepw_0(jpi,jpj,jpk) , gde3w_0(jpi,jpj,jpk) , & 271 & gdept (jpi,jpj,jpk,jpt) , gdepw (jpi,jpj,jpk,jpt) , gde3w (jpi,jpj,jpk) , STAT=ierr(4) ) 272 ! 273 ALLOCATE( e3t_0(jpi,jpj,jpk) , e3u_0(jpi,jpj,jpk) , e3v_0(jpi,jpj,jpk) , e3f_0(jpi,jpj,jpk) , e3w_0(jpi,jpj,jpk) , & 274 & e3t (jpi,jpj,jpk,jpt) , e3u (jpi,jpj,jpk,jpt) , e3v (jpi,jpj,jpk,jpt) , e3f (jpi,jpj,jpk) , e3w (jpi,jpj,jpk,jpt) , & 275 & e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , & 276 & e3uw (jpi,jpj,jpk,jpt) , e3vw (jpi,jpj,jpk,jpt) , STAT=ierr(5) ) 277 ! 278 ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) , hf_0(jpi,jpj) , & 279 & r1_ht_0(jpi,jpj) , r1_hu_0(jpi,jpj) , r1_hv_0(jpi,jpj) , r1_hf_0(jpi,jpj) , & 280 & ht (jpi,jpj) , hu (jpi,jpj,jpt) , hv (jpi,jpj,jpt) , & 281 & r1_hu (jpi,jpj,jpt) , r1_hv (jpi,jpj,jpt) , STAT=ierr(6) ) 282 ! 283 ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(7) ) 284 ! 285 ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(8) ) 286 ! 287 ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) , & 285 & gdept (jpi,jpj,jpk,jpt) , gdepw (jpi,jpj,jpk,jpt) , gde3w (jpi,jpj,jpk) , STAT=ierr(ii) ) 286 ! 287 !!st4 288 #if ! defined key_qco 289 ii = ii+1 290 ALLOCATE( e3t(jpi,jpj,jpk,jpt) , e3u (jpi,jpj,jpk,jpt) , e3v (jpi,jpj,jpk,jpt) , e3f(jpi,jpj,jpk) , & 291 & e3w(jpi,jpj,jpk,jpt) , e3uw(jpi,jpj,jpk,jpt) , e3vw(jpi,jpj,jpk,jpt) , STAT=ierr(ii) ) 292 #endif 293 !!st4 294 ! 295 ii = ii+1 296 ALLOCATE( r3t (jpi,jpj,jpt) , r3u (jpi,jpj,jpt) , r3v (jpi,jpj,jpt) , r3f (jpi,jpj) , & 297 & r3t_f(jpi,jpj) , r3u_f(jpi,jpj) , r3v_f(jpi,jpj) , STAT=ierr(ii) ) 298 ! 299 ii = ii+1 300 ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) , hf_0(jpi,jpj) , & 301 & r1_ht_0(jpi,jpj) , r1_hu_0(jpi,jpj) , r1_hv_0(jpi,jpj), r1_hf_0(jpi,jpj) , STAT=ierr(ii) ) 302 ! 303 #if ! defined key_qco 304 ii = ii+1 305 ALLOCATE( ht (jpi,jpj) , hu (jpi,jpj,jpt), hv (jpi,jpj,jpt) , & 306 & r1_hu (jpi,jpj,jpt), r1_hv (jpi,jpj,jpt) , STAT=ierr(ii) ) 307 #endif 308 ! 309 ii = ii+1 310 ALLOCATE( risfdep(jpi,jpj) , bathy(jpi,jpj) , STAT=ierr(ii) ) 311 ! 312 ii = ii+1 313 ALLOCATE( gdept_1d(jpk) , gdepw_1d(jpk) , e3t_1d(jpk) , e3w_1d(jpk) , STAT=ierr(ii) ) 314 ! 315 ii = ii+1 316 ALLOCATE( tmask_i(jpi,jpj) , tmask_h(jpi,jpj) , & 288 317 & ssmask (jpi,jpj) , ssumask(jpi,jpj) , ssvmask(jpi,jpj) , ssfmask(jpi,jpj) , & 289 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj) , STAT=ierr(9) ) 290 ! 291 ALLOCATE( mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), mikf(jpi,jpj), STAT=ierr(10) ) 292 ! 293 ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) , & 294 & vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , STAT=ierr(11) ) 295 ! 296 ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(12) ) 318 & mbkt (jpi,jpj) , mbku (jpi,jpj) , mbkv (jpi,jpj) , STAT=ierr(ii) ) 319 ! 320 ii = ii+1 321 ALLOCATE( mikt(jpi,jpj), miku(jpi,jpj), mikv(jpi,jpj), mikf(jpi,jpj), STAT=ierr(ii) ) 322 ! 323 ii = ii+1 324 ALLOCATE( tmask(jpi,jpj,jpk) , umask(jpi,jpj,jpk) , & 325 & vmask(jpi,jpj,jpk) , fmask(jpi,jpj,jpk) , STAT=ierr(ii) ) 326 ! 327 ii = ii+1 328 ALLOCATE( wmask(jpi,jpj,jpk) , wumask(jpi,jpj,jpk), wvmask(jpi,jpj,jpk) , STAT=ierr(ii) ) 297 329 ! 298 330 dom_oce_alloc = MAXVAL(ierr) -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domain.F90
r12822 r13151 34 34 USE dommsk ! domain: set the mask system 35 35 USE domwri ! domain: write the meshmask file 36 !!st5 37 #if ! defined key_qco 36 38 USE domvvl ! variable volume 39 #else 40 USE domqco ! variable volume 41 #endif 42 !!st5 37 43 USE c1d ! 1D configuration 38 44 USE dyncor_c1d ! 1D configuration: Coriolis term (cor_c1d routine) … … 78 84 CHARACTER (len=*), INTENT(in) :: cdstr ! model: NEMO or SAS. Determines core restart variables 79 85 ! 80 INTEGER :: ji, jj, jk, ik ! dummy loop indices 86 !!st6 87 INTEGER :: ji, jj, jk, jt ! dummy loop indices 88 !!st6 81 89 INTEGER :: iconf = 0 ! local integers 82 90 CHARACTER (len=64) :: cform = "(A12, 3(A13, I7))" … … 114 122 CASE( 7 ) ; WRITE(numout,*) ' (i.e. cyclic east-west and north-south)' 115 123 CASE DEFAULT 116 CALL ctl_stop( ' jperio is out of range' )124 CALL ctl_stop( 'dom_init: jperio is out of range' ) 117 125 END SELECT 118 126 WRITE(numout,*) ' Ocean model configuration used:' … … 144 152 IF( ln_closea ) CALL dom_clo ! Read in masks to define closed seas and lakes 145 153 146 CALL dom_zgr( ik_top, ik_bot ) ! Vertical mesh and bathymetry 154 CALL dom_zgr( ik_top, ik_bot ) ! Vertical mesh and bathymetry (return top and bottom ocean t-level indices) 147 155 148 156 CALL dom_msk( ik_top, ik_bot ) ! Masks 149 150 ht_0(:,:) = 0._wp 157 ! 158 ht_0(:,:) = 0._wp ! Reference ocean thickness 151 159 hu_0(:,:) = 0._wp 152 160 hv_0(:,:) = 0._wp … … 190 198 ! r1_e1e2t(ji,jj) = r1_e1e2t(ji,jj) / zcoeff 191 199 !!an45 200 !!st7 : make it easier to use key_qco condition (gm stuff) 201 #if defined key_qco 202 ! !== initialisation of time varying coordinate ==! Quasi-Euerian coordinate case 203 ! 204 IF( .NOT.l_offline ) CALL dom_qco_init( Kbb, Kmm, Kaa ) 205 ! 206 IF( ln_linssh ) CALL ctl_stop('STOP','domain: key_qco and ln_linssh = T are incompatible') 207 ! 208 #else 192 209 ! !== time varying part of coordinate system ==! 193 210 ! 194 211 IF( ln_linssh ) THEN != Fix in time : set to the reference one for all 195 !196 ! before ! now ! after !197 gdept(:,:,:, Kbb) = gdept_0 ; gdept(:,:,:,Kmm) = gdept_0 ; gdept(:,:,:,Kaa) = gdept_0 ! depth of grid-points198 gdepw(:,:,:, Kbb) = gdepw_0 ; gdepw(:,:,:,Kmm) = gdepw_0 ; gdepw(:,:,:,Kaa) = gdepw_0 !199 gde3w = gde3w_0 ! --- !200 !201 e3t(:,:,:,Kbb) = e3t_0 ; e3t(:,:,:,Kmm) = e3t_0 ; e3t(:,:,:,Kaa) = e3t_0 ! scale factors202 e3u(:,:,:,Kbb) = e3u_0 ; e3u(:,:,:,Kmm) = e3u_0 ; e3u(:,:,:,Kaa) = e3u_0 !203 e3v(:,:,:,Kbb) = e3v_0 ; e3v(:,:,:,Kmm) = e3v_0 ; e3v(:,:,:,Kaa) = e3v_0 !204 e3f = e3f_0 ! --- !205 e3w(:,:,:,Kbb) = e3w_0 ; e3w(:,:,:,Kmm) = e3w_0 ; e3w(:,:,:,Kaa) = e3w_0 !206 e3uw(:,:,:,Kbb) = e3uw_0 ; e3uw(:,:,:,Kmm) = e3uw_0 ; e3uw(:,:,:,Kaa) = e3uw_0 !207 e3vw(:,:,:,Kbb) = e3vw_0 ; e3vw(:,:,:,Kmm) = e3vw_0 ; e3vw(:,:,:,Kaa) = e3vw_0 !208 !209 z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF210 z1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:))211 ! 212 ! before ! now ! after !213 ht = ht_0 ! ! water column thickness214 hu(:,:,Kbb) = hu_0 ; hu(:,:,Kmm) = hu_0 ; hu(:,:,Kaa) = hu_0 !215 hv(:,:,Kbb) = hv_0 ; hv(:,:,Kmm) = hv_0 ; hv(:,:,Kaa) = hv_0 !216 r1_h u(:,:,Kbb) = z1_hu_0 ; r1_hu(:,:,Kmm) = z1_hu_0 ; r1_hu(:,:,Kaa) = z1_hu_0 ! inverse of water column thickness217 r1_hv(:,:,Kbb) = z1_hv_0 ; r1_hv(:,:,Kmm) = z1_hv_0 ; r1_hv(:,:,Kaa) = z1_hv_0 !218 !212 ! 213 DO jt = 1, jpt ! depth of t- and w-grid-points 214 gdept(:,:,:,jt) = gdept_0(:,:,:) 215 gdepw(:,:,:,jt) = gdepw_0(:,:,:) 216 END DO 217 gde3w(:,:,:) = gde3w_0(:,:,:) ! = gdept as the sum of e3t 218 ! 219 DO jt = 1, jpt ! vertical scale factors 220 e3t(:,:,:,jt) = e3t_0(:,:,:) 221 e3u(:,:,:,jt) = e3u_0(:,:,:) 222 e3v(:,:,:,jt) = e3v_0(:,:,:) 223 e3w(:,:,:,jt) = e3w_0(:,:,:) 224 e3uw(:,:,:,jt) = e3uw_0(:,:,:) 225 e3vw(:,:,:,jt) = e3vw_0(:,:,:) 226 END DO 227 e3f(:,:,:) = e3f_0(:,:,:) 228 ! 229 DO jt = 1, jpt ! water column thickness and its inverse 230 hu(:,:,jt) = hu_0(:,:) 231 hv(:,:,jt) = hv_0(:,:) 232 r1_hu(:,:,jt) = r1_hu_0(:,:) 233 r1_hv(:,:,jt) = r1_hv_0(:,:) 234 END DO 235 ht(:,:) = ht_0(:,:) 219 236 ! 220 237 ELSE != time varying : initialize before/now/after variables 221 238 ! 222 IF( .NOT.l_offline ) CALL dom_vvl_init( Kbb, Kmm, Kaa ) 223 ! 224 ENDIF 239 IF( .NOT.l_offline ) CALL dom_vvl_init( Kbb, Kmm, Kaa ) 240 ! 241 ENDIF 242 #endif 243 !!st7 225 244 ! 226 245 IF( lk_c1d ) CALL cor_c1d ! 1D configuration: Coriolis set at T-point … … 238 257 WRITE(numout,*) 239 258 ENDIF 240 241 259 ! 242 260 END SUBROUTINE dom_init -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domvvl.F90
r13005 r13151 1 1 2 MODULE domvvl 2 3 !!====================================================================== … … 11 12 !!---------------------------------------------------------------------- 12 13 13 !!----------------------------------------------------------------------14 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness15 !! dom_vvl_sf_nxt : Compute next vertical scale factors16 !! dom_vvl_sf_update : Swap vertical scale factors and update the vertical grid17 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another18 !! dom_vvl_rst : read/write restart file19 !! dom_vvl_ctl : Check the vvl options20 !!----------------------------------------------------------------------21 14 USE oce ! ocean dynamics and tracers 22 15 USE phycst ! physical constant … … 36 29 PRIVATE 37 30 38 PUBLIC dom_vvl_init ! called by domain.F9039 PUBLIC dom_vvl_zgr ! called by isfcpl.F9040 PUBLIC dom_vvl_sf_nxt ! called by step.F9041 PUBLIC dom_vvl_sf_update ! called by step.F9042 PUBLIC dom_vvl_interpol ! called by dynnxt.F9043 44 31 ! !!* Namelist nam_vvl 45 32 LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate … … 62 49 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors 63 50 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 64 51 !!stoops 52 #if defined key_qco 53 !!---------------------------------------------------------------------- 54 !! 'key_qco' EMPTY MODULE Quasi-Eulerian vertical coordonate 55 !!---------------------------------------------------------------------- 56 #else 57 !!---------------------------------------------------------------------- 58 !! Default key Old management of time varying vertical coordinate 59 !!---------------------------------------------------------------------- 60 !!st 61 !!---------------------------------------------------------------------- 62 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness 63 !! dom_vvl_sf_nxt : Compute next vertical scale factors 64 !! dom_vvl_sf_update : Swap vertical scale factors and update the vertical grid 65 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 66 !! dom_vvl_rst : read/write restart file 67 !! dom_vvl_ctl : Check the vvl options 68 !!---------------------------------------------------------------------- 69 70 PUBLIC dom_vvl_init ! called by domain.F90 71 PUBLIC dom_vvl_zgr ! called by isfcpl.F90 72 PUBLIC dom_vvl_sf_nxt ! called by step.F90 73 PUBLIC dom_vvl_sf_update ! called by step.F90 74 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 75 PUBLIC dom_vvl_interpol_st! called by dynnxt.F90 76 PUBLIC dom_vvl_sf_nxt_st ! called by step.F90 77 PUBLIC dom_vvl_sf_update_st 78 !!st 79 65 80 !! * Substitutions 66 81 # include "do_loop_substitute.h90" … … 132 147 e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 133 148 ! 134 CALL dom_vvl_zgr (Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column149 CALL dom_vvl_zgr_st(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column 135 150 ! 136 151 END SUBROUTINE dom_vvl_init … … 290 305 END SUBROUTINE dom_vvl_zgr 291 306 307 308 SUBROUTINE dom_vvl_zgr_st(Kbb, Kmm, Kaa) 309 !!---------------------------------------------------------------------- 310 !! *** ROUTINE dom_vvl_init *** 311 !! 312 !! ** Purpose : Interpolation of all scale factors, 313 !! depths and water column heights 314 !! 315 !! ** Method : - interpolate scale factors 316 !! 317 !! ** Action : - e3t_(n/b) and tilde_e3t_(n/b) 318 !! - Regrid: e3(u/v)_n 319 !! e3(u/v)_b 320 !! e3w_n 321 !! e3(u/v)w_b 322 !! e3(u/v)w_n 323 !! gdept_n, gdepw_n and gde3w_n 324 !! - h(t/u/v)_0 325 !! - frq_rst_e3t and frq_rst_hdv 326 !! 327 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 328 !!---------------------------------------------------------------------- 329 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 330 !!---------------------------------------------------------------------- 331 INTEGER :: ji, jj, jk 332 INTEGER :: ii0, ii1, ij0, ij1 333 REAL(wp):: zcoef 334 !!---------------------------------------------------------------------- 335 ! 336 ! !== Set of all other vertical scale factors ==! (now and before) 337 ! ! Horizontal interpolation of e3t 338 CALL dom_vvl_interpol_st( r3u(:,:,Kbb), e3u(:,:,:,Kbb), 'U' ) ! from T to U 339 CALL dom_vvl_interpol_st( r3u(:,:,Kmm), e3u(:,:,:,Kmm), 'U' ) 340 CALL dom_vvl_interpol_st( r3v(:,:,Kbb), e3v(:,:,:,Kbb), 'V' ) ! from T to V 341 CALL dom_vvl_interpol_st( r3v(:,:,Kmm), e3v(:,:,:,Kmm), 'V' ) 342 CALL dom_vvl_interpol_st( r3f(:,:), e3f(:,:,:), 'F' ) ! from U to F 343 ! ! Vertical interpolation of e3t,u,v 344 CALL dom_vvl_interpol_st( r3t(:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) ! from T to W 345 CALL dom_vvl_interpol_st( r3t(:,:,Kbb), e3w (:,:,:,Kbb), 'W' ) 346 CALL dom_vvl_interpol_st( r3u(:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) ! from U to UW 347 CALL dom_vvl_interpol_st( r3u(:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 348 CALL dom_vvl_interpol_st( r3v(:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) ! from V to UW 349 CALL dom_vvl_interpol_st( r3v(:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 350 351 ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...) 352 e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm) 353 e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) 354 e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm) 355 ! 356 DO_3D_11_11( 1, jpk ) 357 gdepw(ji,jj,jk,Kmm) = gdepw_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kmm)) 358 gdept(ji,jj,jk,Kmm) = gdept_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kmm)) 359 gde3w(ji,jj,jk ) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 360 gdepw(ji,jj,jk,Kbb) = gdepw_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kbb)) 361 gdept(ji,jj,jk,Kbb) = gdept_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kbb)) 362 END_3D 363 ! 364 ! !== thickness of the water column !! (ocean portion only) 365 ht(:,:) = ht_0(:,:) + ssh(:,:,Kmm) 366 hu(:,:,Kbb) = (hu_0(:,:)*(1._wp+r3u(:,:,Kbb))) 367 hv(:,:,Kbb) = (hv_0(:,:)*(1._wp+r3v(:,:,Kbb))) 368 hu(:,:,Kbb) = (hu_0(:,:)*(1._wp+r3u(:,:,Kmm))) 369 hv(:,:,Kbb) = (hv_0(:,:)*(1._wp+r3v(:,:,Kmm))) 370 ! !== inverse of water column thickness ==! (u- and v- points) 371 r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF 372 r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) ) 373 r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) ) 374 r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) ) 375 ! 376 IF(lwxios) THEN 377 ! define variables in restart file when writing with XIOS 378 CALL iom_set_rstw_var_active('e3t_b') 379 CALL iom_set_rstw_var_active('e3t_n') 380 ! 381 ENDIF 382 ! 383 END SUBROUTINE dom_vvl_zgr_st 384 292 385 293 386 SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) … … 572 665 573 666 667 668 SUBROUTINE dom_vvl_sf_nxt_st( kt, Kbb, Kmm, Kaa, kcall ) 669 !!---------------------------------------------------------------------- 670 !! *** ROUTINE dom_vvl_sf_nxt *** 671 !! 672 !! ** Purpose : - compute the after scale factors used in tra_zdf, dynnxt, 673 !! tranxt and dynspg routines 674 !! 675 !! ** Method : - z_star case: Repartition of ssh INCREMENT proportionnaly to the level thickness. 676 !! - z_tilde_case: after scale factor increment = 677 !! high frequency part of horizontal divergence 678 !! + retsoring towards the background grid 679 !! + thickness difusion 680 !! Then repartition of ssh INCREMENT proportionnaly 681 !! to the "baroclinic" level thickness. 682 !! 683 !! ** Action : - hdiv_lf : restoring towards full baroclinic divergence in z_tilde case 684 !! - tilde_e3t_a: after increment of vertical scale factor 685 !! in z_tilde case 686 !! - e3(t/u/v)_a 687 !! 688 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 689 !!---------------------------------------------------------------------- 690 INTEGER, INTENT( in ) :: kt ! time step 691 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time step 692 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 693 ! 694 INTEGER :: ji, jj, jk ! dummy loop indices 695 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 696 REAL(wp) :: z_tmin, z_tmax ! local scalars 697 LOGICAL :: ll_do_bclinic ! local logical 698 REAL(wp), DIMENSION(jpi,jpj) :: zht, z_scale, zwu, zwv, zhdiv 699 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t 700 !!---------------------------------------------------------------------- 701 ! 702 IF( ln_linssh ) RETURN ! No calculation in linear free surface 703 ! 704 IF( ln_timing ) CALL timing_start('dom_vvl_sf_nxt') 705 ! 706 IF( kt == nit000 ) THEN 707 IF(lwp) WRITE(numout,*) 708 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors' 709 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' 710 ENDIF 711 712 ll_do_bclinic = .TRUE. 713 IF( PRESENT(kcall) ) THEN 714 IF( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE. 715 ENDIF 716 717 ! ******************************* ! 718 ! After acale factors at t-points ! 719 ! ******************************* ! 720 ! 721 DO jk = 1, jpkm1 722 e3t(:,:,jk,Kaa) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kaa) ) 723 e3u(:,:,jk,Kaa) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kaa) ) 724 e3v(:,:,jk,Kaa) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kaa) ) 725 END DO 726 ! 727 ! *********************************** ! 728 ! After scale factors at u- v- points ! 729 ! *********************************** ! 730 731 !!st CALL dom_vvl_interpol_st( r3u(:,:,Kaa), e3u(:,:,:,Kaa), 'U' ) 732 !!st CALL dom_vvl_interpol_st( r3v(:,:,Kaa), e3v(:,:,:,Kaa), 'V' ) 733 734 ! *********************************** ! 735 ! After depths at u- v points ! 736 ! *********************************** ! 737 738 !!st hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1) 739 !!st hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1) 740 !!st DO jk = 2, jpkm1 741 !!st hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk) 742 !!st hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk) 743 !!st 744 !!st END DO 745 hu(:,:,Kaa) = (hu_0(:,:)*(1._wp+r3u(:,:,Kaa))) 746 hv(:,:,Kaa) = (hv_0(:,:)*(1._wp+r3v(:,:,Kaa))) 747 ! ! Inverse of the local depth 748 !!gm BUG ? don't understand the use of umask_i here ..... 749 r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) ) 750 r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) ) 751 ! 752 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt') 753 ! 754 END SUBROUTINE dom_vvl_sf_nxt_st 755 756 757 574 758 SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa ) 575 759 !!---------------------------------------------------------------------- … … 672 856 ! 673 857 END SUBROUTINE dom_vvl_sf_update 674 858 859 860 SUBROUTINE dom_vvl_sf_update_st( kt, Kbb, Kmm, Kaa ) 861 !!---------------------------------------------------------------------- 862 !! *** ROUTINE dom_vvl_sf_update *** 863 !! 864 !! ** Purpose : for z tilde case: compute time filter and swap of scale factors 865 !! compute all depths and related variables for next time step 866 !! write outputs and restart file 867 !! 868 !! ** Method : - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE) 869 !! - reconstruct scale factor at other grid points (interpolate) 870 !! - recompute depths and water height fields 871 !! 872 !! ** Action : - tilde_e3t_(b/n) ready for next time step 873 !! - Recompute: 874 !! e3(u/v)_b 875 !! e3w(:,:,:,Kmm) 876 !! e3(u/v)w_b 877 !! e3(u/v)w_n 878 !! gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w 879 !! h(u/v) and h(u/v)r 880 !! 881 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 882 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 883 !!---------------------------------------------------------------------- 884 INTEGER, INTENT( in ) :: kt ! time step 885 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices 886 ! 887 INTEGER :: ji, jj, jk ! dummy loop indices 888 REAL(wp) :: zcoef ! local scalar 889 !!---------------------------------------------------------------------- 890 ! 891 IF( ln_linssh ) RETURN ! No calculation in linear free surface 892 ! 893 IF( ln_timing ) CALL timing_start('dom_vvl_sf_update') 894 ! 895 IF( kt == nit000 ) THEN 896 IF(lwp) WRITE(numout,*) 897 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step' 898 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~' 899 ENDIF 900 ! 901 902 ! Compute all missing vertical scale factor and depths 903 ! ==================================================== 904 ! Horizontal scale factor interpolations 905 ! -------------------------------------- 906 ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt 907 ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 908 909 CALL dom_vvl_interpol_st( r3f(:,:), e3f(:,:,:), 'F' ) 910 911 ! Vertical scale factor interpolations 912 CALL dom_vvl_interpol_st( r3t(:,:,Kmm), e3w(:,:,:,Kmm), 'W' ) 913 CALL dom_vvl_interpol_st( r3u(:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' ) 914 CALL dom_vvl_interpol_st( r3v(:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' ) 915 CALL dom_vvl_interpol_st( r3t(:,:,Kbb), e3w(:,:,:,Kbb), 'W' ) 916 CALL dom_vvl_interpol_st( r3u(:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' ) 917 CALL dom_vvl_interpol_st( r3v(:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' ) 918 919 ! t- and w- points depth (set the isf depth as it is in the initial step) 920 DO_3D_11_11( 1, jpk ) 921 gdepw(ji,jj,jk,Kmm) = gdepw_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kmm)) 922 gdept(ji,jj,jk,Kmm) = gdept_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kmm)) 923 gde3w(ji,jj,jk ) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 924 END_3D 925 926 ! Local depth and Inverse of the local depth of the water 927 ! ------------------------------------------------------- 928 ! 929 ht(:,:) = ht_0(:,:) + ssh(:,:,Kmm) 930 931 ! write restart file 932 ! ================== 933 IF( lrst_oce ) CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' ) 934 ! 935 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_update') 936 ! 937 END SUBROUTINE dom_vvl_sf_update_st 938 939 940 941 SUBROUTINE dom_vvl_interpol_st( rc3, pe3, cdp ) 942 !!--------------------------------------------------------------------- 943 !! *** ROUTINE dom_vvl__interpol *** 944 !! 945 !! ** Purpose : interpolate scale factors from one grid point to another 946 !! 947 !! ** Method : e3_out = e3_0 + interpolation(e3_in - e3_0) 948 !! - horizontal interpolation: grid cell surface averaging 949 !! - vertical interpolation: simple averaging 950 !!---------------------------------------------------------------------- 951 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: rc3 ! input e3 NOT used here (ssh is used instead) 952 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe3 ! scale factor e3 to be updated [m] 953 CHARACTER(LEN=*) , INTENT(in ) :: cdp ! grid point of the scale factor ( 'U', 'V', 'W, 'F', 'UW' or 'VW' ) 954 ! 955 INTEGER :: ji, jj, jk ! dummy loop indices 956 REAL(wp), DIMENSION(jpi,jpj) :: zc3 ! 2D workspace 957 !!---------------------------------------------------------------------- 958 ! 959 SELECT CASE ( cdp ) !== type of interpolation ==! 960 ! 961 CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean 962 DO jk = 1, jpkm1 963 pe3(:,:,jk) = e3u_0(:,:,jk) * ( 1.0_wp + rc3(:,:) ) 964 END DO 965 ! 966 CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean 967 DO jk = 1, jpkm1 968 pe3(:,:,jk) = e3v_0(:,:,jk) * ( 1.0_wp + rc3(:,:) ) 969 END DO 970 ! 971 CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean 972 DO jk = 1, jpkm1 ! Horizontal interpolation of e3f from ssh 973 e3f(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + rc3(:,:) ) 974 END DO 975 ! 976 CASE( 'W' ) !* from T- to W-point : vertical simple mean 977 DO jk = 1, jpk 978 pe3(:,:,jk) = e3w_0(:,:,jk) * ( 1.0_wp + rc3(:,:) ) 979 END DO 980 ! 981 CASE( 'UW' ) !* from U- to UW-point 982 DO jk = 1, jpk 983 pe3(:,:,jk) = e3uw_0(:,:,jk) * ( 1.0_wp + rc3(:,:) ) 984 END DO 985 CASE( 'VW' ) !* from U- to UW-point : vertical simple mean 986 DO jk = 1, jpk 987 pe3(:,:,jk) = e3vw_0(:,:,jk) * ( 1.0_wp + rc3(:,:) ) 988 END DO 989 ! 990 END SELECT 991 ! 992 END SUBROUTINE dom_vvl_interpol_st 993 675 994 676 995 SUBROUTINE dom_vvl_interpol( pssh, pe3, cdp ) … … 723 1042 & + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1) ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj) 724 1043 END_2D 725 !!an dans le cas tourné, hf augmente et trend VOR diminue726 ! DO_2D_10_10727 ! zc3(ji,jj) = ( e1e2t(ji ,jj ) * pssh(ji ,jj ) &728 ! & + e1e2t(ji+1,jj ) * pssh(ji+1,jj ) &729 ! & + e1e2t(ji ,jj+1) * pssh(ji ,jj+1) &730 ! & + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1) ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj) &731 ! & / MAX( tmask(ji,jj) + tmask(ji+1,jj) + tmask(ji,jj+1) + tmask(ji+1,jj+1), 1._wp )732 ! END_2D733 734 1044 CALL lbc_lnk( 'domvvl', zc3(:,:), 'F', 1._wp ) 735 1045 ! … … 771 1081 ! 772 1082 END SUBROUTINE dom_vvl_interpol 773 1083 774 1084 775 1085 SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw ) … … 1036 1346 END SUBROUTINE dom_vvl_ctl 1037 1347 1348 #endif 1349 !!stoops 1350 1038 1351 !!====================================================================== 1039 1352 END MODULE domvvl -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90
r13005 r13151 74 74 CASE( np_VEC_c2 ) 75 75 CALL dyn_keg ( kt, nn_dynkeg, Kmm, puu, pvv, Krhs ) ! vector form : horizontal gradient of kinetic energy 76 !!an SWE : w = 077 76 CALL dyn_zad ( kt, Kmm, puu, pvv, Krhs ) ! vector form : vertical advection 78 77 CASE( np_FLX_c2 ) -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynatf.F90
r12614 r13151 58 58 59 59 PUBLIC dyn_atf ! routine called by step.F90 60 !!st22 61 #if defined key_qco 62 !!---------------------------------------------------------------------- 63 !! 'key_qco' EMPTY ROUTINE Quasi-Eulerian vertical coordonate 64 !!---------------------------------------------------------------------- 65 CONTAINS 66 67 SUBROUTINE dyn_atf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v ) 68 INTEGER , INTENT(in ) :: kt ! ocean time-step index 69 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! before and after time level indices 70 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities to be time filtered 71 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: pe3t, pe3u, pe3v ! scale factors to be time filtered 72 73 WRITE(*,*) 'dyn_atf: You should not have seen this print! error?', kt 74 END SUBROUTINE dyn_atf 75 76 #else 60 77 61 78 !! * Substitutions … … 312 329 END SUBROUTINE dyn_atf 313 330 331 #endif 332 !!st22 314 333 !!========================================================================= 315 334 END MODULE dynatf -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90
r13005 r13151 29 29 30 30 PUBLIC dyn_keg ! routine called by step module 31 PUBLIC dyn_kegAD ! routine called by step module32 31 33 32 INTEGER, PARAMETER, PUBLIC :: nkeg_C2 = 0 !: 2nd order centered scheme (standard scheme) … … 156 155 ! 157 156 END SUBROUTINE dyn_keg 158 159 160 SUBROUTINE dyn_kegAD( kt, kscheme, puu, pvv, pu_rhs, pv_rhs )161 !!----------------------------------------------------------------------162 !! *** ROUTINE dyn_kegAD ***163 !!164 !! ** Purpose : Compute the now momentum trend due to the horizontal165 !! gradient of the horizontal kinetic energy and add it to the166 !! general momentum trend.167 !!168 !! ** Method : * kscheme = nkeg_C2 : 2nd order centered scheme that169 !! conserve kinetic energy. Compute the now horizontal kinetic energy170 !! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]171 !! * kscheme = nkeg_HW : Hollingsworth correction following172 !! Arakawa (2001). The now horizontal kinetic energy is given by:173 !! zhke = 1/6 [ mi-1( 2 * un^2 + ((u(j+1)+u(j-1))/2)^2 )174 !! + mj-1( 2 * vn^2 + ((v(i+1)+v(i-1))/2)^2 ) ]175 !!176 !! Take its horizontal gradient and add it to the general momentum177 !! trend.178 !! u(rhs) = u(rhs) - 1/e1u di[ zhke ]179 !! v(rhs) = v(rhs) - 1/e2v dj[ zhke ]180 !!181 !! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend182 !! - send this trends to trd_dyn (l_trddyn=T) for post-processing183 !!184 !! ** References : Arakawa, A., International Geophysics 2001.185 !! Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.186 !!----------------------------------------------------------------------187 !188 INTEGER , INTENT( in ) :: kt ! ocean time-step index189 INTEGER , INTENT( in ) :: kscheme ! =0/1/2 type of KEG scheme190 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: puu, pvv ! ocean velocities at Kmm191 REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) :: pu_rhs, pv_rhs ! RHS192 !193 INTEGER :: ji, jj, jk ! dummy loop indices194 REAL(wp) :: zu, zv ! local scalars195 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhke196 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv197 !!----------------------------------------------------------------------198 !199 IF( ln_timing ) CALL timing_start('dyn_keg')200 !201 IF( kt == nit000 ) THEN202 IF(lwp) WRITE(numout,*)203 IF(lwp) WRITE(numout,*) 'dyn_kegAD : kinetic energy gradient trend, scheme number=', kscheme204 IF(lwp) WRITE(numout,*) '~~~~~~~~~'205 ENDIF206 207 zhke(:,:,jpk) = 0._wp208 157 209 SELECT CASE ( kscheme ) !== Horizontal kinetic energy at T-point ==!210 !!an45 to be ADDED : que cas C2 - "wet points only" il suffit de x2 le terme quadratic a la coast (nn_dynkeg_adv = 2)211 CASE ( nkeg_C2_wpo ) !-- Standard scheme --!212 DO_3D_01_01( 1, jpkm1 )213 zu = ( puu(ji-1,jj ,jk) * puu(ji-1,jj ,jk) &214 & + puu(ji ,jj ,jk) * puu(ji ,jj ,jk) ) * ( 2._wp - umask(ji-1,jj,jk) * umask(ji,jj,jk) )215 zv = ( pvv(ji ,jj-1,jk) * pvv(ji ,jj-1,jk) &216 & + pvv(ji ,jj ,jk) * pvv(ji ,jj ,jk) ) * ( 2._wp - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )217 zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )218 END_3D219 !!an45220 !221 CASE ( nkeg_C2 ) !-- Standard scheme --!222 DO_3D_01_01( 1, jpkm1 )223 zu = puu(ji-1,jj ,jk) * puu(ji-1,jj ,jk) &224 & + puu(ji ,jj ,jk) * puu(ji ,jj ,jk)225 zv = pvv(ji ,jj-1,jk) * pvv(ji ,jj-1,jk) &226 & + pvv(ji ,jj ,jk) * pvv(ji ,jj ,jk)227 zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )228 END_3D229 !!an 00_00 ?230 CASE ( nkeg_HW ) !-- Hollingsworth scheme --!231 DO_3D_00_00( 1, jpkm1 )232 zu = 8._wp * ( puu(ji-1,jj ,jk) * puu(ji-1,jj ,jk) &233 & + puu(ji ,jj ,jk) * puu(ji ,jj ,jk) ) &234 & + ( puu(ji-1,jj-1,jk) + puu(ji-1,jj+1,jk) ) * ( puu(ji-1,jj-1,jk) + puu(ji-1,jj+1,jk) ) &235 & + ( puu(ji ,jj-1,jk) + puu(ji ,jj+1,jk) ) * ( puu(ji ,jj-1,jk) + puu(ji ,jj+1,jk) )236 !237 zv = 8._wp * ( pvv(ji ,jj-1,jk) * pvv(ji ,jj-1,jk) &238 & + pvv(ji ,jj ,jk) * pvv(ji ,jj ,jk) ) &239 & + ( pvv(ji-1,jj-1,jk) + pvv(ji+1,jj-1,jk) ) * ( pvv(ji-1,jj-1,jk) + pvv(ji+1,jj-1,jk) ) &240 & + ( pvv(ji-1,jj ,jk) + pvv(ji+1,jj ,jk) ) * ( pvv(ji-1,jj ,jk) + pvv(ji+1,jj ,jk) )241 zhke(ji,jj,jk) = r1_48 * ( zv + zu )242 END_3D243 CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )244 !245 END SELECT246 !247 IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** NO alternating direction ***!248 !249 DO_3D_00_00( 1, jpkm1 )250 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)251 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)252 END_3D253 !254 ELSEIF( PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : i-component ***!255 !256 DO_3D_00_00( 1, jpkm1 )257 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)258 END_3D259 !260 ELSEIF( .NOT. PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : j-component ***!261 !262 DO_3D_00_00( 1, jpkm1 )263 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)264 END_3D265 !266 ENDIF267 IF( ln_timing ) CALL timing_stop('dyn_kegAD')268 !269 END SUBROUTINE dyn_kegAD270 158 !!====================================================================== 271 159 END MODULE dynkeg -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynldf_lap_blp.F90
r13005 r13151 20 20 USE in_out_manager ! I/O manager 21 21 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 22 USE lib_mpp ! MPP library 23 ! 24 USE usrdef_nam , ONLY : nn_dynldf_lap_typ ! use laplacian parameter 25 ! 22 26 23 IMPLICIT NONE 27 24 PRIVATE … … 34 31 INTEGER, PUBLIC, PARAMETER :: np_dynldf_lap_symN = 3 ! symmetric laplacian (cartesian) 35 32 36 !INTEGER, PUBLIC, PARAMETER :: nn_dynldf_lap_typ = 1 ! choose type of laplacian (ideally from namelist)33 INTEGER, PUBLIC, PARAMETER :: ln_dynldf_lap_typ = 1 ! choose type of laplacian (ideally from namelist) 37 34 !!anSYM 38 35 !! * Substitutions 39 36 # include "do_loop_substitute.h90" 37 !!st21 38 # include "domzgr_substitute.h90" 40 39 !!---------------------------------------------------------------------- 41 40 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 80 79 WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass 81 80 WRITE(numout,*) '~~~~~~~ ' 82 WRITE(numout,*) ' nn_dynldf_lap_typ = ', nn_dynldf_lap_typ83 SELECT CASE( nn_dynldf_lap_typ ) ! print the choice of operator81 WRITE(numout,*) ' ln_dynldf_lap_typ = ', ln_dynldf_lap_typ 82 SELECT CASE( ln_dynldf_lap_typ ) ! print the choice of operator 84 83 CASE( np_dynldf_lap_rot ) ; WRITE(numout,*) ' ==>>> div-rot laplacian' 85 84 CASE( np_dynldf_lap_sym ) ; WRITE(numout,*) ' ==>>> symmetric laplacian (covariant form)' 86 CASE( np_dynldf_lap_symN) ; WRITE(numout,*) ' ==>>> symmetric laplacian ( cartesianform)'85 CASE( np_dynldf_lap_symN) ; WRITE(numout,*) ' ==>>> symmetric laplacian (simple form)' 87 86 END SELECT 88 87 ENDIF … … 92 91 ENDIF 93 92 ! 94 SELECT CASE( nn_dynldf_lap_typ )93 SELECT CASE( ln_dynldf_lap_typ ) 95 94 ! 96 95 CASE ( np_dynldf_lap_rot ) !== Vorticity-Divergence form ==! … … 102 101 !!gm open question here : e3f at before or now ? probably now... 103 102 !!gm note that ahmf has already been multiplied by fmask 104 zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & 105 & * ( e2v(ji ,jj-1) * pv(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk) & 106 & - e1u(ji-1,jj ) * pu(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk) ) 107 ! ! ahm * div (computed from 2 to jpi/jpj) 103 zcur(ji-1,jj-1) = & 104 & ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & 105 & * ( e2v(ji ,jj-1) * pv(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk) & 106 & - e1u(ji-1,jj ) * pu(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk) ) 107 ! ! ahm * div (computed from 2 to jpi/jpj) 108 108 !!gm note that ahmt has already been multiplied by tmask 109 109 zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb) & … … 160 160 END DO ! End of slab 161 161 ! 162 CASE ( np_dynldf_lap_symN ) !== Symmetric form ==! ( cartesianway)162 CASE ( np_dynldf_lap_symN ) !== Symmetric form ==! (naive way) 163 163 ! 164 164 DO jk = 1, jpkm1 ! Horizontal slab … … 193 193 ! 194 194 CASE DEFAULT ! error 195 CALL ctl_stop('STOP','dyn_ldf_lap: wrong value for nn_dynldf_lap_typ' )195 CALL ctl_stop('STOP','dyn_ldf_lap: wrong value for ln_dynldf_lap_typ' ) 196 196 END SELECT 197 197 ! -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90
r13005 r13151 22 22 !! - ! 2018-04 (G. Madec) add pre-computed gradient for metric term calculation 23 23 !! 4.x ! 2020-03 (G. Madec, A. Nasser) make ln_dynvor_msk truly efficient on relative vorticity 24 !! 4.x ! 2020-03 (G. Madec, A. Nasser) alternate direction computation of vorticity tendancy25 !! ! for ENS, ENE26 24 !!---------------------------------------------------------------------- 27 25 … … 47 45 USE lib_mpp ! MPP library 48 46 USE timing ! Timing 49 !!anAD only50 USE dynkeg, ONLY : dyn_kegAD51 USE dynadv, ONLY : nn_dynkeg52 47 53 48 IMPLICIT NONE … … 58 53 59 54 ! !!* Namelist namdyn_vor: vorticity term 60 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS) 61 LOGICAL, PUBLIC :: ln_dynvor_ens_adVO = .FALSE. !: AD enstrophy conserving scheme (ENS_ad) 62 LOGICAL, PUBLIC :: ln_dynvor_ens_adKE = .FALSE. !: AD enstrophy conserving scheme (ENS_ad) 63 LOGICAL, PUBLIC :: ln_dynvor_ens_adKEVO = .FALSE. !: AD enstrophy conserving scheme (ENS_ad) 64 LOGICAL, PUBLIC :: ln_dynvor_ene !: f-point energy conserving scheme (ENE) 65 LOGICAL, PUBLIC :: ln_dynvor_ene_adVO = .FALSE. !: f-point AD energy conserving scheme (ENE_ad) 66 LOGICAL, PUBLIC :: ln_dynvor_ene_adKE = .FALSE. !: f-point AD energy conserving scheme (ENE_ad) 67 LOGICAL, PUBLIC :: ln_dynvor_ene_adKEVO = .FALSE. !: f-point AD energy conserving scheme (ENE_ad) 68 LOGICAL, PUBLIC :: ln_dynvor_enT !: t-point energy conserving scheme (ENT) 69 LOGICAL, PUBLIC :: ln_dynvor_eeT !: t-point energy conserving scheme (EET) 70 LOGICAL, PUBLIC :: ln_dynvor_een !: energy & enstrophy conserving scheme (EEN) 71 INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 72 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX) 73 LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 55 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS) 56 LOGICAL, PUBLIC :: ln_dynvor_ene !: f-point energy conserving scheme (ENE) 57 LOGICAL, PUBLIC :: ln_dynvor_enT !: t-point energy conserving scheme (ENT) 58 LOGICAL, PUBLIC :: ln_dynvor_eeT !: t-point energy conserving scheme (EET) 59 LOGICAL, PUBLIC :: ln_dynvor_een !: energy & enstrophy conserving scheme (EEN) 60 INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 61 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX) 62 LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 74 63 75 64 INTEGER, PUBLIC :: nvor_scheme !: choice of the type of advection scheme … … 81 70 INTEGER, PUBLIC, PARAMETER :: np_EEN = 4 ! EEN scheme 82 71 INTEGER, PUBLIC, PARAMETER :: np_MIX = 5 ! MIX scheme 83 !!an 84 INTEGER, PUBLIC, PARAMETER :: np_ENS_adKE = 11 ! ENS scheme - AD scheme (KE only) 85 INTEGER, PUBLIC, PARAMETER :: np_ENS_adVO = 12 ! ENS scheme - AD scheme (VOR only) 86 INTEGER, PUBLIC, PARAMETER :: np_ENS_adKEVO = 13 ! ENS scheme - AD scheme (KE+VOR) 87 INTEGER, PUBLIC, PARAMETER :: np_ENE_adKE = 21 ! ENE scheme - AD scheme (KE only) 88 INTEGER, PUBLIC, PARAMETER :: np_ENE_adVO = 22 ! ENE scheme - AD scheme (VOR only) 89 INTEGER, PUBLIC, PARAMETER :: np_ENE_adKEVO = 23 ! ENE scheme - AD scheme (KE+VOR) 90 !!an 91 92 !!an ds step on pourra spécifier la valeur de ntot = np_COR ou np_COR + np_RVO 93 INTEGER, PUBLIC :: ncor, nrvm, ntot ! choice of calculated vorticity 72 73 INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity 94 74 ! ! associated indices: 95 75 INTEGER, PUBLIC, PARAMETER :: np_COR = 1 ! Coriolis (planetary) … … 110 90 !! * Substitutions 111 91 # include "do_loop_substitute.h90" 92 !!st23 93 # include "domzgr_substitute.h90" 112 94 !!---------------------------------------------------------------------- 113 95 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 117 99 CONTAINS 118 100 119 SUBROUTINE dyn_vor( kt, K bb, Kmm, puu, pvv, Krhs)101 SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs ) 120 102 !!---------------------------------------------------------------------- 121 103 !! … … 127 109 !! for futher diagnostics (l_trddyn=T) 128 110 !!---------------------------------------------------------------------- 129 INTEGER :: ji, jj, jk ! dummy loop indice 130 INTEGER , INTENT( in ) :: kt ! ocean time-step index 131 INTEGER , INTENT( in ) :: Kmm, Krhs, Kbb ! ocean time level indices 132 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation 111 INTEGER , INTENT( in ) :: kt ! ocean time-step index 112 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 113 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation 133 114 ! 134 115 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 135 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuu, zvv136 116 !!---------------------------------------------------------------------- 137 117 ! … … 187 167 IF( ln_stcor ) CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 188 168 CASE( np_ENE ) !* energy conserving scheme 189 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 190 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 169 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 191 170 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 192 193 CASE( np_ENE_adVO ) !* energy conserving scheme with alternating direction194 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components195 196 !== Alternative direction - VOR only ==!197 198 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add uu-vorticity trend199 200 ALLOCATE( zuu(jpi,jpj,jpk) )201 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend202 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)203 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )204 CALL vor_ene( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend205 DEALLOCATE( zuu )206 ELSE207 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add uu-vorticity trend208 209 ALLOCATE( zvv(jpi,jpj,jpk) )210 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend211 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)212 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )213 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend214 DEALLOCATE( zvv )215 ENDIF216 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend217 CASE( np_ENE_adKE ) !* energy conserving scheme with alternating direction218 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components219 220 !== Alternative direction - KEG only ==!221 222 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend223 224 ALLOCATE( zuu(jpi,jpj,jpk) )225 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend226 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)227 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )228 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend229 DEALLOCATE( zuu )230 ELSE231 232 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend233 234 ALLOCATE( zvv(jpi,jpj,jpk) )235 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend236 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)237 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )238 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend239 DEALLOCATE( zvv )240 ENDIF241 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend242 243 CASE( np_ENE_adKEVO ) !* energy conserving scheme with alternating direction244 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components245 246 !== Alternative direction - KE + VOR ==!247 248 ALLOCATE( zuu(jpi,jpj,jpk) )249 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend250 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) !251 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)252 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )253 CALL vor_ene( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend254 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )255 DEALLOCATE( zuu )256 ELSE257 258 ALLOCATE( zvv(jpi,jpj,jpk) )259 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend260 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )261 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)262 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )263 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend264 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )265 DEALLOCATE( zvv )266 ENDIF267 171 CASE( np_ENS ) !* enstrophy conserving scheme 268 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )269 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend270 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend271 CASE( np_ENS_adVO ) !* enstrophy conserving scheme with alternating direction272 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components273 274 !== Alternative direction - VOR only ==!275 276 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )277 278 ALLOCATE( zuu(jpi,jpj,jpk) )279 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend280 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)281 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )282 CALL vor_ens( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend283 DEALLOCATE( zuu )284 ELSE285 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) )286 287 ALLOCATE( zvv(jpi,jpj,jpk) )288 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend289 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)290 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )291 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , zvv(:,:,:), pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend292 DEALLOCATE( zvv )293 ENDIF294 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend295 CASE( np_ENS_adKE ) !* enstrophy conserving scheme with alternating direction296 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components297 298 !== Alternative direction - KEG only ==!299 300 172 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 301 302 ALLOCATE( zuu(jpi,jpj,jpk) )303 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend304 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)305 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )306 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend307 DEALLOCATE( zuu )308 ELSE309 310 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend311 312 ALLOCATE( zvv(jpi,jpj,jpk) )313 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend314 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)315 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )316 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend317 DEALLOCATE( zvv )318 ENDIF319 CASE( np_ENS_adKEVO ) !* enstrophy conserving scheme with alternating direction320 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components321 322 !== Alternative direction - KE + VOR ==!323 324 ALLOCATE( zuu(jpi,jpj,jpk) )325 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend326 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) !327 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:)328 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp )329 CALL vor_ens( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend330 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )331 DEALLOCATE( zuu )332 ELSE333 334 ALLOCATE( zvv(jpi,jpj,jpk) )335 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend336 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) )337 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:)338 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp )339 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend340 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) )341 DEALLOCATE( zvv )342 ENDIF343 173 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 344 174 CASE( np_MIX ) !* mixed ene-ens scheme … … 427 257 DO_2D_01_01 428 258 zwt(ji,jj) = r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) & 429 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 259 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) & 260 & * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 430 261 END_2D 431 262 CASE ( np_MET ) !* metric term 432 263 DO_2D_01_01 433 zwt(ji,jj) = ( ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 434 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) * e3t(ji,jj,jk,Kmm) 264 zwt(ji,jj) = ( ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 265 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) & 266 & * e3t(ji,jj,jk,Kmm) 435 267 END_2D 436 268 CASE ( np_CRV ) !* Coriolis + relative vorticity 437 269 DO_2D_01_01 438 zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) & 439 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 270 zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) & 271 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) ) & 272 & * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 440 273 END_2D 441 274 CASE ( np_CME ) !* Coriolis + metric 442 275 DO_2D_01_01 443 zwt(ji,jj) = ( ff_t(ji,jj) * e1e2t(ji,jj) & 444 & + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 445 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) * e3t(ji,jj,jk,Kmm) 276 zwt(ji,jj) = ( ff_t(ji,jj) * e1e2t(ji,jj) & 277 & + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 278 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) & 279 & * e3t(ji,jj,jk,Kmm) 446 280 END_2D 447 281 CASE DEFAULT ! error … … 485 319 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 486 320 !!---------------------------------------------------------------------- 487 INTEGER 488 INTEGER 489 INTEGER 490 REAL(wp), DIMENSION(jpi,jpj,jpk) 491 REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL,INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend321 INTEGER , INTENT(in ) :: kt ! ocean time-step index 322 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 323 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 324 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities 325 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend 492 326 ! 493 327 INTEGER :: ji, jj, jk ! dummy loop indices … … 543 377 END SELECT 544 378 ! 545 IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** NO alternating direction ***! 546 ! 547 ! !== horizontal fluxes and potential vorticity ==! 548 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 549 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 550 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 551 ! 552 ! !== compute and add the vorticity term trend =! 553 DO_2D_00_00 554 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 555 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 556 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 557 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 558 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 559 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 560 END_2D 561 ! 562 ! 563 ELSEIF( PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : i-component ***! 564 ! 565 ! 566 ! !== horizontal fluxes and potential vorticity ==! 567 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 568 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 569 ! 570 ! !== compute and add the vorticity term trend =! 571 DO_2D_00_00 572 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 573 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 574 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 575 END_2D 576 ! 577 ELSEIF( .NOT. PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : j-component ***! 578 ! 579 ! 580 ! !== horizontal fluxes and potential vorticity ==! 581 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 582 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 583 ! 584 ! !== compute and add the vorticity term trend =! 585 DO_2D_00_00 586 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 587 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 588 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 589 END_2D 590 ! 591 ENDIF 379 ! !== horizontal fluxes and potential vorticity ==! 380 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 381 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 382 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 383 ! 384 ! !== compute and add the vorticity term trend =! 385 DO_2D_00_00 386 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 387 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 388 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 389 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 390 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 391 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 392 END_2D 592 393 ! ! =============== 593 394 END DO ! End of slab … … 616 417 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 617 418 !!---------------------------------------------------------------------- 618 INTEGER 619 INTEGER 620 INTEGER 621 REAL(wp), DIMENSION(jpi,jpj,jpk) 622 REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL,INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend419 INTEGER , INTENT(in ) :: kt ! ocean time-step index 420 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 421 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 422 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities 423 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend 623 424 ! 624 425 INTEGER :: ji, jj, jk ! dummy loop indices … … 674 475 ! 675 476 ! 676 !!an wut ? v et u 677 IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** NO alternating direction ***! 678 ! 679 ! !== horizontal fluxes and potential vorticity ==! 680 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 681 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 682 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 683 ! 684 ! !== compute and add the vorticity term trend =! 685 DO_2D_00_00 686 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 687 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 688 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 689 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 690 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 691 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 692 END_2D 693 ! 694 ELSEIF( PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : i-component ***! 695 ! 696 ! 697 ! !== horizontal fluxes and potential vorticity ==! 698 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 699 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 700 ! 701 ! !== compute and add the vorticity term trend =! 702 DO_2D_00_00 703 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 704 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 705 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 706 END_2D 707 ! 708 ELSEIF( .NOT. PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : j-component ***! 709 ! 710 ! 711 ! !== horizontal fluxes and potential vorticity ==! 712 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 713 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 714 ! 715 ! !== compute and add the vorticity term trend =! 716 DO_2D_00_00 717 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 718 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 719 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 720 END_2D 721 ! 722 ENDIF 477 ! !== horizontal fluxes and potential vorticity ==! 478 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 479 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 480 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 481 ! 482 ! !== compute and add the vorticity term trend =! 483 DO_2D_00_00 484 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 485 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 486 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 487 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 488 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 489 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 490 END_2D 723 491 ! ! =============== 724 492 END DO ! End of slab … … 772 540 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 773 541 DO_2D_10_10 774 ze3f = ( e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 775 & + e3t(ji,jj ,jk,Kmm)*tmask(ji,jj ,jk) + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 542 ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & 543 & + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 544 & + e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) & 545 & + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 776 546 IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4._wp / ze3f 777 547 ELSE ; z1_e3f(ji,jj) = 0._wp … … 780 550 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 781 551 DO_2D_10_10 782 ze3f = ( e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 783 & + e3t(ji,jj ,jk,Kmm)*tmask(ji,jj ,jk) + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 552 ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & 553 & + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 554 & + e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) & 555 & + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 784 556 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 785 557 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) … … 1000 772 !! 1001 773 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT, & 1002 & ln_dynvor_een, nn_een_e3f , ln_dynvor_mix, ln_dynvor_msk, & 1003 & ln_dynvor_ens_adVO, ln_dynvor_ens_adKE, ln_dynvor_ens_adKEVO, & ! Alternative direction parameters 1004 & ln_dynvor_ene_adVO, ln_dynvor_ene_adKE, ln_dynvor_ene_adKEVO 774 & ln_dynvor_een, nn_een_e3f , ln_dynvor_mix, ln_dynvor_msk 1005 775 !!---------------------------------------------------------------------- 1006 776 ! … … 1019 789 IF(lwp) THEN ! Namelist print 1020 790 WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme' 1021 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens 1022 WRITE(numout,*) ' f-point energy conserving scheme ln_dynvor_ene 1023 WRITE(numout,*) ' t-point energy conserving scheme ln_dynvor_enT 1024 WRITE(numout,*) ' energy conserving scheme (een using e3t) ln_dynvor_eeT 1025 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een 1026 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f 1027 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix 1028 WRITE(numout,*) ' masked (=T) or unmasked(=F) vorticity ln_dynvor_msk 791 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens 792 WRITE(numout,*) ' f-point energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 793 WRITE(numout,*) ' t-point energy conserving scheme ln_dynvor_enT = ', ln_dynvor_enT 794 WRITE(numout,*) ' energy conserving scheme (een using e3t) ln_dynvor_eeT = ', ln_dynvor_eeT 795 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 796 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f 797 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 798 WRITE(numout,*) ' masked (=T) or unmasked(=F) vorticity ln_dynvor_msk = ', ln_dynvor_msk 1029 799 ENDIF 1030 800 … … 1039 809 DO_3D_10_10( 1, jpk ) 1040 810 IF( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 1041 & + tmask(ji,jj ,jk) + tmask(ji+1,jj 811 & + tmask(ji,jj ,jk) + tmask(ji+1,jj+1,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp 1042 812 END_3D 1043 813 ! … … 1049 819 ioptio = 0 ! type of scheme for vorticity (set nvor_scheme) 1050 820 IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF 1051 IF( ln_dynvor_ens_adVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adVO ; ENDIF1052 IF( ln_dynvor_ens_adKE ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adKE ; ENDIF1053 IF( ln_dynvor_ens_adKEVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adKEVO ; ENDIF1054 821 IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF 1055 IF( ln_dynvor_ene_adVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adVO ; ENDIF1056 IF( ln_dynvor_ene_adKE ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adKE ; ENDIF1057 IF( ln_dynvor_ene_adKEVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adKEVO ; ENDIF1058 822 IF( ln_dynvor_enT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENT ; ENDIF 1059 823 IF( ln_dynvor_eeT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EET ; ENDIF … … 1072 836 CASE( np_VEC_c2 ) 1073 837 IF(lwp) WRITE(numout,*) ' ==>>> vector form dynamics : total vorticity = Coriolis + relative vorticity' 1074 nrvm = np_RVO ! relative vorticity 838 nrvm = np_RVO ! relative vorticity 1075 839 ntot = np_CRV ! relative + planetary vorticity 1076 840 CASE( np_FLX_c2 , np_FLX_ubs ) … … 1102 866 WRITE(numout,*) 1103 867 SELECT CASE( nvor_scheme ) 1104 1105 CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme (ENS)' 1106 CASE( np_ENS_adVO ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adVO) on vorticity only' 1107 CASE( np_ENS_adKE ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adKE) on kinetic energy only' 1108 CASE( np_ENS_adKEVO ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adKEVO) on kinetic energy and vorticity' 1109 1110 CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at F-points) (ENE)' 1111 CASE( np_ENE_adVO ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adVO) on vorticity only' 1112 CASE( np_ENE_adKE ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adKE) on kinetic energy only' 1113 CASE( np_ENE_adKEVO ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adKEVO) on kinetic energy and vorticity' 1114 1115 CASE( np_ENT ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at T-points) (ENT)' 1116 CASE( np_EET ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3t) (EET)' 1117 CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)' 1118 CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)' 868 CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme (ENS)' 869 CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at F-points) (ENE)' 870 CASE( np_ENT ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at T-points) (ENT)' 871 CASE( np_EET ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3t) (EET)' 872 CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)' 873 CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)' 1119 874 END SELECT 1120 875 ENDIF -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/ldfdyn.F90
r13005 r13151 25 25 USE lib_mpp ! distribued memory computing library 26 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 27 ! 28 USE usrdef_nam , ONLY : ln_dynldf_lap_PM 29 ! 27 30 28 IMPLICIT NONE 31 29 PRIVATE … … 62 60 INTEGER , PUBLIC :: nldf_dyn !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 63 61 LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef. 64 !!an 65 !LOGICAL , PUBLIC :: ll_dynldf_lap_PM !: flag for P.Marchand modification on viscosity 66 !!an 62 67 63 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahmt, ahmf !: eddy viscosity coef. at T- and F-points [m2/s or m4/s] 68 64 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dtensq !: horizontal tension squared (Smagorinsky only) … … 327 323 IF( .NOT.l_ldfdyn_time ) THEN !* No time variation 328 324 IF( ln_dynldf_lap ) THEN ! laplacian operator (mask only) 329 !!an ! 330 WRITE(numout,*) ' ln_dynldf_lap_PM = ',ln_dynldf_lap_PM 331 IF( ln_dynldf_lap_PM ) THEN ! laplacian operator (mask only) 325 ahmt(:,:,1:jpkm1) = ahmt(:,:,1:jpkm1) * tmask(:,:,1:jpkm1) 326 WRITE(numout,*) ' ahmt tmask ' 332 327 !! mask tension at the coast (equivalent of ghostpoints at T) 333 DO jk = 1, jpk 334 DO jj = 1, jpjm1 335 DO ji = 1, jpim1 ! NO vector opt. 336 ! si sum(fmask)==3 = mouillé (on touche pas) 337 ! si sum = 2 alors on met a 0 zsum = fmask + fmask + fmask,.. et si zsum < 2 -> 0 sinon = 1 338 zsum = fmask(ji,jj ,jk) + fmask(ji+1,jj ,jk) & 339 & + fmask(ji,jj+1,jk) + fmask(ji+1,jj+1,jk) 340 IF ( zsum < 2._wp ) ahmt(ji,jj,jk) = 0 341 ! 342 !ahmt(ji,jj,jk) = ahmt(ji,jj,jk) * fmask(ji,jj ,jk) * fmask(ji+1,jj ,jk) & 343 ! & * fmask(ji,jj+1,jk) * fmask(ji+1,jj+1,jk) 344 END DO 345 END DO 346 END DO 347 ahmt(jpi,:,1:jpkm1) = 0._wp 348 ahmt(:,jpj,1:jpkm1) = 0._wp 349 WRITE(numout,*) ' ahmt x0' 350 !! apply no slip at the coast (ssfmask = 1 within the domain and at the coast contrary to fmask in free slip) 351 DO jk = 1, jpkm1 352 ahmf(:,:,jk) = ahmf(:,:,jk) * ( 2._wp * ssfmask(:,:) - fmask(:,:,jk) ) 353 END DO 354 WRITE(numout,*) ' ahmf x2' 355 ELSE 356 ! classic boundary condition on the viscosity coefficient 357 ahmt(:,:,1:jpkm1) = ahmt(:,:,1:jpkm1) * tmask(:,:,1:jpkm1) 358 WRITE(numout,*) ' ahmt tmasked ' 359 ahmf(:,:,1:jpkm1) = ahmf(:,:,1:jpkm1) * fmask(:,:,1:jpkm1) 360 WRITE(numout,*) ' ahmf fmasked ' 361 ENDIF 362 !!an ! 328 ! DO jk = 1, jpk 329 ! DO jj = 1, jpjm1 330 ! DO ji = 1, jpim1 ! NO vector opt. 331 ! ! si sum(fmask)==3 = mouillé (on touche pas) 332 ! ! si sum = 2 alors on met a 0 zsum = fmask + fmask + fmask,.. et si zsum < 2 -> 0 sinon = 1 333 ! zsum = fmask(ji,jj ,jk) + fmask(ji+1,jj ,jk) & 334 ! & + fmask(ji,jj+1,jk) + fmask(ji+1,jj+1,jk) 335 ! IF ( zsum < 2._wp ) ahmt(ji,jj,jk) = 0 336 ! ! 337 ! !ahmt(ji,jj,jk) = ahmt(ji,jj,jk) * fmask(ji,jj ,jk) * fmask(ji+1,jj ,jk) & 338 ! ! & * fmask(ji,jj+1,jk) * fmask(ji+1,jj+1,jk) 339 ! END DO 340 ! END DO 341 ! END DO 342 ! ahmt(jpi,:,1:jpkm1) = 0._wp 343 ! ahmt(:,jpj,1:jpkm1) = 0._wp 344 ! WRITE(numout,*) ' an45 ahmt x0' 345 346 ahmf(:,:,1:jpkm1) = ahmf(:,:,1:jpkm1) * fmask(:,:,1:jpkm1) 347 WRITE(numout,*) ' ahmf fmask ' 348 !!an apply no slip at the coast (ssfmask = 1 within the domain and at the coast contrary to fmask in free slip) 349 ! DO jk = 1, jpkm1 350 ! ahmf(:,:,jk) = ahmf(:,:,jk) * ( 2._wp * ssfmask(:,:) - fmask(:,:,jk) ) 351 ! END DO 352 ! WRITE(numout,*) ' an45 ahmf x2' 353 363 354 ELSEIF( ln_dynldf_blp ) THEN ! bilaplacian operator (square root + mask) 364 355 ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1) -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/nemogcm.F90
r12614 r13151 60 60 USE diacfl ! CFL diagnostics (dia_cfl_init routine) 61 61 USE diamlr ! IOM context management for multiple-linear-regression analysis 62 #if defined key_RK3 63 USE stpRK3 64 #elif defined key_qco 65 USE stpLF 66 #else 62 67 USE step ! NEMO time-stepping (stp routine) 68 #endif 63 69 USE isfstp ! ice shelf (isf_stp_init routine) 64 70 USE icbini ! handle bergs, initialisation … … 175 181 IF ( istp == nitend ) elapsed_time = zstptiming - elapsed_time 176 182 ENDIF 177 178 CALL stp ( istp ) 183 #if defined key_RK3 184 CALL stp_RK3 ( istp ) 185 #elif defined key_qco 186 CALL stp_LF ( istp ) 187 #else 188 CALL stp ( istp ) 189 #endif 179 190 istp = istp + 1 180 191 -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/sbcice_cice.F90
r12614 r13151 12 12 USE oce ! ocean dynamics and tracers 13 13 USE dom_oce ! ocean space and time domain 14 !!st8 15 # if ! defined key_qco 14 16 USE domvvl 17 # else 18 USE domqco 19 # endif 20 !!st8 15 21 USE phycst, only : rcp, rho0, r1_rho0, rhos, rhoi 16 22 USE in_out_manager ! I/O manager … … 233 239 !!gm This should be put elsewhere.... (same remark for limsbc) 234 240 !!gm especially here it is assumed zstar coordinate, but it can be ztilde.... 241 !!st9 242 #if defined key_qco 243 IF( .NOT.ln_linssh ) CALL dom_qco_zgr( Kbb, Kmm, Kaa ) ! interpolation scale factor, depth and water column 244 #else 235 245 IF( .NOT.ln_linssh ) THEN 236 246 ! 237 247 DO jk = 1,jpkm1 ! adjust initial vertical scale factors 238 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)* tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )239 e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)* tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )240 END DO248 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)*r1_ht_0(:,:)*tmask(:,:,jk) ) 249 e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)*r1_ht_0(:,:)*tmask(:,:,jk) ) 250 END DO 241 251 e3t(:,:,:,Krhs) = e3t(:,:,:,Kbb) 242 252 ! Reconstruction of all vertical scale factors at now and before time-steps … … 267 277 END DO 268 278 ENDIF 279 #endif 280 !!st9 269 281 ENDIF 270 282 ENDIF -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90
r13005 r13151 6 6 !! History : NEMO ! 2020-03 (A. Nasser, G. Madec) Original code from 4.0.2 7 7 !!---------------------------------------------------------------------- 8 8 #if defined key_qco 9 !!---------------------------------------------------------------------- 10 !! 'key_qco' EMPTY MODULE Quasi-Eulerian vertical coordonate 11 !!---------------------------------------------------------------------- 12 #else 9 13 !!---------------------------------------------------------------------- 10 14 !! stp : Shallow Water time-stepping … … 13 17 USE phycst ! physical constants 14 18 USE usrdef_nam 15 USE lib_mpp ! MPP library16 USE dynvor , ONLY : ln_dynvor_ens_adVO, ln_dynvor_ens_adKE, ln_dynvor_ens_adKEVO, &17 & ln_dynvor_ene_adVO, ln_dynvor_ene_adKE, ln_dynvor_ene_adKEVO18 19 ! 19 20 USE iom ! xIOs server … … 122 123 ! LATERAL PHYSICS 123 124 ! ! eddy diffusivity coeff. 124 IF( l_ldfdyn_time ) CALL ldf_dyn( kstp, Nbb )! eddy viscosity coeff.125 IF( l_ldfdyn_time ) CALL ldf_dyn( kstp, Nbb ) ! eddy viscosity coeff. 125 126 126 127 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 129 130 130 131 CALL ssh_nxt ( kstp, Nbb, Nnn, ssh, Naa ) ! after ssh (includes call to div_hor) 131 132 IF( .NOT.ln_linssh ) CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn, Naa ) ! after vertical scale factors133 134 132 uu(:,:,:,Nrhs) = 0._wp ! set dynamics trends to zero 135 133 vv(:,:,:,Nrhs) = 0._wp 136 134 137 IF( ln_bdy ) CALL bdy_dyn3d_dmp ( kstp, Nbb, uu, vv, Nrhs ) ! bdy damping trends 135 IF( .NOT.ln_linssh ) CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn, Naa ) ! after vertical scale factors 136 137 IF( ln_bdy ) CALL bdy_dyn3d_dmp ( kstp, Nbb, uu, vv, Nrhs ) ! bdy damping trends 138 138 139 139 #if defined key_agrif 140 140 IF(.NOT. Agrif_Root()) & 141 & CALL Agrif_Sponge_dyn ! momentum sponge 142 #endif 141 & CALL Agrif_Sponge_dyn ! momentum sponge 142 #endif 143 CALL dyn_adv( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS 144 145 CALL dyn_vor( kstp, Nnn , uu, vv, Nrhs ) ! vorticity ==> RHS 146 147 CALL dyn_ldf( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! lateral mixing 143 148 144 149 !!an - calcul du gradient de pression horizontal (explicit) … … 148 153 END_3D 149 154 ! 150 151 ! IF( kstp == nit000 .AND. lwp ) THEN152 ! WRITE(numout,*)153 ! WRITE(numout,*) 'step.F90 : classic script used'154 ! WRITE(numout,*) '~~~~~~~'155 ! IF( ln_dynvor_ens_adVO .OR. ln_dynvor_ens_adKE .OR. ln_dynvor_ens_adKEVO &156 ! & .OR. ln_dynvor_ene_adVO .OR. ln_dynvor_ene_adKE .OR. ln_dynvor_ene_adKEVO ) THEN157 ! CALL ctl_stop('STOP','step : alternative direction asked but classis step' )158 ! ENDIF159 ! ENDIF160 !!an161 ! CALL dyn_adv( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS162 !163 ! CALL dyn_vor( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! vorticity ==> RHS164 !165 !!an In dynvor, dynkegAD is called even if not AD, so we keep the same step.F90166 167 CALL dyn_vor( kstp, Nbb, Nnn , uu, vv, Nrhs) ! vorticity ==> RHS168 169 CALL dyn_ldf( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! lateral mixing170 171 155 ! add wind stress forcing and layer linear friction to the RHS 172 156 z1_2rho0 = 0.5_wp * r1_rho0 … … 175 159 & - rn_rfr * uu(ji,jj,jk,Nbb) 176 160 vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn) & 177 & - rn_rfr * vv(ji,jj,jk,Nbb) 161 & - rn_rfr * vv(ji,jj,jk,Nbb) 178 162 END_3D 179 163 !!an … … 182 166 ! Leap-Frog time splitting + Robert-Asselin time filter on u,v,e3 183 167 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 184 168 169 !! what about IF( .NOT.ln_linssh ) ? 185 170 !!an futur module dyn_nxt (a la place de dyn_atf) 186 171 … … 209 194 uu(ji,jj,jk,Naa) = zua 210 195 vv(ji,jj,jk,Naa) = zva 211 END_3D196 END_3D 212 197 ENDIF 213 198 ! … … 220 205 zue3a = zue3b + rDt * e3u(ji,jj,jk,Nrhs) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk) 221 206 zve3a = zve3b + rDt * e3v(ji,jj,jk,Nrhs) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk) 222 ! 207 ! 223 208 uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa) 224 vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa) 225 END_3D226 ELSE ! Leap Frog time stepping + Asselin filter 209 vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa) 210 END_3D 211 ELSE ! Leap Frog time stepping + Asselin filter 227 212 DO_3D_11_11(1,jpkm1) 228 213 zue3n = e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nnn) … … 239 224 ! ! Asselin time filter on u,v (Nnn) 240 225 uu(ji,jj,jk,Nnn) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_tf 241 vv(ji,jj,jk,Nnn) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_tf 226 vv(ji,jj,jk,Nnn) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_tf 242 227 ! 243 228 e3u(ji,jj,jk,Nnn) = ze3u_tf … … 246 231 ! 247 232 uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa) 248 vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa) 249 END_3D233 vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa) 234 END_3D 250 235 ENDIF 251 236 ENDIF 252 237 238 253 239 CALL lbc_lnk_multi( 'stp', uu(:,:,:,Nnn), 'U', -1., vv(:,:,:,Nnn), 'V', -1., & !* local domain boundaries 254 240 & uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1. ) … … 263 249 !! CALL dyn_atf ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v ) ! time filtering of "now" velocities and scale factors 264 250 !!an TO BE ADDED : a simplifier 265 ! CALL ssh_atf ( kstp, Nbb, Nnn, Naa, ssh ) ! time filtering of "now" sea surface height 266 251 !! CALL ssh_atf ( kstp, Nbb, Nnn, Naa, ssh ) ! time filtering of "now" sea surface height 267 252 IF ( .NOT.( l_1st_euler ) ) THEN ! Only do time filtering for leapfrog timesteps 268 253 ! ! filtering "now" field 269 254 ssh(:,:,Nnn) = ssh(:,:,Nnn) + rn_atfp * ( ssh(:,:,Nbb) - 2 * ssh(:,:,Nnn) + ssh(:,:,Naa) ) 270 255 ENDIF 271 272 256 !!an 273 257 … … 280 264 ! 281 265 CALL dom_vvl_sf_update( kstp, Nbb, Nnn, Naa ) ! recompute vertical scale factors 282 283 266 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 284 267 ! diagnostics and outputs … … 287 270 IF( ln_diacfl ) CALL dia_cfl ( kstp, Nnn ) ! Courant number diagnostics 288 271 289 CALL dia_wri ( kstp, Nnn ) ! ocean model: outputs 290 272 CALL dia_wri ( kstp, Nnn ) ! ocean model: outputs 291 273 ! 292 274 IF( lrst_oce ) CALL rst_write ( kstp, Nbb, Nnn ) ! write output ocean restart file … … 335 317 ! 336 318 END SUBROUTINE stp 319 #endif 337 320 ! 338 321 !!====================================================================== -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/stpctl.F90
r12614 r13151 35 35 INTEGER :: idrun, idtime, idssh, idu, ids1, ids2, idt1, idt2, idc1, idw1, istatus 36 36 LOGICAL :: lsomeoce 37 !!stoops 38 # include "domzgr_substitute.h90" 37 39 !!---------------------------------------------------------------------- 38 40 !! NEMO/OCE 4.0 , NEMO Consortium (2018)
Note: See TracChangeset
for help on using the changeset viewer.
