- Timestamp:
- 2015-10-26T10:08:06+01:00 (9 years ago)
- Location:
- branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN
- Files:
-
- 23 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90
r4147 r5832 17 17 !! 3.3 ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module 18 18 !! - ! 2010-10 (R. Furner, G. Madec) runoff and cla added directly here 19 !! 3.6 ! 2014-11 (P. Mathiot) isf added directly here 19 20 !!---------------------------------------------------------------------- 20 21 … … 25 26 USE oce ! ocean dynamics and tracers 26 27 USE dom_oce ! ocean space and time domain 27 USE sbc_oce, ONLY : ln_rnf 28 USE sbc_oce, ONLY : ln_rnf, nn_isf ! surface boundary condition: ocean 28 29 USE sbcrnf ! river runoff 30 USE sbcisf ! ice shelf 29 31 USE cla ! cross land advection (cla_div routine) 30 32 USE in_out_manager ! I/O manager … … 66 68 !! - compute the now divergence given by : 67 69 !! hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 68 !! correct hdiv with runoff inflow (div_rnf) and cross land flow (div_cla) 70 !! correct hdiv with runoff inflow (div_rnf), ice shelf melting (div_isf) 71 !! and cross land flow (div_cla) 69 72 !! II. vorticity : 70 73 !! - save the curl computed at the previous time-step … … 95 98 ! 96 99 CALL wrk_alloc( jpi , jpj+2, zwu ) 97 CALL wrk_alloc( jpi+4, jpj , zwv, k jstart = -1 )100 CALL wrk_alloc( jpi+4, jpj , zwv, kistart = -1 ) 98 101 ! 99 102 IF( kt == nit000 ) THEN … … 225 228 ! ! =============== 226 229 227 IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 228 IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 ) CALL cla_div ( kt ) ! Cross Land Advection (Update Hor. divergence) 230 IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 231 IF( ln_divisf .AND. (nn_isf /= 0) ) CALL sbc_isf_div( hdivn ) ! ice shelf (update hdivn field) 232 IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (Update Hor. divergence) 229 233 230 234 ! 4. Lateral boundary conditions on hdivn and rotn … … 233 237 ! 234 238 CALL wrk_dealloc( jpi , jpj+2, zwu ) 235 CALL wrk_dealloc( jpi+4, jpj , zwv, k jstart = -1 )239 CALL wrk_dealloc( jpi+4, jpj , zwv, kistart = -1 ) 236 240 ! 237 241 IF( nn_timing == 1 ) CALL timing_stop('div_cur') … … 323 327 ! ! =============== 324 328 325 IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 329 IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 330 IF( ln_divisf .AND. (nn_isf .GT. 0) ) CALL sbc_isf_div( hdivn ) ! ice shelf (update hdivn field) 326 331 IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field) 327 332 ! -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
r4624 r5832 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 36 LOGICAL, PUBLIC :: ln_dynzad_zts !: vertical advection with sub-timestepping (requires vector form) 32 37 33 38 INTEGER :: nadv ! choice of the formulation and scheme for the advection … … 37 42 # include "vectopt_loop_substitute.h90" 38 43 !!---------------------------------------------------------------------- 39 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)44 !! NEMO/OPA 3.6 , NEMO Consortium (2015) 40 45 !! $Id$ 41 46 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 62 67 SELECT CASE ( nadv ) ! compute advection trend and add it to general trend 63 68 CASE ( 0 ) 64 CALL dyn_keg ( kt ) ! vector form : horizontal gradient of kinetic energy 65 CALL dyn_zad ( kt ) ! vector form : vertical advection 66 CASE ( 1 ) 67 CALL dyn_adv_cen2( kt ) ! 2nd order centered scheme 68 CASE ( 2 ) 69 CALL dyn_adv_ubs ( kt ) ! 3rd order UBS scheme 69 CALL dyn_keg ( kt, nn_dynkeg ) ! vector form : horizontal gradient of kinetic energy 70 CALL dyn_zad ( kt ) ! vector form : vertical advection 71 CASE ( 1 ) 72 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 74 CASE ( 2 ) 75 CALL dyn_adv_cen2( kt ) ! 2nd order centered scheme 76 CASE ( 3 ) 77 CALL dyn_adv_ubs ( kt ) ! 3rd order UBS scheme 70 78 ! 71 CASE (-1 ) ! esopa: test all possibility with control print72 CALL dyn_keg ( kt )79 CASE (-1 ) ! esopa: test all possibility with control print 80 CALL dyn_keg ( kt, nn_dynkeg ) 73 81 CALL dyn_zad ( kt ) 74 82 CALL dyn_adv_cen2( kt ) … … 88 96 !! momentum advection formulation & scheme and set nadv 89 97 !!---------------------------------------------------------------------- 90 INTEGER :: ioptio 91 INTEGER :: ios ! Local integer output status for namelist read 92 !! 93 NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs 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 94 101 !!---------------------------------------------------------------------- 95 102 ! 96 103 REWIND( numnam_ref ) ! Namelist namdyn_adv in reference namelist : Momentum advection scheme 97 104 READ ( numnam_ref, namdyn_adv, IOSTAT = ios, ERR = 901) … … 108 115 WRITE(numout,*) '~~~~~~~~~~~' 109 116 WRITE(numout,*) ' Namelist namdyn_adv : chose a advection formulation & scheme for momentum' 110 WRITE(numout,*) ' Vector/flux form (T/F) ln_dynadv_vec = ', ln_dynadv_vec 111 WRITE(numout,*) ' 2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2 112 WRITE(numout,*) ' 3rd order UBS advection scheme ln_dynadv_ubs = ', ln_dynadv_ubs 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 113 122 ENDIF 114 123 … … 120 129 121 130 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' ) 131 IF( ln_dynzad_zts .AND. .NOT. ln_dynadv_vec ) & 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' ) 122 135 123 136 ! ! Set nadv 124 137 IF( ln_dynadv_vec ) nadv = 0 125 IF( ln_dynadv_cen2 ) nadv = 1 126 IF( ln_dynadv_ubs ) nadv = 2 138 IF( ln_dynzad_zts ) nadv = 1 139 IF( ln_dynadv_cen2 ) nadv = 2 140 IF( ln_dynadv_ubs ) nadv = 3 127 141 IF( lk_esopa ) nadv = -1 128 142 129 143 IF(lwp) THEN ! Print the choice 130 144 WRITE(numout,*) 131 IF( nadv == 0 ) WRITE(numout,*) ' vector form : keg + zad + vor is used' 132 IF( nadv == 1 ) WRITE(numout,*) ' flux form : 2nd order scheme is used' 133 IF( nadv == 2 ) WRITE(numout,*) ' flux form : UBS scheme is used' 145 IF( nadv == 0 ) WRITE(numout,*) ' vector form : keg + zad + vor is used' 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 151 IF( nadv == 2 ) WRITE(numout,*) ' flux form : 2nd order scheme is used' 152 IF( nadv == 3 ) WRITE(numout,*) ' flux form : UBS scheme is used' 134 153 IF( nadv == -1 ) WRITE(numout,*) ' esopa test: use all advection formulation' 135 154 ENDIF -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90
r3294 r5832 15 15 USE oce ! ocean dynamics and tracers 16 16 USE dom_oce ! ocean space and time domain 17 USE trdmod_oce ! ocean variables trends 18 USE trdmod ! ocean dynamics trends 17 USE trd_oce ! trends: ocean variables 18 USE trddyn ! trend manager: dynamics 19 ! 19 20 USE in_out_manager ! I/O manager 20 21 USE lib_mpp ! MPP library 21 22 USE prtctl ! Print control 22 USE wrk_nemo 23 USE timing 23 USE wrk_nemo ! Memory Allocation 24 USE timing ! Timing 24 25 25 26 IMPLICIT NONE … … 103 104 zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 104 105 zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 105 CALL trd_ mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt )106 CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt ) 106 107 zfu_t(:,:,:) = ua(:,:,:) 107 108 zfv_t(:,:,:) = va(:,:,:) … … 153 154 zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 154 155 zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 155 CALL trd_ mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt )156 CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 156 157 ENDIF 157 158 ! ! Control print -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90
r4153 r5832 16 16 USE oce ! ocean dynamics and tracers 17 17 USE dom_oce ! ocean space and time domain 18 USE trdmod ! ocean dynamics trends 19 USE trdmod_oce ! ocean variables trends 18 USE trd_oce ! trends: ocean variables 19 USE trddyn ! trend manager: dynamics 20 ! 20 21 USE in_out_manager ! I/O manager 21 22 USE prtctl ! Print control 22 23 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 23 24 USE lib_mpp ! MPP library 24 USE wrk_nemo 25 USE timing 25 USE wrk_nemo ! Memory Allocation 26 USE timing ! Timing 26 27 27 28 IMPLICIT NONE … … 115 116 DO jj = 2, jpjm1 ! laplacian 116 117 DO ji = fs_2, fs_jpim1 ! vector opt. 117 zlu_uu(ji,jj,jk,1) = ( ub (ji+1,jj,jk)-2.*ub (ji,jj,jk)+ub (ji-1,jj,jk) ) * umask(ji,jj,jk) 118 zlv_vv(ji,jj,jk,1) = ( vb (ji,jj+1,jk)-2.*vb (ji,jj,jk)+vb (ji,jj-1,jk) ) * vmask(ji,jj,jk) 119 zlu_uv(ji,jj,jk,1) = ( ub (ji,jj+1,jk)-2.*ub (ji,jj,jk)+ub (ji,jj-1,jk) ) * umask(ji,jj,jk) 120 zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj,jk)-2.*vb (ji,jj,jk)+vb (ji-1,jj,jk) ) * vmask(ji,jj,jk) 121 ! 122 zlu_uu(ji,jj,jk,2) = ( zfu(ji+1,jj,jk)-2.*zfu(ji,jj,jk)+zfu(ji-1,jj,jk) ) * umask(ji,jj,jk) 123 zlv_vv(ji,jj,jk,2) = ( zfv(ji,jj+1,jk)-2.*zfv(ji,jj,jk)+zfv(ji,jj-1,jk) ) * vmask(ji,jj,jk) 124 zlu_uv(ji,jj,jk,2) = ( zfu(ji,jj+1,jk)-2.*zfu(ji,jj,jk)+zfu(ji,jj-1,jk) ) * umask(ji,jj,jk) 125 zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj,jk)-2.*zfv(ji,jj,jk)+zfv(ji-1,jj,jk) ) * vmask(ji,jj,jk) 126 END DO 127 END DO 128 END DO 129 !!gm BUG !!! just below this should be +1 in all the communications 130 ! CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', -1.) ; CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', -1.) 131 ! CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', -1.) ; CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', -1.) 132 ! CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', -1.) ; CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', -1.) 133 ! CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', -1.) ; CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', -1.) 134 ! 135 !!gm corrected: 118 ! 119 zlu_uu(ji,jj,jk,1) = ( ub (ji+1,jj ,jk) - 2.*ub (ji,jj,jk) + ub (ji-1,jj ,jk) ) * umask(ji,jj,jk) 120 zlv_vv(ji,jj,jk,1) = ( vb (ji ,jj+1,jk) - 2.*vb (ji,jj,jk) + vb (ji ,jj-1,jk) ) * vmask(ji,jj,jk) 121 zlu_uv(ji,jj,jk,1) = ( ub (ji ,jj+1,jk) - ub (ji ,jj ,jk) ) * fmask(ji ,jj ,jk) & 122 & - ( ub (ji ,jj ,jk) - ub (ji ,jj-1,jk) ) * fmask(ji ,jj-1,jk) 123 zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj ,jk) - vb (ji ,jj ,jk) ) * fmask(ji ,jj ,jk) & 124 & - ( vb (ji ,jj ,jk) - vb (ji-1,jj ,jk) ) * fmask(ji-1,jj ,jk) 125 ! 126 zlu_uu(ji,jj,jk,2) = ( zfu(ji+1,jj ,jk) - 2.*zfu(ji,jj,jk) + zfu(ji-1,jj ,jk) ) * umask(ji,jj,jk) 127 zlv_vv(ji,jj,jk,2) = ( zfv(ji ,jj+1,jk) - 2.*zfv(ji,jj,jk) + zfv(ji ,jj-1,jk) ) * vmask(ji,jj,jk) 128 zlu_uv(ji,jj,jk,2) = ( zfu(ji ,jj+1,jk) - zfu(ji ,jj ,jk) ) * fmask(ji ,jj ,jk) & 129 & - ( zfu(ji ,jj ,jk) - zfu(ji ,jj-1,jk) ) * fmask(ji ,jj-1,jk) 130 zlv_vu(ji,jj,jk,2) = ( zfv(ji+1,jj ,jk) - zfv(ji ,jj ,jk) ) * fmask(ji ,jj ,jk) & 131 & - ( zfv(ji ,jj ,jk) - zfv(ji-1,jj ,jk) ) * fmask(ji-1,jj ,jk) 132 END DO 133 END DO 134 END DO 136 135 CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', 1. ) ; CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', 1. ) 137 136 CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', 1. ) ; CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', 1. ) 138 137 CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. ) ; CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. ) 139 138 CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. ) ; CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. ) 140 !!gm end141 139 142 140 ! ! ====================== ! … … 196 194 zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 197 195 zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 198 CALL trd_ mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt )196 CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt ) 199 197 zfu_t(:,:,:) = ua(:,:,:) 200 198 zfv_t(:,:,:) = va(:,:,:) … … 245 243 zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 246 244 zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 247 CALL trd_ mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt )245 CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 248 246 ENDIF 249 247 ! ! Control print -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
r3294 r5832 10 10 11 11 !!---------------------------------------------------------------------- 12 !! dyn_bfr : Update the momentum trend with the bottom friction contribution12 !! dyn_bfr : Update the momentum trend with the bottom friction contribution 13 13 !!---------------------------------------------------------------------- 14 USE oce 15 USE dom_oce 16 USE zdf_oce 17 USE zdfbfr 18 USE trd mod ! ocean active dynamics and tracers trends19 USE trd mod_oce ! ocean variables trends20 USE in_out_manager 21 USE prtctl 22 USE timing 23 USE wrk_nemo 14 USE oce ! ocean dynamics and tracers variables 15 USE dom_oce ! ocean space and time domain variables 16 USE zdf_oce ! ocean vertical physics variables 17 USE zdfbfr ! ocean bottom friction variables 18 USE trd_oce ! trends: ocean variables 19 USE trddyn ! trend manager: dynamics 20 USE in_out_manager ! I/O manager 21 USE prtctl ! Print control 22 USE timing ! Timing 23 USE wrk_nemo ! Memory Allocation 24 24 25 25 IMPLICIT NONE 26 26 PRIVATE 27 27 28 PUBLIC dyn_bfr 28 PUBLIC dyn_bfr ! routine called by step.F90 29 29 30 30 !! * Substitutions … … 57 57 IF( nn_timing == 1 ) CALL timing_start('dyn_bfr') 58 58 ! 59 !!gm issue: better to put the logical in step to control the call of zdf_bfr 60 !! ==> change the logical from ln_bfrimp to ln_bfr_exp !! 59 61 IF( .NOT.ln_bfrimp) THEN ! only for explicit bottom friction form 60 62 ! implicit bfr is implemented in dynzdf_imp 61 63 64 !!gm bug : time step is only rdt (not 2 rdt if euler start !) 62 65 zm1_2dt = - 1._wp / ( 2._wp * rdt ) 63 66 … … 69 72 70 73 71 # if defined key_vectopt_loop72 DO jj = 1, 173 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)74 # else75 74 DO jj = 2, jpjm1 76 75 DO ji = 2, jpim1 77 # endif78 76 ikbu = mbku(ji,jj) ! deepest ocean u- & v-levels 79 77 ikbv = mbkv(ji,jj) … … 84 82 END DO 85 83 END DO 84 85 IF ( ln_isfcav ) THEN 86 DO jj = 2, jpjm1 87 DO ji = 2, jpim1 88 ! (ISF) stability criteria for top friction 89 ikbu = miku(ji,jj) ! first wet ocean u- & v-levels 90 ikbv = mikv(ji,jj) 91 ! 92 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 93 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( tfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu) & 94 & * (1.-umask(ji,jj,1)) 95 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( tfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv) & 96 & * (1.-vmask(ji,jj,1)) 97 ! (ISF) 98 END DO 99 END DO 100 END IF 86 101 87 102 ! … … 89 104 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 90 105 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 91 CALL trd_ mod( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_trd_bfr, 'DYN', kt )106 CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt ) 92 107 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 93 108 ENDIF -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r4624 r5832 16 16 !! 3.4 ! 2011-11 (H. Liu) hpg_prj: Original code for s-coordinates 17 17 !! ! (A. Coward) suppression of hel, wdj and rot options 18 !! 3.6 ! 2014-11 (P. Mathiot) hpg_isf: original code for ice shelf cavity 18 19 !!---------------------------------------------------------------------- 19 20 … … 25 26 !! hpg_zps : z-coordinate plus partial steps (interpolation) 26 27 !! hpg_sco : s-coordinate (standard jacobian formulation) 28 !! hpg_isf : s-coordinate (sco formulation) adapted to ice shelf 27 29 !! hpg_djc : s-coordinate (Density Jacobian with Cubic polynomial) 28 30 !! hpg_prj : s-coordinate (Pressure Jacobian with Cubic polynomial) 29 31 !!---------------------------------------------------------------------- 30 32 USE oce ! ocean dynamics and tracers 33 USE sbc_oce ! surface variable (only for the flag with ice shelf) 31 34 USE dom_oce ! ocean space and time domain 32 35 USE phycst ! physical constants 33 USE trdmod ! ocean dynamics trends 34 USE trdmod_oce ! ocean variables trends 36 USE trd_oce ! trends: ocean variables 37 USE trddyn ! trend manager: dynamics 38 ! 35 39 USE in_out_manager ! I/O manager 36 40 USE prtctl ! Print control 37 USE lbclnk ! lateral boundary condition 41 USE lbclnk ! lateral boundary condition 38 42 USE lib_mpp ! MPP library 43 USE eosbn2 ! compute density 39 44 USE wrk_nemo ! Memory Allocation 40 45 USE timing ! Timing … … 52 57 LOGICAL , PUBLIC :: ln_hpg_djc !: s-coordinate (Density Jacobian with Cubic polynomial) 53 58 LOGICAL , PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) 59 LOGICAL , PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf 54 60 LOGICAL , PUBLIC :: ln_dynhpg_imp !: semi-implicite hpg flag 55 61 … … 74 80 !! 75 81 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 76 !! - Save the trend(l_trddyn=T)82 !! - send trends to trd_dyn for futher diagnostics (l_trddyn=T) 77 83 !!---------------------------------------------------------------------- 78 84 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 94 100 CASE ( 3 ) ; CALL hpg_djc ( kt ) ! s-coordinate (Density Jacobian with Cubic polynomial) 95 101 CASE ( 4 ) ; CALL hpg_prj ( kt ) ! s-coordinate (Pressure Jacobian scheme) 102 CASE ( 5 ) ; CALL hpg_isf ( kt ) ! s-coordinate similar to sco modify for ice shelf 96 103 END SELECT 97 104 ! … … 99 106 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 100 107 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 101 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt )108 CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt ) 102 109 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 103 110 ENDIF … … 125 132 !! 126 133 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 127 & ln_hpg_djc, ln_hpg_prj, ln_ dynhpg_imp134 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf, ln_dynhpg_imp 128 135 !!---------------------------------------------------------------------- 129 136 ! … … 145 152 WRITE(numout,*) ' z-coord. - partial steps (interpolation) ln_hpg_zps = ', ln_hpg_zps 146 153 WRITE(numout,*) ' s-coord. (standard jacobian formulation) ln_hpg_sco = ', ln_hpg_sco 154 WRITE(numout,*) ' s-coord. (standard jacobian formulation) for isf ln_hpg_isf = ', ln_hpg_isf 147 155 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 148 156 WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj … … 155 163 & either ln_hpg_sco or ln_hpg_prj instead') 156 164 ! 157 IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj ) ) &165 IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) ) & 158 166 & CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:& 159 167 & the standard jacobian formulation hpg_sco or & 160 168 & the pressure jacobian formulation hpg_prj') 169 170 IF( ln_hpg_isf .AND. .NOT. ln_isfcav ) & 171 & CALL ctl_stop( ' hpg_isf not available if ln_isfcav = false ' ) 172 IF( .NOT. ln_hpg_isf .AND. ln_isfcav ) & 173 & CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' ) 161 174 ! 162 175 ! ! Set nhpg from ln_hpg_... flags … … 166 179 IF( ln_hpg_djc ) nhpg = 3 167 180 IF( ln_hpg_prj ) nhpg = 4 181 IF( ln_hpg_isf ) nhpg = 5 168 182 ! 169 183 ! ! Consistency check … … 174 188 IF( ln_hpg_djc ) ioptio = ioptio + 1 175 189 IF( ln_hpg_prj ) ioptio = ioptio + 1 190 IF( ln_hpg_isf ) ioptio = ioptio + 1 176 191 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 192 ! 193 ! initialisation of ice load 194 riceload(:,:)=0.0 177 195 ! 178 196 END SUBROUTINE dyn_hpg_init … … 315 333 316 334 ! partial steps correction at the last level (use gru & grv computed in zpshde.F90) 317 # if defined key_vectopt_loop318 jj = 1319 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)320 # else321 335 DO jj = 2, jpjm1 322 336 DO ji = 2, jpim1 323 # endif324 337 iku = mbku(ji,jj) 325 338 ikv = mbkv(ji,jj) … … 338 351 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 339 352 ENDIF 340 # if ! defined key_vectopt_loop 341 END DO 342 # endif 353 END DO 343 354 END DO 344 355 ! … … 346 357 ! 347 358 END SUBROUTINE hpg_zps 348 349 359 350 360 SUBROUTINE hpg_sco( kt ) … … 434 444 END SUBROUTINE hpg_sco 435 445 446 SUBROUTINE hpg_isf( kt ) 447 !!--------------------------------------------------------------------- 448 !! *** ROUTINE hpg_sco *** 449 !! 450 !! ** Method : s-coordinate case. Jacobian scheme. 451 !! The now hydrostatic pressure gradient at a given level, jk, 452 !! is computed by taking the vertical integral of the in-situ 453 !! density gradient along the model level from the suface to that 454 !! level. s-coordinates (ln_sco): a corrective term is added 455 !! to the horizontal pressure gradient : 456 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] 457 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ] 458 !! add it to the general momentum trend (ua,va). 459 !! ua = ua - 1/e1u * zhpi 460 !! va = va - 1/e2v * zhpj 461 !! iceload is added and partial cell case are added to the top and bottom 462 !! 463 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 464 !!---------------------------------------------------------------------- 465 INTEGER, INTENT(in) :: kt ! ocean time-step index 466 !! 467 INTEGER :: ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j ! dummy loop indices 468 REAL(wp) :: zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1 ! temporary scalars 469 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj, zrhd 470 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop 471 REAL(wp), POINTER, DIMENSION(:,:) :: ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj 472 !!---------------------------------------------------------------------- 473 ! 474 CALL wrk_alloc( jpi,jpj, 2, ztstop) 475 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 476 CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 477 ! 478 IF( kt == nit000 ) THEN 479 IF(lwp) WRITE(numout,*) 480 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 481 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 482 ENDIF 483 484 ! Local constant initialization 485 zcoef0 = - grav * 0.5_wp 486 ! To use density and not density anomaly 487 ! IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume 488 ! ELSE ; znad = 0._wp ! Fixed volume 489 ! ENDIF 490 znad=1._wp 491 ! iniitialised to 0. zhpi zhpi 492 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 493 494 !================================================================================== 495 !=====Compute iceload and contribution of the half first wet layer ================= 496 !=================================================================================== 497 498 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 499 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 500 501 ! compute density of the water displaced by the ice shelf 502 zrhd = rhd ! save rhd 503 DO jk = 1, jpk 504 zdept(:,:)=gdept_1d(jk) 505 CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk)) 506 END DO 507 WHERE ( tmask(:,:,:) == 1._wp) 508 rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd 509 END WHERE 510 511 ! compute rhd at the ice/oce interface (ice shelf side) 512 CALL eos(ztstop,risfdep,zrhdtop_isf) 513 514 ! compute rhd at the ice/oce interface (ocean side) 515 DO ji=1,jpi 516 DO jj=1,jpj 517 ikt=mikt(ji,jj) 518 ztstop(ji,jj,1)=tsn(ji,jj,ikt,1) 519 ztstop(ji,jj,2)=tsn(ji,jj,ikt,2) 520 END DO 521 END DO 522 CALL eos(ztstop,risfdep,zrhdtop_oce) 523 ! 524 ! Surface value + ice shelf gradient 525 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 526 ziceload = 0._wp 527 DO jj = 1, jpj 528 DO ji = 1, jpi ! vector opt. 529 ikt=mikt(ji,jj) 530 ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 531 DO jk=2,ikt-1 532 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 533 & * (1._wp - tmask(ji,jj,jk)) 534 END DO 535 IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) & 536 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 537 END DO 538 END DO 539 riceload(:,:) = 0.0_wp ; riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5 540 ! compute zp from z=0 to first T wet point (correction due to zps not yet applied) 541 DO jj = 2, jpjm1 542 DO ji = fs_2, fs_jpim1 ! vector opt. 543 ikt=mikt(ji,jj) ; iktp1i=mikt(ji+1,jj); iktp1j=mikt(ji,jj+1) 544 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 545 ! we assume ISF is in isostatic equilibrium 546 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj ,iktp1i) & 547 & * ( 2._wp * znad + rhd(ji+1,jj ,iktp1i) + zrhdtop_oce(ji+1,jj ) ) & 548 & - 0.5_wp * fse3w(ji ,jj ,ikt ) & 549 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) & 550 & + ( ziceload(ji+1,jj) - ziceload(ji,jj)) ) 551 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji ,jj+1,iktp1j) & 552 & * ( 2._wp * znad + rhd(ji ,jj+1,iktp1j) + zrhdtop_oce(ji ,jj+1) ) & 553 & - 0.5_wp * fse3w(ji ,jj ,ikt ) & 554 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) & 555 & + ( ziceload(ji,jj+1) - ziceload(ji,jj) ) ) 556 ! s-coordinate pressure gradient correction (=0 if z coordinate) 557 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 558 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 559 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 560 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 561 ! add to the general momentum trend 562 ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 563 va(ji,jj,1) = va(ji,jj,1) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 564 END DO 565 END DO 566 !================================================================================== 567 !===== Compute partial cell contribution for the top cell ========================= 568 !================================================================================== 569 DO jj = 2, jpjm1 570 DO ji = fs_2, fs_jpim1 ! vector opt. 571 iku = miku(ji,jj) ; 572 zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp 573 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 574 ! u direction 575 IF ( iku .GT. 1 ) THEN 576 ! case iku 577 zhpi(ji,jj,iku) = zcoef0 / e1u(ji,jj) * ze3wu & 578 & * ( rhd (ji+1,jj,iku) + rhd (ji,jj,iku) & 579 & + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad ) 580 ! corrective term ( = 0 if z coordinate ) 581 zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj) 582 ! zhpi will be added in interior loop 583 ua(ji,jj,iku) = ua(ji,jj,iku) + zuap 584 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure 585 IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj) = zhpi(ji,jj,iku) 586 587 ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps) 588 zhpiint = zcoef0 / e1u(ji,jj) & 589 & * ( fse3w(ji+1,jj ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad) & 590 & + (rhd(ji+1,jj,iku ) + znad) ) * tmask(ji+1,jj,iku) & 591 & - fse3w(ji ,jj ,iku+1) * ( (rhd(ji ,jj,iku+1) + znad) & 592 & + (rhd(ji ,jj,iku ) + znad) ) * tmask(ji ,jj,iku) ) 593 zhpi(ji,jj,iku+1) = zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint 594 END IF 595 596 ! v direction 597 ikv = mikv(ji,jj) 598 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 599 IF ( ikv .GT. 1 ) THEN 600 ! case ikv 601 zhpj(ji,jj,ikv) = zcoef0 / e2v(ji,jj) * ze3wv & 602 & * ( rhd(ji,jj+1,ikv) + rhd (ji,jj,ikv) & 603 & + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad ) 604 ! corrective term ( = 0 if z coordinate ) 605 zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj) 606 ! zhpi will be added in interior loop 607 va(ji,jj,ikv) = va(ji,jj,ikv) + zvap 608 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure 609 IF (mbkv(ji,jj) == ikv + 1) zpshpj(ji,jj) = zhpj(ji,jj,ikv) 610 611 ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps) 612 zhpjint = zcoef0 / e2v(ji,jj) & 613 & * ( fse3w(ji ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad) & 614 & + (rhd(ji,jj+1,ikv ) + znad) ) * tmask(ji,jj+1,ikv) & 615 & - fse3w(ji ,jj ,ikv+1) * ( (rhd(ji,jj ,ikv+1) + znad) & 616 & + (rhd(ji,jj ,ikv ) + znad) ) * tmask(ji,jj ,ikv) ) 617 zhpj(ji,jj,ikv+1) = zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint 618 END IF 619 END DO 620 END DO 621 622 !================================================================================== 623 !===== Compute interior value ===================================================== 624 !================================================================================== 625 626 DO jj = 2, jpjm1 627 DO ji = fs_2, fs_jpim1 ! vector opt. 628 iku=miku(ji,jj); ikv=mikv(ji,jj) 629 DO jk = 2, jpkm1 630 ! hydrostatic pressure gradient along s-surfaces 631 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 632 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1) & 633 & + zcoef0 / e1u(ji,jj) & 634 & * ( fse3w(ji+1,jj ,jk) * ( (rhd(ji+1,jj,jk ) + znad) & 635 & + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1) & 636 & - fse3w(ji ,jj ,jk) * ( (rhd(ji ,jj,jk ) + znad) & 637 & + (rhd(ji ,jj,jk-1) + znad) ) * tmask(ji ,jj,jk-1) ) 638 ! s-coordinate pressure gradient correction 639 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 640 zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 641 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1) 642 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 643 644 ! hydrostatic pressure gradient along s-surfaces 645 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 646 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1) & 647 & + zcoef0 / e2v(ji,jj) & 648 & * ( fse3w(ji ,jj+1,jk) * ( (rhd(ji,jj+1,jk ) + znad) & 649 & + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1) & 650 & - fse3w(ji ,jj ,jk) * ( (rhd(ji,jj ,jk ) + znad) & 651 & + (rhd(ji,jj ,jk-1) + znad) ) * tmask(ji,jj ,jk-1) ) 652 ! s-coordinate pressure gradient correction 653 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 654 zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 655 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1) 656 ! add to the general momentum trend 657 va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 658 END DO 659 END DO 660 END DO 661 662 !================================================================================== 663 !===== Compute bottom cell contribution (partial cell) ============================ 664 !================================================================================== 665 666 DO jj = 2, jpjm1 667 DO ji = 2, jpim1 668 iku = mbku(ji,jj) 669 ikv = mbkv(ji,jj) 670 671 IF (iku .GT. 1) THEN 672 ! remove old value (interior case) 673 zuap = -zcoef0 * ( rhd (ji+1,jj ,iku) + rhd (ji,jj,iku) + 2._wp * znad ) & 674 & * ( fsde3w(ji+1,jj ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj) 675 ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 676 ! put new value 677 ! -zpshpi to avoid double contribution of the partial step in the top layer 678 zuap = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj) / e1u(ji,jj) 679 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 680 ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 681 END IF 682 ! v direction 683 IF (ikv .GT. 1) THEN 684 ! remove old value (interior case) 685 zvap = -zcoef0 * ( rhd (ji ,jj+1,ikv) + rhd (ji,jj,ikv) + 2._wp * znad ) & 686 & * ( fsde3w(ji ,jj+1,ikv) - fsde3w(ji,jj,ikv) ) / e2v(ji,jj) 687 va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 688 ! put new value 689 ! -zpshpj to avoid double contribution of the partial step in the top layer 690 zvap = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) / e2v(ji,jj) 691 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj) 692 va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 693 END IF 694 END DO 695 END DO 696 697 ! set back to original density value into the ice shelf cell (maybe useless because it is masked) 698 rhd = zrhd 699 ! 700 CALL wrk_dealloc( jpi,jpj,2, ztstop) 701 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 702 CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 703 ! 704 END SUBROUTINE hpg_isf 705 706 436 707 SUBROUTINE hpg_djc( kt ) 437 708 !!--------------------------------------------------------------------- … … 671 942 !! 672 943 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 673 !! - Save the trend (l_trddyn=T)674 !!675 944 !!---------------------------------------------------------------------- 676 945 INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear … … 687 956 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdept, zrhh 688 957 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 958 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_n, zsshv_n 689 959 !!---------------------------------------------------------------------- 690 960 ! 691 961 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 692 962 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 963 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 693 964 ! 694 965 IF( kt == nit000 ) THEN … … 724 995 725 996 ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)" 726 DO jj = 1, jpj; DO ji = 1, jpi 727 zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad 728 END DO ; END DO 729 730 DO jk = 2, jpk; DO jj = 1, jpj; DO ji = 1, jpi 731 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk) 732 END DO ; END DO ; END DO 733 734 fsp(:,:,:) = zrhh(:,:,:) 997 DO jj = 1, jpj 998 DO ji = 1, jpi 999 zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad 1000 END DO 1001 END DO 1002 1003 DO jk = 2, jpk 1004 DO jj = 1, jpj 1005 DO ji = 1, jpi 1006 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk) 1007 END DO 1008 END DO 1009 END DO 1010 1011 fsp(:,:,:) = zrhh (:,:,:) 735 1012 xsp(:,:,:) = zdept(:,:,:) 736 1013 … … 765 1042 766 1043 ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1) 1044 1045 ! Prepare zsshu_n and zsshv_n 767 1046 DO jj = 2, jpjm1 768 1047 DO ji = 2, jpim1 769 zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshn(ji,jj) * znad) ! probable bug: changed from sshu_n for ztilde compilation 770 zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshn(ji,jj) * znad) ! probable bug: changed from sshv_n for ztilde compilation 1048 zsshu_n(ji,jj) = (e12u(ji,jj) * sshn(ji,jj) + e12u(ji+1, jj) * sshn(ji+1,jj)) * & 1049 & r1_e12u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1050 zsshv_n(ji,jj) = (e12v(ji,jj) * sshn(ji,jj) + e12v(ji+1, jj) * sshn(ji,jj+1)) * & 1051 & r1_e12v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1052 END DO 1053 END DO 1054 1055 DO jj = 2, jpjm1 1056 DO ji = 2, jpim1 1057 zu(ji,jj,1) = - ( fse3u(ji,jj,1) - zsshu_n(ji,jj) * znad) 1058 zv(ji,jj,1) = - ( fse3v(ji,jj,1) - zsshv_n(ji,jj) * znad) 771 1059 END DO 772 1060 END DO … … 930 1218 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 931 1219 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1220 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 932 1221 ! 933 1222 END SUBROUTINE hpg_prj 934 1223 1224 935 1225 SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type) 936 1226 !!---------------------------------------------------------------------- … … 940 1230 !! 941 1231 !! ** Method : f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 1232 !! 942 1233 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 943 !!944 1234 !!---------------------------------------------------------------------- 945 1235 IMPLICIT NONE … … 949 1239 INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline 950 1240 ! 2: Linear 951 952 ! Local Variables 1241 ! 953 1242 INTEGER :: ji, jj, jk ! dummy loop indices 954 1243 INTEGER :: jpi, jpj, jpkm1 … … 1040 1329 ENDIF 1041 1330 1042 1043 1331 END SUBROUTINE cspline 1044 1332 … … 1050 1338 !! ** Purpose : 1-d linear interpolation 1051 1339 !! 1052 !! ** Method : 1053 !! interpolation is straight forward 1340 !! ** Method : interpolation is straight forward 1054 1341 !! extrapolation is also permitted (no value limit) 1055 !!1056 1342 !!---------------------------------------------------------------------- 1057 1343 IMPLICIT NONE … … 1070 1356 END FUNCTION interp1 1071 1357 1358 1072 1359 FUNCTION interp2(x, a, b, c, d) RESULT(f) 1073 1360 !!---------------------------------------------------------------------- … … 1133 1420 END FUNCTION integ_spline 1134 1421 1135 1136 1422 !!====================================================================== 1137 1423 END MODULE dynhpg -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
r3294 r5832 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 … … 14 15 USE oce ! ocean dynamics and tracers 15 16 USE dom_oce ! ocean space and time domain 16 USE trdmod ! ocean dynamics trends 17 USE trdmod_oce ! ocean variables trends 17 USE trd_oce ! trends: ocean variables 18 USE trddyn ! trend manager: dynamics 19 ! 18 20 USE in_out_manager ! I/O manager 21 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 19 22 USE lib_mpp ! MPP library 20 23 USE prtctl ! Print control … … 27 30 PUBLIC dyn_keg ! routine called by step module 28 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 29 37 !! * Substitutions 30 38 # include "vectopt_loop_substitute.h90" 31 39 !!---------------------------------------------------------------------- 32 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)40 !! NEMO/OPA 3.6 , NEMO Consortium (2015) 33 41 !! $Id$ 34 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 36 44 CONTAINS 37 45 38 SUBROUTINE dyn_keg( kt )46 SUBROUTINE dyn_keg( kt, kscheme ) 39 47 !!---------------------------------------------------------------------- 40 48 !! *** ROUTINE dyn_keg *** … … 44 52 !! general momentum trend. 45 53 !! 46 !! ** 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 47 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 !! 48 62 !! Take its horizontal gradient and add it to the general momentum 49 63 !! trend (ua,va). … … 52 66 !! 53 67 !! ** Action : - Update the (ua, va) with the hor. ke gradient trend 54 !! - save this trends (l_trddyn=T) for post-processing 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. 55 72 !!---------------------------------------------------------------------- 56 INTEGER, INTENT( in ) :: kt ! ocean time-step index 57 !! 73 INTEGER, INTENT( in ) :: kt ! ocean time-step index 74 INTEGER, INTENT( in ) :: kscheme ! =0/1 type of KEG scheme 75 ! 58 76 INTEGER :: ji, jj, jk ! dummy loop indices 59 77 REAL(wp) :: zu, zv ! temporary scalars … … 62 80 !!---------------------------------------------------------------------- 63 81 ! 64 IF( nn_timing == 1 ) CALL timing_start('dyn_keg')82 IF( nn_timing == 1 ) CALL timing_start('dyn_keg') 65 83 ! 66 CALL wrk_alloc( jpi, jpj, jpk,zhke )84 CALL wrk_alloc( jpi,jpj,jpk, zhke ) 67 85 ! 68 86 IF( kt == nit000 ) THEN 69 87 IF(lwp) WRITE(numout,*) 70 IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend '88 IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme 71 89 IF(lwp) WRITE(numout,*) '~~~~~~~' 72 90 ENDIF 73 91 74 92 IF( l_trddyn ) THEN ! Save ua and va trends 75 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )93 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 76 94 ztrdu(:,:,:) = ua(:,:,:) 77 95 ztrdv(:,:,:) = va(:,:,:) 78 96 ENDIF 79 97 80 ! ! =============== 81 DO jk = 1, jpkm1 ! Horizontal slab 82 ! ! =============== 83 DO jj = 2, jpj ! Horizontal kinetic energy at T-point 84 DO ji = fs_2, jpi ! vector opt. 85 zu = 0.25 * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 86 & + un(ji ,jj ,jk) * un(ji ,jj ,jk) ) 87 zv = 0.25 * ( vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 88 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) 89 zhke(ji,jj,jk) = zv + zu 90 !!gm simplier coding ==>> ~ faster 91 ! don't forget to suppress local zu zv scalars 92 ! zhke(ji,jj,jk) = 0.25 * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 93 ! & + un(ji ,jj ,jk) * un(ji ,jj ,jk) & 94 ! & + vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 95 ! & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) 96 !!gm end <<== 97 END DO 98 END DO 99 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 100 138 DO ji = fs_2, fs_jpim1 ! vector opt. 101 139 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) … … 103 141 END DO 104 142 END DO 105 !!gm idea to be tested ==>> is it faster on scalar computers ? 106 ! DO jj = 2, jpjm1 ! add the gradient of kinetic energy to the general momentum trends 107 ! DO ji = fs_2, fs_jpim1 ! vector opt. 108 ! ua(ji,jj,jk) = ua(ji,jj,jk) - 0.25 * ( + un(ji+1,jj ,jk) * un(ji+1,jj ,jk) & 109 ! & + vn(ji+1,jj-1,jk) * vn(ji+1,jj-1,jk) & 110 ! & + vn(ji+1,jj ,jk) * vn(ji+1,jj ,jk) & 111 ! ! 112 ! & - un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 113 ! & - vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 114 ! & - vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) / e1u(ji,jj) 115 ! ! 116 ! va(ji,jj,jk) = va(ji,jj,jk) - 0.25 * ( un(ji-1,jj+1,jk) * un(ji-1,jj+1,jk) & 117 ! & + un(ji ,jj+1,jk) * un(ji ,jj+1,jk) & 118 ! & + vn(ji ,jj+1,jk) * vn(ji ,jj+1,jk) & 119 ! ! 120 ! & - un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 121 ! & - un(ji ,jj ,jk) * un(ji ,jj ,jk) & 122 ! & - vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) / e2v(ji,jj) 123 ! END DO 124 ! END DO 125 !!gm en idea <<== 126 ! ! =============== 127 END DO ! End of slab 128 ! ! =============== 129 130 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 131 146 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 132 147 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 133 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_keg, 'DYN', kt )134 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )148 CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt ) 149 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 135 150 ENDIF 136 151 ! … … 138 153 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 139 154 ! 140 CALL wrk_dealloc( jpi, jpj, jpk,zhke )155 CALL wrk_dealloc( jpi,jpj,jpk, zhke ) 141 156 ! 142 IF( nn_timing == 1 ) CALL timing_stop('dyn_keg')157 IF( nn_timing == 1 ) CALL timing_stop('dyn_keg') 143 158 ! 144 159 END SUBROUTINE dyn_keg -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
r4522 r5832 15 15 USE phycst ! physical constants 16 16 USE ldfdyn_oce ! ocean dynamics lateral physics 17 USE ldftra_oce ! ocean tracers lateral physics 17 18 USE ldfslp ! lateral mixing: slopes of mixing orientation 18 19 USE dynldf_bilapg ! lateral mixing (dyn_ldf_bilapg routine) … … 20 21 USE dynldf_iso ! lateral mixing (dyn_ldf_iso routine) 21 22 USE dynldf_lap ! lateral mixing (dyn_ldf_lap routine) 22 USE ldftra_oce, ONLY: ln_traldf_hor ! ocean tracers lateral physics23 USE trd mod ! ocean dynamics and tracer trends24 USE trdmod_oce ! ocean variables trends23 USE trd_oce ! trends: ocean variables 24 USE trddyn ! trend manager: dynamics (trd_dyn routine) 25 ! 25 26 USE prtctl ! Print control 26 27 USE in_out_manager ! I/O manager … … 30 31 USE timing ! Timing 31 32 32 33 33 IMPLICIT NONE 34 34 PRIVATE … … 55 55 !! ** Purpose : compute the lateral ocean dynamics physics. 56 56 !!---------------------------------------------------------------------- 57 !58 57 INTEGER, INTENT(in) :: kt ! ocean time-step index 59 58 ! … … 107 106 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 108 107 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 109 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_ldf, 'DYN', kt )108 CALL trd_dyn( ztrdu, ztrdv, jpdyn_ldf, kt ) 110 109 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 111 110 ENDIF -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90
r3634 r5832 19 19 USE dom_oce ! ocean space and time domain 20 20 USE ldfdyn_oce ! ocean dynamics: lateral physics 21 ! 21 22 USE in_out_manager ! I/O manager 22 USE trdmod ! ocean dynamics trends23 USE trdmod_oce ! ocean variables trends24 23 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 25 24 USE wrk_nemo ! Memory Allocation … … 70 69 !! Add this before trend to the general trend (ua,va): 71 70 !! (ua,va) = (ua,va) + (diffu,diffv) 72 !! 'key_trddyn' defined: the two components of the horizontal73 !! diffusion trend are saved.74 71 !! 75 72 !! ** Action : - Update (ua,va) with the before iso-level biharmonic 76 73 !! mixing trend. 77 74 !!---------------------------------------------------------------------- 78 !79 75 INTEGER, INTENT(in) :: kt ! ocean time-step index 80 76 ! -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90
r4488 r5832 19 19 USE dom_oce ! ocean space and time domain 20 20 USE ldfdyn_oce ! ocean dynamics lateral physics 21 USE zdf_oce ! ocean vertical physics 22 USE ldfslp ! iso-neutral slopes available 21 23 USE ldftra_oce, ONLY: ln_traldf_iso 22 USE zdf_oce ! ocean vertical physics 23 USE trdmod ! ocean dynamics trends 24 USE trdmod_oce ! ocean variables trends 25 USE ldfslp ! iso-neutral slopes available 24 ! 26 25 USE in_out_manager ! I/O manager 27 26 USE lib_mpp ! MPP library … … 81 80 !! -3- Add this trend to the general trend (ta,sa): 82 81 !! (ua,va) = (ua,va) + (zwk3,zwk4) 83 !! 'key_trddyn' defined: the trend is saved for diagnostics.84 82 !! 85 83 !! ** Action : - Update (ua,va) arrays with the before geopotential 86 84 !! biharmonic mixing trend. 87 !! - save the trend in (zwk3,zwk4) ('key_trddyn')88 85 !!---------------------------------------------------------------------- 89 86 INTEGER, INTENT( in ) :: kt ! ocean time-step index … … 201 198 !! pu and pv (all the components except 202 199 !! second order vertical derivative term) 203 !! 'key_trddyn' defined: the trend is saved for diagnostics. 204 !!---------------------------------------------------------------------- 205 !! 200 !!---------------------------------------------------------------------- 206 201 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu , pv ! 1st call: before horizontal velocity 207 202 ! ! 2nd call: ahm x these fields -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
r4488 r5832 22 22 USE ldftra_oce ! ocean tracer lateral physics 23 23 USE zdf_oce ! ocean vertical physics 24 USE trdmod ! ocean dynamics trends25 USE trdmod_oce ! ocean variables trends26 24 USE ldfslp ! iso-neutral slopes 25 ! 27 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 28 27 USE in_out_manager ! I/O manager -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap.F90
r3294 r5832 19 19 USE ldfdyn_oce ! ocean dynamics: lateral physics 20 20 USE zdf_oce ! ocean vertical physics 21 ! 21 22 USE in_out_manager ! I/O manager 22 USE trdmod ! ocean dynamics trends23 USE trdmod_oce ! ocean variables trends24 USE ldfslp ! iso-neutral slopes25 23 USE timing ! Timing 26 24 … … 57 55 !! Add this before trend to the general trend (ua,va): 58 56 !! (ua,va) = (ua,va) + (diffu,diffv) 59 !! 'key_trddyn' activated: the two components of the horizontal60 !! diffusion trend are saved.61 57 !! 62 !! ** Action : - Update (ua,va) with the before iso-level harmonic 63 !! mixing trend. 58 !! ** Action : - Update (ua,va) with the iso-level harmonic mixing trend 64 59 !!---------------------------------------------------------------------- 65 60 INTEGER, INTENT( in ) :: kt ! ocean time-step index -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynnept.F90
- Property svn:keywords set to Id
r4624 r5832 69 69 !!---------------------------------------------------------------------- 70 70 71 !! $Id$ 71 72 CONTAINS 72 73 -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r4370 r5832 18 18 !! 3.3 ! 2011-03 (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL 19 19 !! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes 20 !! 3.7 ! 2014-04 (G. Madec) add the diagnostic of the time filter trends 20 21 !!------------------------------------------------------------------------- 21 22 … … 34 35 USE bdydyn ! ocean open boundary conditions 35 36 USE bdyvol ! ocean open boundary condition (bdy_vol routines) 37 USE trd_oce ! trends: ocean variables 38 USE trddyn ! trend manager: dynamics 39 USE trdken ! trend manager: kinetic energy 40 ! 36 41 USE in_out_manager ! I/O manager 42 USE iom ! I/O manager library 37 43 USE lbclnk ! lateral boundary condition (or mpp link) 38 44 USE lib_mpp ! MPP library 39 45 USE wrk_nemo ! Memory Allocation 40 46 USE prtctl ! Print control 41 47 USE timing ! Timing 42 48 #if defined key_agrif 43 49 USE agrif_opa_interp 44 50 #endif 45 USE timing ! Timing46 51 47 52 IMPLICIT NONE … … 79 84 !! at the local domain boundaries through lbc_lnk call, 80 85 !! at the one-way open boundaries (lk_bdy=T), 81 !! at the AGRIF zoom 86 !! at the AGRIF zoom boundaries (lk_agrif=T) 82 87 !! 83 88 !! * Apply the time filter applied and swap of the dynamics … … 99 104 REAL(wp) :: z2dt ! temporary scalar 100 105 #endif 101 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars102 REAL(wp) :: zve3a, zve3n, zve3b, zvf 103 REAL(wp), POINTER, DIMENSION(:,:) :: zu a, zva104 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3u_f, ze3v_f 106 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars 107 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - 108 REAL(wp), POINTER, DIMENSION(:,:) :: zue, zve 109 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3u_f, ze3v_f, zua, zva 105 110 !!---------------------------------------------------------------------- 106 111 ! 107 IF( nn_timing == 1 ) CALL timing_start('dyn_nxt')108 ! 109 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f)110 IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva)112 IF( nn_timing == 1 ) CALL timing_start('dyn_nxt') 113 ! 114 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 115 IF( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zue, zve ) 111 116 ! 112 117 IF( kt == nit000 ) THEN … … 152 157 153 158 # if defined key_dynspg_ts 159 !!gm IF ( lk_dynspg_ts ) THEN .... 154 160 ! Ensure below that barotropic velocities match time splitting estimate 155 161 ! Compute actual transport and replace it with ts estimate at "after" time step 156 zu a(:,:) = 0._wp157 zv a(:,:) = 0._wp158 DO jk = 1, jpkm1159 zu a(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)160 zv a(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)162 zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 163 zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 164 DO jk = 2, jpkm1 165 zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 166 zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 161 167 END DO 162 168 DO jk = 1, jpkm1 163 ua(:,:,jk) = ( ua(:,:,jk) - zu a(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)164 va(:,:,jk) = ( va(:,:,jk) - zv a(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)169 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 170 va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 165 171 END DO 166 172 … … 175 181 END DO 176 182 ENDIF 183 !!gm ENDIF 177 184 # endif 178 185 … … 195 202 # endif 196 203 #endif 204 205 IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics 206 z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step 207 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1._wp / rdt 208 ! 209 ! ! Kinetic energy and Conversion 210 IF( ln_KE_trd ) CALL trd_dyn( ua, va, jpdyn_ken, kt ) 211 ! 212 IF( ln_dyn_trd ) THEN ! 3D output: total momentum trends 213 zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt 214 zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt 215 CALL iom_put( "utrd_tot", zua ) ! total momentum trends, except the asselin time filter 216 CALL iom_put( "vtrd_tot", zva ) 217 ENDIF 218 ! 219 zua(:,:,:) = un(:,:,:) ! save the now velocity before the asselin filter 220 zva(:,:,:) = vn(:,:,:) ! (caution: there will be a shift by 1 timestep in the 221 ! ! computation of the asselin filter trends) 222 ENDIF 197 223 198 224 ! Time filter and swap of dynamics arrays … … 217 243 DO jj = 1, jpj 218 244 DO ji = 1, jpi 219 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2. e0_wp * un(ji,jj,jk) + ua(ji,jj,jk) )220 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2. e0_wp * vn(ji,jj,jk) + va(ji,jj,jk) )245 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 246 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 221 247 ! 222 248 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity … … 240 266 ! Add volume filter correction: compatibility with tracer advection scheme 241 267 ! => time filter + conservation correction (only at the first level) 242 fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)243 !268 fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) & 269 & -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 244 270 ENDIF 245 271 ! … … 301 327 ! Revert "before" velocities to time split estimate 302 328 ! Doing it here also means that asselin filter contribution is removed 303 zu a(:,:) = 0._wp304 zv a(:,:) = 0._wp305 DO jk = 1, jpkm1306 zu a(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)307 zv a(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)329 zue(:,:) = fse3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 330 zve(:,:) = fse3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1) 331 DO jk = 2, jpkm1 332 zue(:,:) = zue(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 333 zve(:,:) = zve(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 308 334 END DO 309 335 DO jk = 1, jpkm1 310 ub(:,:,jk) = ub(:,:,jk) - (zu a(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)311 vb(:,:,jk) = vb(:,:,jk) - (zv a(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)336 ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk) 337 vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk) 312 338 END DO 313 339 ENDIF … … 327 353 hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 328 354 END DO 329 hur_b(:,:) = umask (:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) )330 hvr_b(:,:) = vmask (:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) )355 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) ) 356 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) ) 331 357 ENDIF 332 358 ! … … 335 361 ! 336 362 DO jk = 1, jpkm1 337 #if defined key_vectopt_loop338 DO jj = 1, 1 !Vector opt. => forced unrolling339 DO ji = 1, jpij340 #else341 363 DO jj = 1, jpj 342 364 DO ji = 1, jpi 343 #endif344 365 un_b(ji,jj) = un_b(ji,jj) + fse3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 345 366 vn_b(ji,jj) = vn_b(ji,jj) + fse3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) … … 358 379 ! 359 380 ! 381 382 IF( l_trddyn ) THEN ! 3D output: asselin filter trends on momentum 383 zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt 384 zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt 385 CALL trd_dyn( zua, zva, jpdyn_atf, kt ) 386 ENDIF 387 ! 360 388 IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & 361 389 & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) 362 390 ! 363 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f)364 IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva)391 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 392 IF( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zue, zve ) 365 393 ! 366 394 IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt') … … 370 398 !!========================================================================= 371 399 END MODULE dynnxt 372 -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r4496 r5832 26 26 USE sbctide 27 27 USE updtide 28 USE trdmod ! ocean dynamics trends 29 USE trdmod_oce ! ocean variables trends 28 USE trd_oce ! trends: ocean variables 29 USE trddyn ! trend manager: dynamics 30 ! 30 31 USE prtctl ! Print control (prt_ctl routine) 31 32 USE in_out_manager ! I/O manager 32 33 USE lib_mpp ! MPP library 33 USE solver 34 USE wrk_nemo 35 USE timing 34 USE solver ! solver initialization 35 USE wrk_nemo ! Memory Allocation 36 USE timing ! Timing 36 37 37 38 … … 163 164 END DO 164 165 END DO 165 END DO 166 END DO 167 168 !!gm add here a call to dyn_trd for ice pressure gradient, the surf pressure trends ???? 169 166 170 ENDIF 167 171 … … 191 195 CASE( 2 ) 192 196 z2dt = 2. * rdt 193 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt197 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 194 198 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / z2dt - ztrdu(:,:,:) 195 199 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / z2dt - ztrdv(:,:,:) 196 200 END SELECT 197 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_spg, 'DYN', kt )201 CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt ) 198 202 ! 199 203 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) … … 246 250 IF( ( ioptio > 1 .AND. .NOT. lk_esopa ) .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) ) & 247 251 & CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' ) 252 IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. ln_isfcav ) & 253 & CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' ) 248 254 ! 249 255 IF( lk_esopa ) nspg = -1 -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90
r4328 r5832 19 19 USE sbc_oce ! surface boundary condition: ocean 20 20 USE phycst ! physical constants 21 ! 21 22 USE in_out_manager ! I/O manager 22 23 USE lib_mpp ! distributed memory computing library -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90
r4328 r5832 13 13 !! - ! 2006-08 (J.Chanut, A.Sellar) Calls to BDY routines. 14 14 !! 3.2 ! 2009-03 (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 15 !! 3.7 ! 2014-04 (F. Roquet, G. Madec) add some trends diag 15 16 !!---------------------------------------------------------------------- 16 17 #if defined key_dynspg_flt || defined key_esopa … … 36 37 USE bdyvol ! ocean open boundary condition (bdy_vol routine) 37 38 USE cla ! cross land advection 39 USE trd_oce ! trends: ocean variables 40 USE trddyn ! trend manager: dynamics 41 ! 38 42 USE in_out_manager ! I/O manager 39 43 USE lib_mpp ! distributed memory computing library … … 43 47 USE iom 44 48 USE lib_fortran 49 USE timing ! Timing 45 50 #if defined key_agrif 46 51 USE agrif_opa_interp 47 52 #endif 48 USE timing ! Timing49 53 50 54 IMPLICIT NONE … … 99 103 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 100 104 !! 101 !! References : Roullet and Madec 1999, JGR.105 !! References : Roullet and Madec, JGR, 2000. 102 106 !!--------------------------------------------------------------------- 103 107 INTEGER, INTENT(in ) :: kt ! ocean time-step index 104 108 INTEGER, INTENT( out) :: kindic ! solver convergence flag (<0 if not converge) 105 ! !109 ! 106 110 INTEGER :: ji, jj, jk ! dummy loop indices 107 111 REAL(wp) :: z2dt, z2dtg, zgcb, zbtd, ztdgu, ztdgv ! local scalars 112 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 113 REAL(wp), POINTER, DIMENSION(:,:) :: zpw 108 114 !!---------------------------------------------------------------------- 109 115 ! 110 116 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_flt') 111 !112 117 ! 113 118 IF( kt == nit000 ) THEN … … 179 184 END DO 180 185 ! 186 IF( l_trddyn ) THEN ! temporary save of spg trends 187 CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 188 DO jk = 1, jpkm1 ! unweighted time stepping 189 DO jj = 2, jpjm1 190 DO ji = fs_2, fs_jpim1 ! vector opt. 191 ztrdu(ji,jj,jk) = spgu(ji,jj) * umask(ji,jj,jk) 192 ztrdv(ji,jj,jk) = spgv(ji,jj) * vmask(ji,jj,jk) 193 END DO 194 END DO 195 END DO 196 CALL trd_dyn( ztrdu, ztrdv, jpdyn_spgexp, kt ) 197 ENDIF 198 ! 181 199 ENDIF 182 200 … … 194 212 DO jj = 2, jpjm1 195 213 DO ji = fs_2, fs_jpim1 ! vector opt. 196 spgu(ji,jj) = 0._wp 197 spgv(ji,jj) = 0._wp 198 END DO 199 END DO 200 201 ! vertical sum 202 !CDIR NOLOOPCHG 203 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 204 DO jk = 1, jpkm1 205 DO ji = 1, jpij 206 spgu(ji,1) = spgu(ji,1) + fse3u_a(ji,1,jk) * ua(ji,1,jk) 207 spgv(ji,1) = spgv(ji,1) + fse3v_a(ji,1,jk) * va(ji,1,jk) 208 END DO 209 END DO 210 ELSE ! No vector opt. 211 DO jk = 1, jpkm1 212 DO jj = 2, jpjm1 213 DO ji = 2, jpim1 214 spgu(ji,jj) = spgu(ji,jj) + fse3u_a(ji,jj,jk) * ua(ji,jj,jk) 215 spgv(ji,jj) = spgv(ji,jj) + fse3v_a(ji,jj,jk) * va(ji,jj,jk) 216 END DO 217 END DO 218 END DO 219 ENDIF 220 221 ! transport: multiplied by the horizontal scale factor 222 DO jj = 2, jpjm1 214 spgu(ji,jj) = fse3u_a(ji,jj,1) * ua(ji,jj,1) 215 spgv(ji,jj) = fse3v_a(ji,jj,1) * va(ji,jj,1) 216 END DO 217 END DO 218 DO jk = 2, jpkm1 ! vertical sum 219 DO jj = 2, jpjm1 220 DO ji = fs_2, fs_jpim1 ! vector opt. 221 spgu(ji,jj) = spgu(ji,jj) + fse3u_a(ji,jj,jk) * ua(ji,jj,jk) 222 spgv(ji,jj) = spgv(ji,jj) + fse3v_a(ji,jj,jk) * va(ji,jj,jk) 223 END DO 224 END DO 225 END DO 226 227 DO jj = 2, jpjm1 ! transport: multiplied by the horizontal scale factor 223 228 DO ji = fs_2, fs_jpim1 ! vector opt. 224 229 spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj) … … 322 327 ENDIF 323 328 #endif 329 330 IF( l_trddyn ) THEN 331 ztrdu(:,:,:) = ua(:,:,:) ! save the after velocity before the filtered SPG 332 ztrdv(:,:,:) = va(:,:,:) 333 ! 334 CALL wrk_alloc( jpi, jpj, zpw ) 335 ! 336 zpw(:,:) = - z2dt * gcx(:,:) 337 CALL iom_put( "ssh_flt" , zpw ) ! output equivalent ssh modification due to implicit filter 338 ! 339 ! ! save surface pressure flux: -pw at z=0 340 zpw(:,:) = - rau0 * grav * sshn(:,:) * wn(:,:,1) * tmask(:,:,1) 341 CALL iom_put( "pw0_exp" , zpw ) 342 zpw(:,:) = wn(:,:,1) 343 CALL iom_put( "w0" , zpw ) 344 zpw(:,:) = rau0 * z2dtg * gcx(:,:) * wn(:,:,1) * tmask(:,:,1) 345 CALL iom_put( "pw0_flt" , zpw ) 346 ! 347 CALL wrk_dealloc( jpi, jpj, zpw ) 348 ! 349 ENDIF 350 324 351 ! Add the trends multiplied by z2dt to the after velocity 325 352 ! ------------------------------------------------------- … … 336 363 END DO 337 364 338 ! write filtered free surface arrays in restart file 339 ! -------------------------------------------------- 340 IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' ) 341 ! 342 ! 343 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_flt') 365 IF( l_trddyn ) THEN ! save the explicit SPG trends for further diagnostics 366 ztrdu(:,:,:) = ( ua(:,:,:) - ztrdu(:,:,:) ) / z2dt 367 ztrdv(:,:,:) = ( va(:,:,:) - ztrdv(:,:,:) ) / z2dt 368 CALL trd_dyn( ztrdu, ztrdv, jpdyn_spgflt, kt ) 369 ! 370 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 371 ENDIF 372 373 IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' ) ! write filtered free surface arrays in restart file 374 ! 375 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_flt') 344 376 ! 345 377 END SUBROUTINE dyn_spg_flt … … 352 384 !! ** Purpose : Read or write filtered free surface arrays in restart file 353 385 !!---------------------------------------------------------------------- 354 INTEGER , INTENT(in) :: kt 355 CHARACTER(len=*), INTENT(in) :: cdrw 386 INTEGER , INTENT(in) :: kt ! ocean time-step 387 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 356 388 !!---------------------------------------------------------------------- 357 389 ! -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r4624 r5832 22 22 USE dom_oce ! ocean space and time domain 23 23 USE sbc_oce ! surface boundary condition: ocean 24 USE sbcisf ! ice shelf variable (fwfisf) 24 25 USE dynspg_oce ! surface pressure gradient variables 25 26 USE phycst ! physical constants … … 44 45 USE agrif_opa_interp ! agrif 45 46 #endif 46 47 #if defined key_asminc 48 USE asminc ! Assimilation increment 49 #endif 47 50 48 51 IMPLICIT NONE … … 76 79 !!---------------------------------------------------------------------- 77 80 !! NEMO/OPA 3.5 , NEMO Consortium (2013) 78 !! $Id : dynspg_ts.F9081 !! $Id$ 79 82 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 80 83 !!---------------------------------------------------------------------- … … 95 98 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) 96 99 97 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &98 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) )100 IF( ln_dynvor_een .or. ln_dynvor_een_old ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 101 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 99 102 100 103 dyn_spg_ts_alloc = MAXVAL(ierr(:)) … … 216 219 ! 217 220 IF ( kt == nit000 .OR. lk_vvl ) THEN 218 IF ( ln_dynvor_een ) THEN 221 IF ( ln_dynvor_een_old ) THEN 222 DO jj = 1, jpjm1 223 DO ji = 1, jpim1 224 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 225 & ht(ji ,jj ) + ht(ji+1,jj ) ) / 4._wp 226 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj) 227 END DO 228 END DO 229 CALL lbc_lnk( zwz, 'F', 1._wp ) 230 zwz(:,:) = ff(:,:) * zwz(:,:) 231 232 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 233 DO jj = 2, jpj 234 DO ji = fs_2, jpi ! vector opt. 235 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 236 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 237 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 238 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 239 END DO 240 END DO 241 ELSE IF ( ln_dynvor_een ) THEN 219 242 DO jj = 1, jpjm1 220 243 DO ji = 1, jpim1 … … 290 313 ! 291 314 DO jk = 1, jpkm1 292 #if defined key_vectopt_loop 293 DO jj = 1, 1 !Vector opt. => forced unrolling 294 DO ji = 1, jpij 295 #else 296 DO jj = 1, jpj 297 DO ji = 1, jpi 298 #endif 299 zu_frc(ji,jj) = zu_frc(ji,jj) + fse3u_n(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 300 zv_frc(ji,jj) = zv_frc(ji,jj) + fse3v_n(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 301 END DO 302 END DO 315 zu_frc(:,:) = zu_frc(:,:) + fse3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 316 zv_frc(:,:) = zv_frc(:,:) + fse3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 303 317 END DO 304 318 ! … … 346 360 END DO 347 361 ! 348 ELSEIF ( ln_dynvor_een ) THEN! enstrophy and energy conserving scheme362 ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN ! enstrophy and energy conserving scheme 349 363 DO jj = 2, jpjm1 350 364 DO ji = fs_2, fs_jpim1 ! vector opt. … … 440 454 ! ! Surface net water flux and rivers 441 455 IF (ln_bt_fw) THEN 442 zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) )456 zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + rdivisf * fwfisf(:,:) ) 443 457 ELSE 444 zssh_frc(:,:) = zraur * z1_2 * (emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)) 458 zssh_frc(:,:) = zraur * z1_2 * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) & 459 & + rdivisf * ( fwfisf(:,:) + fwfisf_b(:,:) ) ) 445 460 ENDIF 446 461 #if defined key_asminc 447 462 ! ! Include the IAU weighted SSH increment 448 463 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 449 zssh_frc(:,:) = zssh_frc(:,:) +ssh_iau(:,:)464 zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:) 450 465 ENDIF 451 466 #endif … … 464 479 ! ! ==================== ! 465 480 ! Initialize barotropic variables: 481 IF( ll_init )THEN 482 sshbb_e(:,:) = 0._wp 483 ubb_e (:,:) = 0._wp 484 vbb_e (:,:) = 0._wp 485 sshb_e (:,:) = 0._wp 486 ub_e (:,:) = 0._wp 487 vb_e (:,:) = 0._wp 488 ENDIF 489 ! 466 490 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 467 491 sshn_e(:,:) = sshn (:,:) … … 533 557 END DO 534 558 END DO 535 CALL lbc_lnk ( zwx, 'U', 1._wp ) ; CALL lbc_lnk(zwy, 'V', 1._wp )559 CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp ) 536 560 ! 537 561 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points … … 611 635 END DO 612 636 END DO 613 CALL lbc_lnk ( zsshu_a, 'U', 1._wp ) ; CALL lbc_lnk(zsshv_a, 'V', 1._wp )637 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) 614 638 ENDIF 615 639 ! … … 685 709 END DO 686 710 ! 687 ELSEIF ( ln_dynvor_een ) THEN!== energy and enstrophy conserving scheme ==!711 ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN !== energy and enstrophy conserving scheme ==! 688 712 DO jj = 2, jpjm1 689 713 DO ji = fs_2, fs_jpim1 ! vector opt. … … 779 803 ! ! ----------------------- 780 804 ! 781 CALL lbc_lnk( ua_e , 'U', -1._wp ) ! local domain boundaries 782 CALL lbc_lnk( va_e , 'V', -1._wp ) 805 CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp ) 783 806 784 807 #if defined key_bdy … … 835 858 END DO 836 859 END DO 837 CALL lbc_lnk ( zsshu_a, 'U', 1._wp ) ; CALL lbc_lnk(zsshv_a, 'V', 1._wp ) ! Boundary conditions860 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 838 861 ENDIF 839 862 ! -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
r4624 r5832 15 15 !! 3.2 ! 2009-04 (R. Benshila) vvl: correction of een scheme 16 16 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 17 !! 3.7 ! 2014-04 (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity 17 18 !!---------------------------------------------------------------------- 18 19 … … 29 30 USE dommsk ! ocean mask 30 31 USE dynadv ! momentum advection (use ln_dynadv_vec value) 31 USE trd mod ! ocean dynamics trends32 USE trd mod_oce ! ocean variables trends32 USE trd_oce ! trends: ocean variables 33 USE trddyn ! trend manager: dynamics 33 34 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 34 35 USE prtctl ! Print control … … 50 51 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme 51 52 LOGICAL, PUBLIC :: ln_dynvor_een !: energy and enstrophy conserving scheme 53 LOGICAL, PUBLIC :: ln_dynvor_een_old !: energy and enstrophy conserving scheme (original formulation) 52 54 53 55 INTEGER :: nvor = 0 ! type of vorticity trend used … … 73 75 !! ** Action : - Update (ua,va) with the now vorticity term trend 74 76 !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative 75 !! and planetary vorticity trends) ('key_trddyn') 77 !! and planetary vorticity trends) and send them to trd_dyn 78 !! for futher diagnostics (l_trddyn=T) 76 79 !!---------------------------------------------------------------------- 77 80 INTEGER, INTENT( in ) :: kt ! ocean time-step index … … 108 111 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 109 112 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 110 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )113 CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 111 114 ztrdu(:,:,:) = ua(:,:,:) 112 115 ztrdv(:,:,:) = va(:,:,:) … … 114 117 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 115 118 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 116 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 117 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 119 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 118 120 ELSE 119 121 CALL vor_ene( kt, ntot, ua, va ) ! total vorticity … … 127 129 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 128 130 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 129 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )131 CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 130 132 ztrdu(:,:,:) = ua(:,:,:) 131 133 ztrdv(:,:,:) = va(:,:,:) … … 133 135 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 134 136 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 135 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 136 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 137 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 137 138 ELSE 138 139 CALL vor_ens( kt, ntot, ua, va ) ! total vorticity … … 146 147 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 147 148 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 148 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )149 CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 149 150 ztrdu(:,:,:) = ua(:,:,:) 150 151 ztrdv(:,:,:) = va(:,:,:) … … 152 153 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 153 154 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 154 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 155 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 155 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 156 156 ELSE 157 157 CALL vor_mix( kt ) ! total vorticity (mix=ens-ene) … … 165 165 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 166 166 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 167 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )167 CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 168 168 ztrdu(:,:,:) = ua(:,:,:) 169 169 ztrdv(:,:,:) = va(:,:,:) … … 171 171 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 172 172 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 173 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 174 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 173 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 175 174 ELSE 176 175 CALL vor_een( kt, ntot, ua, va ) ! total vorticity … … 211 210 !! 212 211 !! ** Action : - Update (ua,va) with the now vorticity term trend 213 !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative214 !! and planetary vorticity trends) ('key_trddyn')215 212 !! 216 213 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. … … 328 325 !! 329 326 !! ** Action : - Update (ua,va) arrays with the now vorticity term trend 330 !! - Save the trends in (ztrdu,ztrdv) in 2 parts (relative331 !! and planetary vorticity trends) ('key_trddyn')332 327 !! 333 328 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. … … 444 439 !! 445 440 !! ** Action : - Update (ua,va) arrays with the now vorticity term trend 446 !! - Save the trends in (ztrdu,ztrdv) in 2 parts (relative447 !! and planetary vorticity trends) ('key_trddyn')448 441 !! 449 442 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. … … 557 550 !! 558 551 !! ** Action : - Update (ua,va) with the now vorticity term trend 559 !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative560 !! and planetary vorticity trends) ('key_trddyn')561 552 !! 562 553 !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 … … 601 592 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' ) 602 593 ENDIF 603 ze3f(:,:,:) = 0. d0594 ze3f(:,:,:) = 0._wp 604 595 #endif 605 596 ENDIF 606 597 607 598 IF( kt == nit000 .OR. lk_vvl ) THEN ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points) 608 DO jk = 1, jpk 609 DO jj = 1, jpjm1 610 DO ji = 1, jpim1 611 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 612 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 613 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 614 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) 615 IF( ze3 /= 0._wp ) ze3f(ji,jj,jk) = zmsk / ze3 616 END DO 617 END DO 618 END DO 599 600 IF( ln_dynvor_een_old ) THEN ! original formulation 601 DO jk = 1, jpk 602 DO jj = 1, jpjm1 603 DO ji = 1, jpim1 604 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 605 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 606 IF( ze3 /= 0._wp ) ze3f(ji,jj,jk) = 4.0_wp / ze3 607 END DO 608 END DO 609 END DO 610 ELSE ! new formulation from NEMO 3.6 611 DO jk = 1, jpk 612 DO jj = 1, jpjm1 613 DO ji = 1, jpim1 614 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 615 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 616 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 617 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) 618 IF( ze3 /= 0._wp ) ze3f(ji,jj,jk) = zmsk / ze3 619 END DO 620 END DO 621 END DO 622 ENDIF 623 619 624 CALL lbc_lnk( ze3f, 'F', 1. ) 620 625 ENDIF … … 715 720 INTEGER :: ios ! Local integer output status for namelist read 716 721 !! 717 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een 722 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old 718 723 !!---------------------------------------------------------------------- 719 724 … … 736 741 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 737 742 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 743 WRITE(numout,*) ' enstrophy and energy conserving scheme (old) ln_dynvor_een_old= ', ln_dynvor_een_old 738 744 ENDIF 739 745 … … 759 765 IF( ln_dynvor_mix ) ioptio = ioptio + 1 760 766 IF( ln_dynvor_een ) ioptio = ioptio + 1 767 IF( ln_dynvor_een_old ) ioptio = ioptio + 1 761 768 IF( lk_esopa ) ioptio = 1 762 769 … … 767 774 IF( ln_dynvor_ens ) nvor = 1 768 775 IF( ln_dynvor_mix ) nvor = 2 769 IF( ln_dynvor_een ) nvor = 3776 IF( ln_dynvor_een .or. ln_dynvor_een_old ) nvor = 3 770 777 IF( lk_esopa ) nvor = -1 771 778 -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
r3294 r5832 16 16 USE dom_oce ! ocean space and time domain 17 17 USE sbc_oce ! surface boundary condition: ocean 18 USE trdmod_oce ! ocean variables trends 19 USE trdmod ! ocean dynamics trends 18 USE trd_oce ! trends: ocean variables 19 USE trddyn ! trend manager: dynamics 20 ! 20 21 USE in_out_manager ! I/O manager 21 USE lib_mpp 22 USE lib_mpp ! MPP library 22 23 USE prtctl ! Print control 23 USE wrk_nemo 24 USE timing 24 USE wrk_nemo ! Memory Allocation 25 USE timing ! Timing 25 26 26 27 IMPLICIT NONE 27 28 PRIVATE 28 29 29 PUBLIC dyn_zad ! routine called by step.F90 30 PUBLIC dyn_zad ! routine called by dynadv.F90 31 PUBLIC dyn_zad_zts ! routine called by dynadv.F90 30 32 31 33 !! * Substitutions … … 53 55 !! 54 56 !! ** Action : - Update (ua,va) with the vert. momentum adv. trends 55 !! - S ave the trends in (ztrdu,ztrdv) ('key_trddyn')57 !! - Send the trends to trddyn for diagnostics (l_trddyn=T) 56 58 !!---------------------------------------------------------------------- 57 59 INTEGER, INTENT(in) :: kt ! ocean time-step inedx … … 83 85 DO jj = 2, jpj ! vertical fluxes 84 86 DO ji = fs_2, jpi ! vector opt. 85 zww(ji,jj) = 0.25 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)87 zww(ji,jj) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) 86 88 END DO 87 89 END DO … … 93 95 END DO 94 96 END DO 95 DO jj = 2, jpjm1 ! Surface and bottom values set to zero 96 DO ji = fs_2, fs_jpim1 ! vector opt. 97 zwuw(ji,jj, 1 ) = 0.e0 98 zwvw(ji,jj, 1 ) = 0.e0 99 zwuw(ji,jj,jpk) = 0.e0 100 zwvw(ji,jj,jpk) = 0.e0 101 END DO 102 END DO 97 ! 98 ! Surface and bottom advective fluxes set to zero 99 IF ( ln_isfcav ) THEN 100 DO jj = 2, jpjm1 101 DO ji = fs_2, fs_jpim1 ! vector opt. 102 zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp 103 zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp 104 zwuw(ji,jj,jpk) = 0._wp 105 zwvw(ji,jj,jpk) = 0._wp 106 END DO 107 END DO 108 ELSE 109 DO jj = 2, jpjm1 110 DO ji = fs_2, fs_jpim1 ! vector opt. 111 zwuw(ji,jj, 1 ) = 0._wp 112 zwvw(ji,jj, 1 ) = 0._wp 113 zwuw(ji,jj,jpk) = 0._wp 114 zwvw(ji,jj,jpk) = 0._wp 115 END DO 116 END DO 117 END IF 103 118 104 119 DO jk = 1, jpkm1 ! Vertical momentum advection at u- and v-points … … 118 133 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 119 134 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 120 CALL trd_ mod(ztrdu, ztrdv, jpdyn_trd_zad, 'DYN', kt)135 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt ) 121 136 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 122 137 ENDIF … … 132 147 END SUBROUTINE dyn_zad 133 148 149 SUBROUTINE dyn_zad_zts ( kt ) 150 !!---------------------------------------------------------------------- 151 !! *** ROUTINE dynzad_zts *** 152 !! 153 !! ** Purpose : Compute the now vertical momentum advection trend and 154 !! add it to the general trend of momentum equation. This version 155 !! uses sub-timesteps for improved numerical stability with small 156 !! vertical grid sizes. This is especially relevant when using 157 !! embedded ice with thin surface boxes. 158 !! 159 !! ** Method : The now vertical advection of momentum is given by: 160 !! w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ] 161 !! w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ] 162 !! Add this trend to the general trend (ua,va): 163 !! (ua,va) = (ua,va) + w dz(u,v) 164 !! 165 !! ** Action : - Update (ua,va) with the vert. momentum adv. trends 166 !! - Save the trends in (ztrdu,ztrdv) ('key_trddyn') 167 !!---------------------------------------------------------------------- 168 INTEGER, INTENT(in) :: kt ! ocean time-step inedx 169 ! 170 INTEGER :: ji, jj, jk, jl ! dummy loop indices 171 INTEGER :: jnzts = 5 ! number of sub-timesteps for vertical advection 172 INTEGER :: jtb, jtn, jta ! sub timestep pointers for leap-frog/euler forward steps 173 REAL(wp) :: zua, zva ! temporary scalars 174 REAL(wp) :: zr_rdt ! temporary scalar 175 REAL(wp) :: z2dtzts ! length of Euler forward sub-timestep for vertical advection 176 REAL(wp) :: zts ! length of sub-timestep for vertical advection 177 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwuw , zwvw, zww 178 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 179 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zus , zvs 180 !!---------------------------------------------------------------------- 181 ! 182 IF( nn_timing == 1 ) CALL timing_start('dyn_zad_zts') 183 ! 184 CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 185 CALL wrk_alloc( jpi,jpj,jpk,3, zus, zvs ) 186 ! 187 IF( kt == nit000 ) THEN 188 IF(lwp)WRITE(numout,*) 189 IF(lwp)WRITE(numout,*) 'dyn_zad_zts : arakawa advection scheme with sub-timesteps' 190 ENDIF 191 192 IF( l_trddyn ) THEN ! Save ua and va trends 193 CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 194 ztrdu(:,:,:) = ua(:,:,:) 195 ztrdv(:,:,:) = va(:,:,:) 196 ENDIF 197 198 IF( neuler == 0 .AND. kt == nit000 ) THEN 199 z2dtzts = rdt / REAL( jnzts, wp ) ! = rdt (restart with Euler time stepping) 200 ELSE 201 z2dtzts = 2._wp * rdt / REAL( jnzts, wp ) ! = 2 rdt (leapfrog) 202 ENDIF 203 204 DO jk = 2, jpkm1 ! Calculate and store vertical fluxes 205 DO jj = 2, jpj 206 DO ji = fs_2, jpi ! vector opt. 207 zww(ji,jj,jk) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) 208 END DO 209 END DO 210 END DO 211 ! 212 ! Surface and bottom advective fluxes set to zero 213 DO jj = 2, jpjm1 214 DO ji = fs_2, fs_jpim1 ! vector opt. 215 zwuw(ji,jj, 1 ) = 0._wp 216 zwvw(ji,jj, 1 ) = 0._wp 217 zwuw(ji,jj,jpk) = 0._wp 218 zwvw(ji,jj,jpk) = 0._wp 219 END DO 220 END DO 221 222 ! Start with before values and use sub timestepping to reach after values 223 224 zus(:,:,:,1) = ub(:,:,:) 225 zvs(:,:,:,1) = vb(:,:,:) 226 227 DO jl = 1, jnzts ! Start of sub timestepping loop 228 229 IF( jl == 1 ) THEN ! Euler forward to kick things off 230 jtb = 1 ; jtn = 1 ; jta = 2 231 zts = z2dtzts 232 ELSEIF( jl == 2 ) THEN ! First leapfrog step 233 jtb = 1 ; jtn = 2 ; jta = 3 234 zts = 2._wp * z2dtzts 235 ELSE ! Shuffle pointers for subsequent leapfrog steps 236 jtb = MOD(jtb,3) + 1 237 jtn = MOD(jtn,3) + 1 238 jta = MOD(jta,3) + 1 239 ENDIF 240 241 DO jk = 2, jpkm1 ! Vertical momentum advection at level w and u- and v- vertical 242 DO jj = 2, jpjm1 ! vertical momentum advection at w-point 243 DO ji = fs_2, fs_jpim1 ! vector opt. 244 zwuw(ji,jj,jk) = ( zww(ji+1,jj ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) !* wumask(ji,jj,jk) 245 zwvw(ji,jj,jk) = ( zww(ji ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) !* wvmask(ji,jj,jk) 246 END DO 247 END DO 248 END DO 249 DO jk = 1, jpkm1 ! Vertical momentum advection at u- and v-points 250 DO jj = 2, jpjm1 251 DO ji = fs_2, fs_jpim1 ! vector opt. 252 ! ! vertical momentum advective trends 253 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 254 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 255 zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts 256 zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts 257 END DO 258 END DO 259 END DO 260 261 END DO ! End of sub timestepping loop 262 263 zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts ) 264 DO jk = 1, jpkm1 ! Recover trends over the outer timestep 265 DO jj = 2, jpjm1 266 DO ji = fs_2, fs_jpim1 ! vector opt. 267 ! ! vertical momentum advective trends 268 ! ! add the trends to the general momentum trends 269 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt 270 va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt 271 END DO 272 END DO 273 END DO 274 275 IF( l_trddyn ) THEN ! save the vertical advection trends for diagnostic 276 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 277 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 278 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt ) 279 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 280 ENDIF 281 ! ! Control print 282 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' zad - Ua: ', mask1=umask, & 283 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 284 ! 285 CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 286 CALL wrk_dealloc( jpi,jpj,jpk,3, zus, zvs ) 287 ! 288 IF( nn_timing == 1 ) CALL timing_stop('dyn_zad_zts') 289 ! 290 END SUBROUTINE dyn_zad_zts 291 134 292 !!====================================================================== 135 293 END MODULE dynzad -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
r3294 r5832 20 20 21 21 USE ldfdyn_oce ! ocean dynamics: lateral physics 22 USE trd mod ! ocean active dynamics and tracers trends23 USE trd mod_oce ! ocean variables trends22 USE trd_oce ! trends: ocean variables 23 USE trddyn ! trend manager: dynamics 24 24 USE in_out_manager ! I/O manager 25 25 USE lib_mpp ! MPP library … … 91 91 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 92 92 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 93 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_zdf, 'DYN', kt ) 94 ! 93 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 95 94 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 96 95 ENDIF -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r4370 r5832 70 70 REAL(wp) :: z1_p2dt, zcoef, zzwi, zzws, zrhs ! local scalars 71 71 REAL(wp) :: ze3ua, ze3va 72 !!----------------------------------------------------------------------73 74 72 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwd, zws 75 73 !!---------------------------------------------------------------------- … … 101 99 102 100 IF( ln_bfrimp ) THEN 103 # if defined key_vectopt_loop104 DO jj = 1, 1105 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)106 # else107 101 DO jj = 2, jpjm1 108 102 DO ji = 2, jpim1 109 # endif110 103 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 111 104 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) … … 114 107 END DO 115 108 END DO 109 IF ( ln_isfcav ) THEN 110 DO jj = 2, jpjm1 111 DO ji = 2, jpim1 112 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 113 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 114 IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu) 115 IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv) 116 END DO 117 END DO 118 END IF 116 119 ENDIF 117 120 … … 138 141 ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk) 139 142 va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk) 140 END DO141 ! Add bottom stress due to barotropic component only:143 END DO 144 ! Add bottom/top stress due to barotropic component only: 142 145 DO jj = 2, jpjm1 143 146 DO ji = fs_2, fs_jpim1 ! vector opt. … … 150 153 END DO 151 154 END DO 155 IF ( ln_isfcav ) THEN 156 DO jj = 2, jpjm1 157 DO ji = fs_2, fs_jpim1 ! vector opt. 158 ikbu = miku(ji,jj) ! top ocean level at u- and v-points 159 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 160 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl * fse3u_a(ji,jj,ikbu) 161 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl * fse3v_a(ji,jj,ikbv) 162 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua 163 va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va 164 END DO 165 END DO 166 END IF 152 167 ENDIF 153 168 #endif … … 164 179 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl * fse3u_a(ji,jj,jk) ! after scale factor at T-point 165 180 zcoef = - p2dt / ze3ua 166 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk )167 zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk)168 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)169 zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1)170 zwd(ji,jj,jk) = 1._wp - z wi(ji,jj,jk)- zzws181 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk ) 182 zwi(ji,jj,jk) = zzwi * wumask(ji,jj,jk ) 183 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 184 zws(ji,jj,jk) = zzws * wumask(ji,jj,jk+1) 185 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 171 186 END DO 172 187 END DO … … 194 209 !----------------------------------------------------------------------- 195 210 ! 196 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 197 DO jj = 2, jpjm1 198 DO ji = fs_2, fs_jpim1 ! vector opt. 199 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 200 END DO 201 END DO 202 END DO 203 ! 204 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 205 DO ji = fs_2, fs_jpim1 ! vector opt. 206 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl * fse3u_a(ji,jj,1) 207 #if defined key_dynspg_ts 208 ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 209 & / ( ze3ua * rau0 ) 210 #else 211 ua(ji,jj,1) = ub(ji,jj,1) + p2dt *(ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 212 & / ( fse3u(ji,jj,1) * rau0 ) ) 213 #endif 214 END DO 215 END DO 211 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 216 212 DO jk = 2, jpkm1 217 213 DO jj = 2, jpjm1 218 214 DO ji = fs_2, fs_jpim1 ! vector opt. 215 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 216 END DO 217 END DO 218 END DO 219 ! 220 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 221 DO ji = fs_2, fs_jpim1 ! vector opt. 222 #if defined key_dynspg_ts 223 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl * fse3u_a(ji,jj,1) 224 ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 225 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 226 #else 227 ua(ji,jj,1) = ub(ji,jj,1) & 228 & + p2dt *(ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 229 & / ( fse3u(ji,jj,1) * rau0 ) * umask(ji,jj,1) ) 230 #endif 231 END DO 232 END DO 233 DO jk = 2, jpkm1 234 DO jj = 2, jpjm1 235 DO ji = fs_2, fs_jpim1 219 236 #if defined key_dynspg_ts 220 237 zrhs = ua(ji,jj,jk) ! zrhs=right hand side … … 233 250 END DO 234 251 DO jk = jpk-2, 1, -1 235 DO jj = 2, jpjm1 236 DO ji = fs_2, fs_jpim1 ! vector opt.252 DO jj = 2, jpjm1 253 DO ji = fs_2, fs_jpim1 237 254 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 238 255 END DO … … 263 280 zcoef = - p2dt / ze3va 264 281 zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk ) 265 zwi(ji,jj,jk) = zzwi * vmask(ji,jj,jk)282 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk) 266 283 zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 267 zws(ji,jj,jk) = zzws * vmask(ji,jj,jk+1)268 zwd(ji,jj,jk) = 1._wp - z wi(ji,jj,jk)- zzws284 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 285 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 269 286 END DO 270 287 END DO … … 292 309 !----------------------------------------------------------------------- 293 310 ! 294 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 311 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 312 DO jk = 2, jpkm1 295 313 DO jj = 2, jpjm1 296 314 DO ji = fs_2, fs_jpim1 ! vector opt. … … 302 320 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 303 321 DO ji = fs_2, fs_jpim1 ! vector opt. 322 #if defined key_dynspg_ts 304 323 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl * fse3v_a(ji,jj,1) 305 #if defined key_dynspg_ts306 324 va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 307 325 & / ( ze3va * rau0 ) 308 326 #else 309 va(ji,jj,1) = vb(ji,jj,1) + p2dt *(va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 327 va(ji,jj,1) = vb(ji,jj,1) & 328 & + p2dt *(va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 310 329 & / ( fse3v(ji,jj,1) * rau0 ) ) 311 330 #endif … … 331 350 END DO 332 351 DO jk = jpk-2, 1, -1 333 DO jj = 2, jpjm1 334 DO ji = fs_2, fs_jpim1 ! vector opt.352 DO jj = 2, jpjm1 353 DO ji = fs_2, fs_jpim1 335 354 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 336 355 END DO … … 352 371 !! restore bottom layer avmu(v) 353 372 IF( ln_bfrimp ) THEN 354 # if defined key_vectopt_loop 355 DO jj = 1, 1 356 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) 357 # else 358 DO jj = 2, jpjm1 359 DO ji = 2, jpim1 360 # endif 361 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 362 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 363 avmu(ji,jj,ikbu+1) = 0.e0 364 avmv(ji,jj,ikbv+1) = 0.e0 365 END DO 366 END DO 373 DO jj = 2, jpjm1 374 DO ji = 2, jpim1 375 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 376 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 377 avmu(ji,jj,ikbu+1) = 0.e0 378 avmv(ji,jj,ikbv+1) = 0.e0 379 END DO 380 END DO 381 IF (ln_isfcav) THEN 382 DO jj = 2, jpjm1 383 DO ji = 2, jpim1 384 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 385 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 386 IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0 387 IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0 388 END DO 389 END DO 390 END IF 367 391 ENDIF 368 392 ! -
branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r4486 r5832 21 21 USE domvvl ! Variable volume 22 22 USE divcur ! hor. divergence and curl (div & cur routines) 23 USE iom ! I/O library24 23 USE restart ! only for lrst_oce 25 24 USE in_out_manager ! I/O manager … … 31 30 USE bdy_par 32 31 USE bdydyn2d ! bdy_ssh routine 33 USE diaar5, ONLY: lk_diaar534 USE iom35 32 #if defined key_agrif 36 33 USE agrif_opa_update … … 111 108 ! 112 109 z1_rau0 = 0.5_wp * r1_rau0 113 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1)110 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 114 111 115 112 #if ! defined key_dynspg_ts … … 138 135 ! ! outputs ! 139 136 ! !------------------------------! 140 CALL iom_put( "ssh" , sshn ) ! sea surface height141 CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height142 137 ! 143 138 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 ) … … 229 224 #endif 230 225 ! 231 ! !------------------------------!232 ! ! outputs !233 ! !------------------------------!234 CALL iom_put( "woce", wn ) ! vertical velocity235 IF( lk_diaar5 ) THEN ! vertical mass transport & its square value236 CALL wrk_alloc( jpi, jpj, z2d )237 CALL wrk_alloc( jpi, jpj, jpk, z3d )238 ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.239 z2d(:,:) = rau0 * e12t(:,:)240 DO jk = 1, jpk241 z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)242 END DO243 CALL iom_put( "w_masstr" , z3d )244 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )245 CALL wrk_dealloc( jpi, jpj, z2d )246 CALL wrk_dealloc( jpi, jpj, jpk, z3d )247 ENDIF248 !249 226 IF( nn_timing == 1 ) CALL timing_stop('wzv') 250 227 … … 291 268 ELSE !** Leap-Frog time-stepping: Asselin filter + swap 292 269 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) ! before <-- now filtered 293 IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)270 IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) - rnf_b(:,:) + rnf(:,:) ) * ssmask(:,:) 294 271 sshn(:,:) = ssha(:,:) ! now <-- after 295 272 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.