Changeset 5254 for branches/2015/dev_r5224_CNRS_OPA_Hollingsworth
- Timestamp:
- 2015-05-05T17:15:46+02:00 (9 years ago)
- Location:
- branches/2015/dev_r5224_CNRS_OPA_Hollingsworth/NEMOGCM
- Files:
-
- 3 edited
-
CONFIG/SHARED/namelist_ref (modified) (1 diff)
-
NEMO/OPA_SRC/DYN/dynadv.F90 (modified) (9 diffs)
-
NEMO/OPA_SRC/DYN/dynkeg.F90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5224_CNRS_OPA_Hollingsworth/NEMOGCM/CONFIG/SHARED/namelist_ref
r5190 r5254 803 803 !----------------------------------------------------------------------- 804 804 ln_dynadv_vec = .true. ! vector form (T) or flux form (F) 805 nn_dynkeg = 0 ! scheme for grad(KE): =0 C2 ; =1 Hollingsworth correction 805 806 ln_dynadv_cen2= .false. ! flux form - 2nd order centered scheme 806 807 ln_dynadv_ubs = .false. ! flux form - 3rd order UBS scheme -
branches/2015/dev_r5224_CNRS_OPA_Hollingsworth/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
r5120 r5254 5 5 !!============================================================================== 6 6 !! History : 1.0 ! 2006-11 (G. Madec) Original code 7 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 7 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 8 !! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option 8 9 !!---------------------------------------------------------------------- 9 10 … … 17 18 USE dynkeg ! kinetic energy gradient (dyn_keg routine) 18 19 USE dynzad ! vertical advection (dyn_zad routine) 20 ! 19 21 USE in_out_manager ! I/O manager 20 22 USE lib_mpp ! MPP library … … 25 27 26 28 PUBLIC dyn_adv ! routine called by step module 27 PUBLIC dyn_adv_init ! routine called by opa module29 PUBLIC dyn_adv_init ! routine called by opa module 28 30 31 ! !* namdyn_adv namelist * 29 32 LOGICAL, PUBLIC :: ln_dynadv_vec !: vector form flag 33 INTEGER, PUBLIC :: nn_dynkeg !: scheme of kinetic energy gradient: =0 C2 ; =1 Hollingsworth 30 34 LOGICAL, PUBLIC :: ln_dynadv_cen2 !: flux form - 2nd order centered scheme flag 31 35 LOGICAL, PUBLIC :: ln_dynadv_ubs !: flux form - 3rd order UBS scheme flag … … 38 42 # include "vectopt_loop_substitute.h90" 39 43 !!---------------------------------------------------------------------- 40 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)44 !! NEMO/OPA 3.6 , NEMO Consortium (2015) 41 45 !! $Id$ 42 46 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 63 67 SELECT CASE ( nadv ) ! compute advection trend and add it to general trend 64 68 CASE ( 0 ) 65 CALL dyn_keg ( kt ) ! vector form : horizontal gradient of kinetic energy66 CALL dyn_zad ( kt ) ! vector form : vertical advection69 CALL dyn_keg ( kt, nn_dynkeg ) ! vector form : horizontal gradient of kinetic energy 70 CALL dyn_zad ( kt ) ! vector form : vertical advection 67 71 CASE ( 1 ) 68 CALL dyn_keg ( kt ) ! vector form : horizontal gradient of kinetic energy69 CALL dyn_zad_zts ( kt ) ! vector form : vertical advection with sub-timestepping72 CALL dyn_keg ( kt, nn_dynkeg ) ! vector form : horizontal gradient of kinetic energy 73 CALL dyn_zad_zts ( kt ) ! vector form : vertical advection with sub-timestepping 70 74 CASE ( 2 ) 71 CALL dyn_adv_cen2( kt ) ! 2nd order centered scheme75 CALL dyn_adv_cen2( kt ) ! 2nd order centered scheme 72 76 CASE ( 3 ) 73 CALL dyn_adv_ubs ( kt ) ! 3rd order UBS scheme77 CALL dyn_adv_ubs ( kt ) ! 3rd order UBS scheme 74 78 ! 75 CASE (-1 ) ! esopa: test all possibility with control print76 CALL dyn_keg ( kt )79 CASE (-1 ) ! esopa: test all possibility with control print 80 CALL dyn_keg ( kt, nn_dynkeg ) 77 81 CALL dyn_zad ( kt ) 78 82 CALL dyn_adv_cen2( kt ) … … 92 96 !! momentum advection formulation & scheme and set nadv 93 97 !!---------------------------------------------------------------------- 94 INTEGER :: ioptio 95 INTEGER :: ios ! Local integer output status for namelist read 96 !! 97 NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts 98 INTEGER :: ioptio, ios ! Local integer 99 ! 100 NAMELIST/namdyn_adv/ ln_dynadv_vec, nn_dynkeg, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts 98 101 !!---------------------------------------------------------------------- 99 102 ! 100 103 REWIND( numnam_ref ) ! Namelist namdyn_adv in reference namelist : Momentum advection scheme 101 104 READ ( numnam_ref, namdyn_adv, IOSTAT = ios, ERR = 901) … … 112 115 WRITE(numout,*) '~~~~~~~~~~~' 113 116 WRITE(numout,*) ' Namelist namdyn_adv : chose a advection formulation & scheme for momentum' 114 WRITE(numout,*) ' Vector/flux form (T/F) ln_dynadv_vec = ', ln_dynadv_vec 115 WRITE(numout,*) ' 2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2 116 WRITE(numout,*) ' 3rd order UBS advection scheme ln_dynadv_ubs = ', ln_dynadv_ubs 117 WRITE(numout,*) ' Sub timestepping of vertical advection ln_dynzad_zts = ', ln_dynzad_zts 117 WRITE(numout,*) ' Vector/flux form (T/F) ln_dynadv_vec = ', ln_dynadv_vec 118 WRITE(numout,*) ' = 0 standard scheme ; =1 Hollingsworth scheme nn_dynkeg = ', nn_dynkeg 119 WRITE(numout,*) ' 2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2 120 WRITE(numout,*) ' 3rd order UBS advection scheme ln_dynadv_ubs = ', ln_dynadv_ubs 121 WRITE(numout,*) ' Sub timestepping of vertical advection ln_dynzad_zts = ', ln_dynzad_zts 118 122 ENDIF 119 123 … … 126 130 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' ) 127 131 IF( ln_dynzad_zts .AND. .NOT. ln_dynadv_vec ) & 128 129 IF( ln_dynzad_zts .AND. ln_isfcav ) &130 CALL ctl_stop( 'Sub timestepping of vertical advection does not work with ln_isfcav = .TRUE.' )132 CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' ) 133 IF( nn_dynkeg /= nkeg_C2 .AND. nn_dynkeg /= nkeg_HW ) & 134 CALL ctl_stop( 'KEG scheme wrong value of nn_dynkeg' ) 131 135 132 136 ! ! Set nadv … … 139 143 IF(lwp) THEN ! Print the choice 140 144 WRITE(numout,*) 141 IF( nadv == 0 ) WRITE(numout,*) ' vector form : keg + zad + vor is used' 145 IF( nadv == 0 ) WRITE(numout,*) ' vector form : keg + zad + vor is used' 142 146 IF( nadv == 1 ) WRITE(numout,*) ' vector form : keg + zad_zts + vor is used' 147 IF( nadv == 0 .OR. nadv == 1 ) THEN 148 IF( nn_dynkeg == nkeg_C2 ) WRITE(numout,*) 'with Centered standard keg scheme' 149 IF( nn_dynkeg == nkeg_HW ) WRITE(numout,*) 'with Hollingsworth keg scheme' 150 ENDIF 143 151 IF( nadv == 2 ) WRITE(numout,*) ' flux form : 2nd order scheme is used' 144 152 IF( nadv == 3 ) WRITE(numout,*) ' flux form : UBS scheme is used' -
branches/2015/dev_r5224_CNRS_OPA_Hollingsworth/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
r4990 r5254 4 4 !! Ocean dynamics: kinetic energy gradient trend 5 5 !!====================================================================== 6 !! History : 1.0 ! 87-09 (P. Andrich, m.-a. Foujols) Original code 7 !! 7.0 ! 97-05 (G. Madec) Split dynber into dynkeg and dynhpg 8 !! 9.0 ! 02-07 (G. Madec) F90: Free form and module 6 !! History : 1.0 ! 1987-09 (P. Andrich, M.-A. Foujols) Original code 7 !! 7.0 ! 1997-05 (G. Madec) Split dynber into dynkeg and dynhpg 8 !! NEMO 1.0 ! 2002-07 (G. Madec) F90: Free form and module 9 !! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option 9 10 !!---------------------------------------------------------------------- 10 11 … … 18 19 ! 19 20 USE in_out_manager ! I/O manager 21 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 20 22 USE lib_mpp ! MPP library 21 23 USE prtctl ! Print control … … 28 30 PUBLIC dyn_keg ! routine called by step module 29 31 32 INTEGER, PARAMETER, PUBLIC :: nkeg_C2 = 0 !: 2nd order centered scheme (standard scheme) 33 INTEGER, PARAMETER, PUBLIC :: nkeg_HW = 1 !: Hollingsworth et al., QJRMS, 1983 34 ! 35 REAL(wp) :: r1_48 = 1._wp / 48._wp !: =1/(4*2*6) 36 30 37 !! * Substitutions 31 38 # include "vectopt_loop_substitute.h90" 32 39 !!---------------------------------------------------------------------- 33 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)40 !! NEMO/OPA 3.6 , NEMO Consortium (2015) 34 41 !! $Id$ 35 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 37 44 CONTAINS 38 45 39 SUBROUTINE dyn_keg( kt )46 SUBROUTINE dyn_keg( kt, kscheme ) 40 47 !!---------------------------------------------------------------------- 41 48 !! *** ROUTINE dyn_keg *** … … 45 52 !! general momentum trend. 46 53 !! 47 !! ** Method : Compute the now horizontal kinetic energy 54 !! ** Method : * kscheme = nkeg_C2 : 2nd order centered scheme that 55 !! conserve kinetic energy. Compute the now horizontal kinetic energy 48 56 !! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ] 57 !! * kscheme = nkeg_HW : Hollingsworth correction following 58 !! Arakawa (2001). The now horizontal kinetic energy is given by: 59 !! zhke = 1/6 [ mi-1( 2 * un^2 + ((un(j+1)+un(j-1))/2)^2 ) 60 !! + mj-1( 2 * vn^2 + ((vn(i+1)+vn(i-1))/2)^2 ) ] 61 !! 49 62 !! Take its horizontal gradient and add it to the general momentum 50 63 !! trend (ua,va). … … 54 67 !! ** Action : - Update the (ua, va) with the hor. ke gradient trend 55 68 !! - send this trends to trd_dyn (l_trddyn=T) for post-processing 69 !! 70 !! ** References : Arakawa, A., International Geophysics 2001. 71 !! Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983. 56 72 !!---------------------------------------------------------------------- 57 INTEGER, INTENT( in ) :: kt ! ocean time-step index 73 INTEGER, INTENT( in ) :: kt ! ocean time-step index 74 INTEGER, INTENT( in ) :: kscheme ! =0/1 type of KEG scheme 58 75 ! 59 76 INTEGER :: ji, jj, jk ! dummy loop indices … … 63 80 !!---------------------------------------------------------------------- 64 81 ! 65 IF( nn_timing == 1 ) CALL timing_start('dyn_keg')82 IF( nn_timing == 1 ) CALL timing_start('dyn_keg') 66 83 ! 67 CALL wrk_alloc( jpi, jpj, jpk,zhke )84 CALL wrk_alloc( jpi,jpj,jpk, zhke ) 68 85 ! 69 86 IF( kt == nit000 ) THEN 70 87 IF(lwp) WRITE(numout,*) 71 IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend '88 IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme 72 89 IF(lwp) WRITE(numout,*) '~~~~~~~' 73 90 ENDIF 74 91 75 92 IF( l_trddyn ) THEN ! Save ua and va trends 76 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )93 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 77 94 ztrdu(:,:,:) = ua(:,:,:) 78 95 ztrdv(:,:,:) = va(:,:,:) 79 96 ENDIF 80 97 81 ! ! =============== 82 DO jk = 1, jpkm1 ! Horizontal slab 83 ! ! =============== 84 DO jj = 2, jpj ! Horizontal kinetic energy at T-point 85 DO ji = fs_2, jpi ! vector opt. 86 zu = 0.25 * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 87 & + un(ji ,jj ,jk) * un(ji ,jj ,jk) ) 88 zv = 0.25 * ( vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 89 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) 90 zhke(ji,jj,jk) = zv + zu 91 !!gm simplier coding ==>> ~ faster 92 ! don't forget to suppress local zu zv scalars 93 ! zhke(ji,jj,jk) = 0.25 * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 94 ! & + un(ji ,jj ,jk) * un(ji ,jj ,jk) & 95 ! & + vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 96 ! & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) 97 !!gm end <<== 98 END DO 99 END DO 100 DO jj = 2, jpjm1 ! add the gradient of kinetic energy to the general momentum trends 98 zhke(:,:,jpk) = 0._wp 99 100 SELECT CASE ( kscheme ) !== Horizontal kinetic energy at T-point ==! 101 ! 102 CASE ( nkeg_C2 ) !-- Standard scheme --! 103 DO jk = 1, jpkm1 104 DO jj = 2, jpj 105 DO ji = fs_2, jpi ! vector opt. 106 zu = un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 107 & + un(ji ,jj ,jk) * un(ji ,jj ,jk) 108 zv = vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 109 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) 110 zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 111 END DO 112 END DO 113 END DO 114 ! 115 CASE ( nkeg_HW ) !-- Hollingsworth scheme --! 116 DO jk = 1, jpkm1 117 DO jj = 2, jpjm1 118 DO ji = fs_2, jpim1 ! vector opt. 119 zu = 8._wp * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 120 & + un(ji ,jj ,jk) * un(ji ,jj ,jk) ) & 121 & + ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) & 122 & + ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) * ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) 123 ! 124 zv = 8._wp * ( vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 125 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) & 126 & + ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) & 127 & + ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) * ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) 128 zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 129 END DO 130 END DO 131 END DO 132 CALL lbc_lnk( zhke, 'T', 1. ) 133 ! 134 END SELECT 135 ! 136 DO jk = 1, jpkm1 !== grad( KE ) added to the general momentum trends ==! 137 DO jj = 2, jpjm1 101 138 DO ji = fs_2, fs_jpim1 ! vector opt. 102 139 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) … … 104 141 END DO 105 142 END DO 106 !!gm idea to be tested ==>> is it faster on scalar computers ? 107 ! DO jj = 2, jpjm1 ! add the gradient of kinetic energy to the general momentum trends 108 ! DO ji = fs_2, fs_jpim1 ! vector opt. 109 ! ua(ji,jj,jk) = ua(ji,jj,jk) - 0.25 * ( + un(ji+1,jj ,jk) * un(ji+1,jj ,jk) & 110 ! & + vn(ji+1,jj-1,jk) * vn(ji+1,jj-1,jk) & 111 ! & + vn(ji+1,jj ,jk) * vn(ji+1,jj ,jk) & 112 ! ! 113 ! & - un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 114 ! & - vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 115 ! & - vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) / e1u(ji,jj) 116 ! ! 117 ! va(ji,jj,jk) = va(ji,jj,jk) - 0.25 * ( un(ji-1,jj+1,jk) * un(ji-1,jj+1,jk) & 118 ! & + un(ji ,jj+1,jk) * un(ji ,jj+1,jk) & 119 ! & + vn(ji ,jj+1,jk) * vn(ji ,jj+1,jk) & 120 ! ! 121 ! & - un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 122 ! & - un(ji ,jj ,jk) * un(ji ,jj ,jk) & 123 ! & - vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) / e2v(ji,jj) 124 ! END DO 125 ! END DO 126 !!gm en idea <<== 127 ! ! =============== 128 END DO ! End of slab 129 ! ! =============== 130 131 IF( l_trddyn ) THEN ! save the Kinetic Energy trends for diagnostic 143 END DO 144 ! 145 IF( l_trddyn ) THEN ! save the Kinetic Energy trends for diagnostic 132 146 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 133 147 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 134 148 CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt ) 135 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )149 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 136 150 ENDIF 137 151 ! … … 139 153 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 140 154 ! 141 CALL wrk_dealloc( jpi, jpj, jpk,zhke )155 CALL wrk_dealloc( jpi,jpj,jpk, zhke ) 142 156 ! 143 IF( nn_timing == 1 ) CALL timing_stop('dyn_keg')157 IF( nn_timing == 1 ) CALL timing_stop('dyn_keg') 144 158 ! 145 159 END SUBROUTINE dyn_keg
Note: See TracChangeset
for help on using the changeset viewer.
