Changeset 10928
- Timestamp:
- 2019-05-03T17:44:56+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src
- Files:
-
- 16 edited
-
OCE/DOM/istate.F90 (modified) (7 diffs)
-
OCE/DYN/divhor.F90 (modified) (3 diffs)
-
OCE/DYN/dynadv_cen2.F90 (modified) (1 diff)
-
OCE/DYN/dynadv_ubs.F90 (modified) (1 diff)
-
OCE/DYN/dynhpg.F90 (modified) (64 diffs)
-
OCE/DYN/dynkeg.F90 (modified) (1 diff)
-
OCE/DYN/dynldf.F90 (modified) (3 diffs)
-
OCE/DYN/dynldf_iso.F90 (modified) (22 diffs)
-
OCE/DYN/dynldf_lap_blp.F90 (modified) (7 diffs)
-
OCE/DYN/dynspg.F90 (modified) (1 diff)
-
OCE/DYN/dynvor.F90 (modified) (1 diff)
-
OCE/DYN/dynzad.F90 (modified) (1 diff)
-
OCE/DYN/dynzdf.F90 (modified) (4 diffs)
-
OCE/nemogcm.F90 (modified) (3 diffs)
-
OCE/step.F90 (modified) (1 diff)
-
OFF/nemogcm.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/istate.F90
r10922 r10928 51 51 CONTAINS 52 52 53 SUBROUTINE istate_init( Kbb, Kmm )53 SUBROUTINE istate_init( Kbb, Kmm, Kaa ) 54 54 !!---------------------------------------------------------------------- 55 55 !! *** ROUTINE istate_init *** … … 57 57 !! ** Purpose : Initialization of the dynamics and tracer fields. 58 58 !!---------------------------------------------------------------------- 59 INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices 59 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! ocean time level indices 60 ! 60 61 INTEGER :: ji, jj, jk ! dummy loop indices 61 62 !!gm see comment further down … … 77 78 rhd (:,:,: ) = 0._wp ; rhop (:,:,: ) = 0._wp ! set one for all to 0 at level jpk 78 79 rn2b (:,:,: ) = 0._wp ; rn2 (:,:,: ) = 0._wp ! set one for all to 0 at levels 1 and jpk 79 ts a (:,:,:,:) = 0._wp ! set one for all to 0 at level jpk80 ts (:,:,:,:,Kaa) = 0._wp ! set one for all to 0 at level jpk 80 81 rab_b(:,:,:,:) = 0._wp ; rab_n(:,:,:,:) = 0._wp ! set one for all to 0 at level jpk 81 82 #if defined key_agrif 82 u a (:,:,:) = 0._wp ! used in agrif_oce_sponge at initialization83 v a (:,:,:) = 0._wp ! used in agrif_oce_sponge at initialization83 uu (:,:,: ,Kaa) = 0._wp ! used in agrif_oce_sponge at initialization 84 vv (:,:,: ,Kaa) = 0._wp ! used in agrif_oce_sponge at initialization 84 85 #endif 85 86 … … 98 99 ! 99 100 IF( ln_tsd_init ) THEN 100 CALL dta_tsd( nit000, ts b) ! read 3D T and S data at nit000101 CALL dta_tsd( nit000, ts(:,:,:,:,Kbb) ) ! read 3D T and S data at nit000 101 102 ! 102 ssh b(:,:) = 0._wp ! set the ocean at rest103 ssh(:,:,Kbb) = 0._wp ! set the ocean at rest 103 104 IF( ll_wd ) THEN 104 ssh b(:,:) = -ssh_ref ! Added in 30 here for bathy that adds 30 as Iterative test CEOD105 ssh(:,:,Kbb) = -ssh_ref ! Added in 30 here for bathy that adds 30 as Iterative test CEOD 105 106 ! 106 107 ! Apply minimum wetdepth criterion … … 108 109 DO jj = 1,jpj 109 110 DO ji = 1,jpi 110 IF( ht_0(ji,jj) + ssh b(ji,jj) < rn_wdmin1 ) THEN111 ssh b(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )111 IF( ht_0(ji,jj) + ssh(ji,jj,Kbb) < rn_wdmin1 ) THEN 112 ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) ) 112 113 ENDIF 113 114 END DO 114 115 END DO 115 116 ENDIF 116 u b (:,:,:) = 0._wp117 v b (:,:,:) = 0._wp117 uu (:,:,:,Kbb) = 0._wp 118 vv (:,:,:,Kbb) = 0._wp 118 119 ! 119 120 ELSE ! user defined initial T and S 120 CALL usr_def_istate( gdept _b, tmask, tsb, ub, vb, sshb)121 CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 121 122 ENDIF 122 ts n (:,:,:,:) = tsb (:,:,:,:) ! set now values from to before ones123 ssh n (:,:) = sshb(:,:)124 u n (:,:,:) = ub (:,:,:)125 v n (:,:,:) = vb (:,:,:)126 hdiv n(:,:,jpk) = 0._wp ! bottom divergence set one for 0 to zero at jpk level127 CALL div_hor( 0, Kmm ) ! compute interior hdivnvalue128 !!gm hdiv n(:,:,:) = 0._wp123 ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones 124 ssh (:,:,Kmm) = ssh(:,:,Kbb) 125 uu (:,:,:,Kmm) = uu (:,:,:,Kbb) 126 vv (:,:,:,Kmm) = vv (:,:,:,Kbb) 127 hdiv(:,:,jpk) = 0._wp ! bottom divergence set one for 0 to zero at jpk level 128 CALL div_hor( 0, Kmm ) ! compute interior hdiv value 129 !!gm hdiv(:,:,:) = 0._wp 129 130 130 131 !!gm POTENTIAL BUG : 131 !!gm ISSUE : if ssh b/= 0 then, in non linear free surface, the e3._n, e3._b should be recomputed132 !!gm ISSUE : if ssh(:,:,Kbb) /= 0 then, in non linear free surface, the e3._n, e3._b should be recomputed 132 133 !! as well as gdept and gdepw.... !!!!! 133 134 !! ===>>>> probably a call to domvvl initialisation here.... … … 139 140 ! ALLOCATE( zuvd(jpi,jpj,jpk,2) ) 140 141 ! CALL dta_uvd( nit000, zuvd ) 141 ! u b(:,:,:) = zuvd(:,:,:,1) ; un(:,:,:) = ub(:,:,:)142 ! v b(:,:,:) = zuvd(:,:,:,2) ; vn(:,:,:) = vb(:,:,:)142 ! uu(:,:,:,Kbb) = zuvd(:,:,:,1) ; uu(:,:,:,Kmm) = uu(:,:,:,Kbb) 143 ! vv(:,:,:,Kbb) = zuvd(:,:,:,2) ; vv(:,:,:,Kmm) = vv(:,:,:,Kbb) 143 144 ! DEALLOCATE( zuvd ) 144 145 ! ENDIF 145 146 ! 146 147 !!gm This is to be changed !!!! 147 ! ! - ML - ssh n could be modified by istate_eel, so that initialization of e3t_bis done here148 ! ! - ML - ssh(:,:,Kmm) could be modified by istate_eel, so that initialization of e3t(:,:,:,Kbb) is done here 148 149 ! IF( .NOT.ln_linssh ) THEN 149 150 ! DO jk = 1, jpk 150 ! e3t _b(:,:,jk) = e3t_n(:,:,jk)151 ! e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm) 151 152 ! END DO 152 153 ! ENDIF … … 158 159 ! Do it whatever the free surface method, these arrays being eventually used 159 160 ! 160 u n_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp161 u b_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp161 uu_b(:,:,Kmm) = 0._wp ; vv_b(:,:,Kmm) = 0._wp 162 uu_b(:,:,Kbb) = 0._wp ; vv_b(:,:,Kbb) = 0._wp 162 163 ! 163 !!gm the use of umsak & vmask is not necessary below as u n, vn, ub, vbare always masked164 !!gm the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked 164 165 DO jk = 1, jpkm1 165 166 DO jj = 1, jpj 166 167 DO ji = 1, jpi 167 u n_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)168 v n_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)168 uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk) 169 vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk) 169 170 ! 170 u b_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)171 v b_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)171 uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk) 172 vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk) 172 173 END DO 173 174 END DO 174 175 END DO 175 176 ! 176 u n_b(:,:) = un_b(:,:) * r1_hu_n(:,:)177 v n_b(:,:) = vn_b(:,:) * r1_hv_n(:,:)177 uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu_n(:,:) 178 vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv_n(:,:) 178 179 ! 179 u b_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)180 v b_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)180 uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu_b(:,:) 181 vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv_b(:,:) 181 182 ! 182 183 END SUBROUTINE istate_init -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/divhor.F90
r10922 r10928 55 55 !! 56 56 !! ** Method : the now divergence is computed as : 57 !! hdiv n= 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )57 !! hdiv = 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 58 58 !! and correct with runoff inflow (div_rnf) and cross land flow (div_cla) 59 59 !! 60 !! ** Action : - update hdiv n, the now horizontal divergence60 !! ** Action : - update hdiv, the now horizontal divergence 61 61 !!---------------------------------------------------------------------- 62 62 INTEGER, INTENT(in) :: kt ! ocean time-step index 63 INTEGER, INTENT(in) :: Kmm ! ocean time -level index63 INTEGER, INTENT(in) :: Kmm ! ocean time level index 64 64 ! 65 65 INTEGER :: ji, jj, jk ! dummy loop indices … … 78 78 DO jj = 2, jpjm1 79 79 DO ji = fs_2, fs_jpim1 ! vector opt. 80 hdiv n(ji,jj,jk) = ( e2u(ji ,jj) * e3u_n(ji ,jj,jk) * un(ji ,jj,jk) &81 & - e2u(ji-1,jj) * e3u _n(ji-1,jj,jk) * un(ji-1,jj,jk) &82 & + e1v(ji,jj ) * e3v _n(ji,jj ,jk) * vn(ji,jj ,jk) &83 & - e1v(ji,jj-1) * e3v _n(ji,jj-1,jk) * vn(ji,jj-1,jk) ) &84 & * r1_e1e2t(ji,jj) / e3t _n(ji,jj,jk)80 hdiv(ji,jj,jk) = ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * uu(ji ,jj,jk,Kmm) & 81 & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm) & 82 & + e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * vv(ji,jj ,jk,Kmm) & 83 & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm) ) & 84 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 85 85 END DO 86 86 END DO … … 88 88 #if defined key_agrif 89 89 IF( .NOT. Agrif_Root() ) THEN 90 IF( nbondi == -1 .OR. nbondi == 2 ) hdiv n( 2 , : ,:) = 0._wp ! west91 IF( nbondi == 1 .OR. nbondi == 2 ) hdiv n( nlci-1, : ,:) = 0._wp ! east92 IF( nbondj == -1 .OR. nbondj == 2 ) hdiv n( : , 2 ,:) = 0._wp ! south93 IF( nbondj == 1 .OR. nbondj == 2 ) hdiv n( : ,nlcj-1,:) = 0._wp ! north90 IF( nbondi == -1 .OR. nbondi == 2 ) hdiv( 2 , : ,:) = 0._wp ! west 91 IF( nbondi == 1 .OR. nbondi == 2 ) hdiv( nlci-1, : ,:) = 0._wp ! east 92 IF( nbondj == -1 .OR. nbondj == 2 ) hdiv( : , 2 ,:) = 0._wp ! south 93 IF( nbondj == 1 .OR. nbondj == 2 ) hdiv( : ,nlcj-1,:) = 0._wp ! north 94 94 ENDIF 95 95 #endif 96 96 ! 97 IF( ln_rnf ) CALL sbc_rnf_div( hdiv n, Kmm ) !== runoffs ==! (update hdivnfield)97 IF( ln_rnf ) CALL sbc_rnf_div( hdiv, Kmm ) !== runoffs ==! (update hdiv field) 98 98 ! 99 99 #if defined key_asminc 100 IF( ln_sshinc .AND. ln_asmiau ) CALL ssh_asm_div( kt, hdiv n ) !== SSH assimilation ==! (update hdivnfield)100 IF( ln_sshinc .AND. ln_asmiau ) CALL ssh_asm_div( kt, hdiv ) !== SSH assimilation ==! (update hdiv field) 101 101 ! 102 102 #endif 103 IF( ln_isf ) CALL sbc_isf_div( hdiv n, Kmm ) !== ice shelf ==! (update hdivnfield)103 IF( ln_isf ) CALL sbc_isf_div( hdiv, Kmm ) !== ice shelf ==! (update hdiv field) 104 104 ! 105 IF( ln_iscpl .AND. ln_hsb ) CALL iscpl_div( hdiv n ) !== ice sheet ==! (update hdivnfield)105 IF( ln_iscpl .AND. ln_hsb ) CALL iscpl_div( hdiv ) !== ice sheet ==! (update hdiv field) 106 106 ! 107 CALL lbc_lnk( 'divhor', hdiv n, 'T', 1. ) ! (no sign change)107 CALL lbc_lnk( 'divhor', hdiv, 'T', 1. ) ! (no sign change) 108 108 ! 109 109 IF( ln_timing ) CALL timing_stop('div_hor') -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_cen2.F90
r10877 r10928 141 141 ENDIF 142 142 ! ! Control print 143 IF(ln_ctl) CALL prt_ctl( tab3d_1= ua, clinfo1=' cen2 adv - Ua: ', mask1=umask, &144 & tab3d_2= va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )143 IF(ln_ctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' cen2 adv - Ua: ', mask1=umask, & 144 & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 145 145 ! 146 146 END SUBROUTINE dyn_adv_cen2 -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_ubs.F90
r10877 r10928 234 234 ENDIF 235 235 ! ! Control print 236 IF(ln_ctl) CALL prt_ctl( tab3d_1= ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask, &237 & tab3d_2= va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )236 IF(ln_ctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask, & 237 & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 238 238 ! 239 239 END SUBROUTINE dyn_adv_ubs -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynhpg.F90
r10491 r10928 81 81 CONTAINS 82 82 83 SUBROUTINE dyn_hpg( kt )83 SUBROUTINE dyn_hpg( kt, Kmm, puu, pvv, Krhs ) 84 84 !!--------------------------------------------------------------------- 85 85 !! *** ROUTINE dyn_hpg *** … … 88 88 !! using the scheme defined in the namelist 89 89 !! 90 !! ** Action : - Update ( ua,va) with the now hydrastatic pressure trend90 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 91 91 !! - send trends to trd_dyn for futher diagnostics (l_trddyn=T) 92 92 !!---------------------------------------------------------------------- 93 INTEGER, INTENT(in) :: kt ! ocean time-step index 93 INTEGER , INTENT( in ) :: kt ! ocean time-step index 94 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 95 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 96 ! 94 97 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 95 98 !!---------------------------------------------------------------------- … … 97 100 IF( ln_timing ) CALL timing_start('dyn_hpg') 98 101 ! 99 IF( l_trddyn ) THEN ! Temporary saving of ua and vatrends (l_trddyn)102 IF( l_trddyn ) THEN ! Temporary saving of puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends (l_trddyn) 100 103 ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 101 ztrdu(:,:,:) = ua(:,:,:)102 ztrdv(:,:,:) = va(:,:,:)104 ztrdu(:,:,:) = puu(:,:,:,Krhs) 105 ztrdv(:,:,:) = pvv(:,:,:,Krhs) 103 106 ENDIF 104 107 ! 105 108 SELECT CASE ( nhpg ) ! Hydrostatic pressure gradient computation 106 CASE ( np_zco ) ; CALL hpg_zco ( kt )! z-coordinate107 CASE ( np_zps ) ; CALL hpg_zps ( kt )! z-coordinate plus partial steps (interpolation)108 CASE ( np_sco ) ; CALL hpg_sco ( kt )! s-coordinate (standard jacobian formulation)109 CASE ( np_djc ) ; CALL hpg_djc ( kt )! s-coordinate (Density Jacobian with Cubic polynomial)110 CASE ( np_prj ) ; CALL hpg_prj ( kt )! s-coordinate (Pressure Jacobian scheme)111 CASE ( np_isf ) ; CALL hpg_isf ( kt )! s-coordinate similar to sco modify for ice shelf109 CASE ( np_zco ) ; CALL hpg_zco ( kt, Kmm, puu, pvv, Krhs ) ! z-coordinate 110 CASE ( np_zps ) ; CALL hpg_zps ( kt, Kmm, puu, pvv, Krhs ) ! z-coordinate plus partial steps (interpolation) 111 CASE ( np_sco ) ; CALL hpg_sco ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (standard jacobian formulation) 112 CASE ( np_djc ) ; CALL hpg_djc ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (Density Jacobian with Cubic polynomial) 113 CASE ( np_prj ) ; CALL hpg_prj ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (Pressure Jacobian scheme) 114 CASE ( np_isf ) ; CALL hpg_isf ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate similar to sco modify for ice shelf 112 115 END SELECT 113 116 ! 114 117 IF( l_trddyn ) THEN ! save the hydrostatic pressure gradient trends for momentum trend diagnostics 115 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)116 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)118 ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 119 ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) 117 120 CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt ) 118 121 DEALLOCATE( ztrdu , ztrdv ) 119 122 ENDIF 120 123 ! 121 IF(ln_ctl) CALL prt_ctl( tab3d_1= ua, clinfo1=' hpg - Ua: ', mask1=umask, &122 & tab3d_2= va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )124 IF(ln_ctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' hpg - Ua: ', mask1=umask, & 125 & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 123 126 ! 124 127 IF( ln_timing ) CALL timing_stop('dyn_hpg') … … 127 130 128 131 129 SUBROUTINE dyn_hpg_init 132 SUBROUTINE dyn_hpg_init( Kmm ) 130 133 !!---------------------------------------------------------------------- 131 134 !! *** ROUTINE dyn_hpg_init *** … … 137 140 !! with the type of vertical coordinate used (zco, zps, sco) 138 141 !!---------------------------------------------------------------------- 142 INTEGER, INTENT( in ) :: Kmm ! ocean time level index 143 ! 139 144 INTEGER :: ioptio = 0 ! temporary integer 140 145 INTEGER :: ios ! Local integer output status for namelist read … … 228 233 ! 229 234 DO jk = 1, jpk !- compute density of the water displaced by the ice shelf 230 CALL eos( zts_top(:,:,:), gdept _n(:,:,jk), zrhd(:,:,jk) )235 CALL eos( zts_top(:,:,:), gdept(:,:,jk,Kmm), zrhd(:,:,jk) ) 231 236 END DO 232 237 ! … … 239 244 DO ji = 1, jpi ! divided by 2 later 240 245 ikt = mikt(ji,jj) 241 ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w _n(ji,jj,1) * (1._wp - tmask(ji,jj,1))246 ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w(ji,jj,1,Kmm) * (1._wp - tmask(ji,jj,1)) 242 247 DO jk = 2, ikt-1 243 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w _n(ji,jj,jk) &248 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w(ji,jj,jk,Kmm) & 244 249 & * (1._wp - tmask(ji,jj,jk)) 245 250 END DO 246 251 IF (ikt >= 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 247 & * ( risfdep(ji,jj) - gdept _n(ji,jj,ikt-1) )252 & * ( risfdep(ji,jj) - gdept(ji,jj,ikt-1,Kmm) ) 248 253 END DO 249 254 END DO … … 256 261 257 262 258 SUBROUTINE hpg_zco( kt )263 SUBROUTINE hpg_zco( kt, Kmm, puu, pvv, Krhs ) 259 264 !!--------------------------------------------------------------------- 260 265 !! *** ROUTINE hpg_zco *** … … 266 271 !! level: zhpi = grav ..... 267 272 !! zhpj = grav ..... 268 !! add it to the general momentum trend (ua,va). 269 !! ua = ua - 1/e1u * zhpi 270 !! va = va - 1/e2v * zhpj 271 !! 272 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 273 !!---------------------------------------------------------------------- 274 INTEGER, INTENT(in) :: kt ! ocean time-step index 273 !! add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 274 !! puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 275 !! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 276 !! 277 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 278 !!---------------------------------------------------------------------- 279 INTEGER , INTENT( in ) :: kt ! ocean time-step index 280 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 281 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 275 282 ! 276 283 INTEGER :: ji, jj, jk ! dummy loop indices … … 290 297 DO jj = 2, jpjm1 291 298 DO ji = fs_2, fs_jpim1 ! vector opt. 292 zcoef1 = zcoef0 * e3w _n(ji,jj,1)299 zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) 293 300 ! hydrostatic pressure gradient 294 301 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 295 302 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 296 303 ! add to the general momentum trend 297 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)298 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)304 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 305 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 299 306 END DO 300 307 END DO … … 305 312 DO jj = 2, jpjm1 306 313 DO ji = fs_2, fs_jpim1 ! vector opt. 307 zcoef1 = zcoef0 * e3w _n(ji,jj,jk)314 zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm) 308 315 ! hydrostatic pressure gradient 309 316 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & … … 315 322 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 316 323 ! add to the general momentum trend 317 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)318 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)324 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 325 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 319 326 END DO 320 327 END DO … … 324 331 325 332 326 SUBROUTINE hpg_zps( kt )333 SUBROUTINE hpg_zps( kt, Kmm, puu, pvv, Krhs ) 327 334 !!--------------------------------------------------------------------- 328 335 !! *** ROUTINE hpg_zps *** … … 330 337 !! ** Method : z-coordinate plus partial steps case. blahblah... 331 338 !! 332 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 333 !!---------------------------------------------------------------------- 334 INTEGER, INTENT(in) :: kt ! ocean time-step index 339 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 340 !!---------------------------------------------------------------------- 341 INTEGER , INTENT( in ) :: kt ! ocean time-step index 342 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 343 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 335 344 !! 336 345 INTEGER :: ji, jj, jk ! dummy loop indices … … 347 356 348 357 ! Partial steps: bottom before horizontal gradient of t, s, rd at the last ocean level 349 !jc CALL zps_hde ( kt, jpts, ts n, gtsu, gtsv, rhd, gru , grv )358 !jc CALL zps_hde ( kt, jpts, ts(:,:,:,:,Kmm), gtsu, gtsv, rhd, gru , grv ) 350 359 351 360 ! Local constant initialization … … 355 364 DO jj = 2, jpjm1 356 365 DO ji = fs_2, fs_jpim1 ! vector opt. 357 zcoef1 = zcoef0 * e3w _n(ji,jj,1)366 zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm) 358 367 ! hydrostatic pressure gradient 359 368 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 360 369 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 361 370 ! add to the general momentum trend 362 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)363 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)371 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 372 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 364 373 END DO 365 374 END DO … … 369 378 DO jj = 2, jpjm1 370 379 DO ji = fs_2, fs_jpim1 ! vector opt. 371 zcoef1 = zcoef0 * e3w _n(ji,jj,jk)380 zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm) 372 381 ! hydrostatic pressure gradient 373 382 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & … … 379 388 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 380 389 ! add to the general momentum trend 381 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)382 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)390 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 391 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 383 392 END DO 384 393 END DO … … 390 399 iku = mbku(ji,jj) 391 400 ikv = mbkv(ji,jj) 392 zcoef2 = zcoef0 * MIN( e3w _n(ji,jj,iku), e3w_n(ji+1,jj ,iku) )393 zcoef3 = zcoef0 * MIN( e3w _n(ji,jj,ikv), e3w_n(ji ,jj+1,ikv) )401 zcoef2 = zcoef0 * MIN( e3w(ji,jj,iku,Kmm), e3w(ji+1,jj ,iku,Kmm) ) 402 zcoef3 = zcoef0 * MIN( e3w(ji,jj,ikv,Kmm), e3w(ji ,jj+1,ikv,Kmm) ) 394 403 IF( iku > 1 ) THEN ! on i-direction (level 2 or more) 395 ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value404 puu (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) - zhpi(ji,jj,iku) ! subtract old value 396 405 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one 397 406 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj) 398 ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend407 puu (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend 399 408 ENDIF 400 409 IF( ikv > 1 ) THEN ! on j-direction (level 2 or more) 401 va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value410 pvv (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) - zhpj(ji,jj,ikv) ! subtract old value 402 411 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one 403 412 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj) 404 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend413 pvv (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 405 414 ENDIF 406 415 END DO … … 410 419 411 420 412 SUBROUTINE hpg_sco( kt )421 SUBROUTINE hpg_sco( kt, Kmm, puu, pvv, Krhs ) 413 422 !!--------------------------------------------------------------------- 414 423 !! *** ROUTINE hpg_sco *** … … 422 431 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] 423 432 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ] 424 !! add it to the general momentum trend (ua,va). 425 !! ua = ua - 1/e1u * zhpi 426 !! va = va - 1/e2v * zhpj 427 !! 428 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 429 !!---------------------------------------------------------------------- 430 INTEGER, INTENT(in) :: kt ! ocean time-step index 433 !! add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 434 !! puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 435 !! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 436 !! 437 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 438 !!---------------------------------------------------------------------- 439 INTEGER , INTENT( in ) :: kt ! ocean time-step index 440 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 441 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 431 442 !! 432 443 INTEGER :: ji, jj, jk, jii, jjj ! dummy loop indices … … 453 464 DO jj = 2, jpjm1 454 465 DO ji = 2, jpim1 455 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji+1,jj) ) > &466 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 456 467 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 457 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) &468 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) & 458 469 & > rn_wdmin1 + rn_wdmin2 459 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( &460 & MAX( ssh n(ji,jj) , sshn(ji+1,jj) ) > &470 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. ( & 471 & MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 461 472 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 462 473 … … 464 475 zcpx(ji,jj) = 1.0_wp 465 476 ELSE IF(ll_tmp2) THEN 466 ! no worries about ssh n(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here467 zcpx(ji,jj) = ABS( (ssh n(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &468 & / (ssh n(ji+1,jj) - sshn(ji ,jj)) )477 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 478 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 479 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 469 480 ELSE 470 481 zcpx(ji,jj) = 0._wp 471 482 END IF 472 483 473 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji,jj+1) ) > &484 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 474 485 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 475 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) &486 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) & 476 487 & > rn_wdmin1 + rn_wdmin2 477 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( &478 & MAX( ssh n(ji,jj) , sshn(ji,jj+1) ) > &488 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. ( & 489 & MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 479 490 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 480 491 … … 482 493 zcpy(ji,jj) = 1.0_wp 483 494 ELSE IF(ll_tmp2) THEN 484 ! no worries about ssh n(ji,jj+1) - sshn(ji,jj) = 0, it won't happen ! here485 zcpy(ji,jj) = ABS( (ssh n(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &486 & / (ssh n(ji,jj+1) - sshn(ji,jj)) )495 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 496 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 497 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 487 498 ELSE 488 499 zcpy(ji,jj) = 0._wp … … 497 508 DO ji = fs_2, fs_jpim1 ! vector opt. 498 509 ! hydrostatic pressure gradient along s-surfaces 499 zhpi(ji,jj,1) = zcoef0 * ( e3w _n(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) &500 & - e3w _n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e1u(ji,jj)501 zhpj(ji,jj,1) = zcoef0 * ( e3w _n(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) &502 & - e3w _n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e2v(ji,jj)510 zhpi(ji,jj,1) = zcoef0 * ( e3w(ji+1,jj ,1,Kmm) * ( znad + rhd(ji+1,jj ,1) ) & 511 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e1u(ji,jj) 512 zhpj(ji,jj,1) = zcoef0 * ( e3w(ji ,jj+1,1,Kmm) * ( znad + rhd(ji ,jj+1,1) ) & 513 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e2v(ji,jj) 503 514 ! s-coordinate pressure gradient correction 504 515 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 505 & * ( gde3w _n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)516 & * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 506 517 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 507 & * ( gde3w _n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)518 & * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 508 519 ! 509 520 IF( ln_wd_il ) THEN … … 515 526 ! 516 527 ! add to the general momentum trend 517 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap518 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap528 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) + zuap 529 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) + zvap 519 530 END DO 520 531 END DO … … 526 537 ! hydrostatic pressure gradient along s-surfaces 527 538 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 528 & * ( e3w _n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) &529 & - e3w _n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) )539 & * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 540 & - e3w(ji ,jj,jk,Kmm) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 530 541 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) & 531 & * ( e3w _n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) &532 & - e3w _n(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) )542 & * ( e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 543 & - e3w(ji,jj ,jk,Kmm) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 533 544 ! s-coordinate pressure gradient correction 534 545 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 535 & * ( gde3w _n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj)546 & * ( gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk) ) * r1_e1u(ji,jj) 536 547 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 537 & * ( gde3w _n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj)548 & * ( gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk) ) * r1_e2v(ji,jj) 538 549 ! 539 550 IF( ln_wd_il ) THEN … … 545 556 ! 546 557 ! add to the general momentum trend 547 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap548 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap558 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) + zuap 559 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) + zvap 549 560 END DO 550 561 END DO … … 556 567 557 568 558 SUBROUTINE hpg_isf( kt )569 SUBROUTINE hpg_isf( kt, Kmm, puu, pvv, Krhs ) 559 570 !!--------------------------------------------------------------------- 560 571 !! *** ROUTINE hpg_isf *** … … 568 579 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] 569 580 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ] 570 !! add it to the general momentum trend ( ua,va).571 !! ua = ua- 1/e1u * zhpi572 !! va = va- 1/e2v * zhpj581 !! add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)). 582 !! puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi 583 !! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 573 584 !! iceload is added and partial cell case are added to the top and bottom 574 585 !! 575 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 576 !!---------------------------------------------------------------------- 577 INTEGER, INTENT(in) :: kt ! ocean time-step index 586 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 587 !!---------------------------------------------------------------------- 588 INTEGER , INTENT( in ) :: kt ! ocean time-step index 589 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 590 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 578 591 !! 579 592 INTEGER :: ji, jj, jk, ikt, iktp1i, iktp1j ! dummy loop indices … … 596 609 DO jj = 1, jpj 597 610 ikt = mikt(ji,jj) 598 zts_top(ji,jj,1) = ts n(ji,jj,ikt,1)599 zts_top(ji,jj,2) = ts n(ji,jj,ikt,2)611 zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm) 612 zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm) 600 613 END DO 601 614 END DO … … 612 625 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 613 626 ! we assume ISF is in isostatic equilibrium 614 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w _n(ji+1,jj,iktp1i) &627 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w(ji+1,jj,iktp1i,Kmm) & 615 628 & * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) ) & 616 & - 0.5_wp * e3w _n(ji,jj,ikt) &629 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 617 630 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 618 631 & + ( riceload(ji+1,jj) - riceload(ji,jj)) ) 619 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w _n(ji,jj+1,iktp1j) &632 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w(ji,jj+1,iktp1j,Kmm) & 620 633 & * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 621 & - 0.5_wp * e3w _n(ji,jj,ikt) &634 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 622 635 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 623 636 & + ( riceload(ji,jj+1) - riceload(ji,jj)) ) 624 637 ! s-coordinate pressure gradient correction (=0 if z coordinate) 625 638 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 626 & * ( gde3w _n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)639 & * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj) 627 640 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 628 & * ( gde3w _n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)641 & * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj) 629 642 ! add to the general momentum trend 630 ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1)631 va(ji,jj,1) = va(ji,jj,1) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1)643 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 644 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 632 645 END DO 633 646 END DO … … 641 654 ! hydrostatic pressure gradient along s-surfaces 642 655 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 643 & * ( e3w _n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) &644 & - e3w _n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) )656 & * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 657 & - e3w(ji ,jj,jk,Kmm) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 645 658 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 646 & * ( e3w _n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) &647 & - e3w _n(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) )659 & * ( e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 660 & - e3w(ji,jj ,jk,Kmm) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 648 661 ! s-coordinate pressure gradient correction 649 662 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 650 & * ( gde3w _n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj)663 & * ( gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk) ) / e1u(ji,jj) 651 664 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 652 & * ( gde3w _n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj)665 & * ( gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk) ) / e2v(ji,jj) 653 666 ! add to the general momentum trend 654 ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)655 va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk)667 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 668 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 656 669 END DO 657 670 END DO … … 661 674 662 675 663 SUBROUTINE hpg_djc( kt )676 SUBROUTINE hpg_djc( kt, Kmm, puu, pvv, Krhs ) 664 677 !!--------------------------------------------------------------------- 665 678 !! *** ROUTINE hpg_djc *** … … 669 682 !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003 670 683 !!---------------------------------------------------------------------- 671 INTEGER, INTENT(in) :: kt ! ocean time-step index 684 INTEGER , INTENT( in ) :: kt ! ocean time-step index 685 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 686 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 672 687 !! 673 688 INTEGER :: ji, jj, jk ! dummy loop indices … … 687 702 DO jj = 2, jpjm1 688 703 DO ji = 2, jpim1 689 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji+1,jj) ) > &704 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 690 705 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 691 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) &706 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) & 692 707 & > rn_wdmin1 + rn_wdmin2 693 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( &694 & MAX( ssh n(ji,jj) , sshn(ji+1,jj) ) > &708 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. ( & 709 & MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 695 710 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 696 711 IF(ll_tmp1) THEN 697 712 zcpx(ji,jj) = 1.0_wp 698 713 ELSE IF(ll_tmp2) THEN 699 ! no worries about ssh n(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here700 zcpx(ji,jj) = ABS( (ssh n(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &701 & / (ssh n(ji+1,jj) - sshn(ji ,jj)) )714 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 715 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 716 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 702 717 ELSE 703 718 zcpx(ji,jj) = 0._wp 704 719 END IF 705 720 706 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji,jj+1) ) > &721 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 707 722 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 708 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) &723 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) & 709 724 & > rn_wdmin1 + rn_wdmin2 710 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( &711 & MAX( ssh n(ji,jj) , sshn(ji,jj+1) ) > &725 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. ( & 726 & MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 712 727 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 713 728 … … 715 730 zcpy(ji,jj) = 1.0_wp 716 731 ELSE IF(ll_tmp2) THEN 717 ! no worries about ssh n(ji,jj+1) - sshn(ji,jj) = 0, it won't happen ! here718 zcpy(ji,jj) = ABS( (ssh n(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &719 & / (ssh n(ji,jj+1) - sshn(ji,jj)) )732 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 733 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 734 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 720 735 ELSE 721 736 zcpy(ji,jj) = 0._wp … … 747 762 DO ji = fs_2, fs_jpim1 ! vector opt. 748 763 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1) 749 dzz (ji,jj,jk) = gde3w _n(ji ,jj ,jk) - gde3w_n(ji,jj,jk-1)764 dzz (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji,jj,jk-1) 750 765 drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 751 dzx (ji,jj,jk) = gde3w _n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk )766 dzx (ji,jj,jk) = gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk ) 752 767 drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 753 dzy (ji,jj,jk) = gde3w _n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk )768 dzy (ji,jj,jk) = gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk ) 754 769 END DO 755 770 END DO … … 838 853 DO jj = 2, jpjm1 839 854 DO ji = fs_2, fs_jpim1 ! vector opt. 840 rho_k(ji,jj,1) = -grav * ( e3w _n(ji,jj,1) - gde3w_n(ji,jj,1) ) &855 rho_k(ji,jj,1) = -grav * ( e3w(ji,jj,1,Kmm) - gde3w(ji,jj,1) ) & 841 856 & * ( rhd(ji,jj,1) & 842 857 & + 0.5_wp * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) & 843 & * ( e3w _n (ji,jj,1) - gde3w_n(ji,jj,1) ) &844 & / ( gde3w _n(ji,jj,2) - gde3w_n(ji,jj,1) ) )858 & * ( e3w (ji,jj,1,Kmm) - gde3w(ji,jj,1) ) & 859 & / ( gde3w(ji,jj,2) - gde3w(ji,jj,1) ) ) 845 860 END DO 846 861 END DO … … 854 869 855 870 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) & 856 & * ( gde3w _n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) ) &871 & * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) ) & 857 872 & - grav * z1_10 * ( & 858 873 & ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) ) & 859 & * ( gde3w _n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) &874 & * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) & 860 875 & - ( dzw (ji,jj,jk) - dzw (ji,jj,jk-1) ) & 861 876 & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) & … … 863 878 864 879 rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) & 865 & * ( gde3w _n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) ) &880 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) ) & 866 881 & - grav* z1_10 * ( & 867 882 & ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) ) & 868 & * ( gde3w _n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) &883 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) & 869 884 & - ( dzu (ji+1,jj,jk) - dzu (ji,jj,jk) ) & 870 885 & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) & … … 872 887 873 888 rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) & 874 & * ( gde3w _n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) ) &889 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) ) & 875 890 & - grav* z1_10 * ( & 876 891 & ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) ) & 877 & * ( gde3w _n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) &892 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) & 878 893 & - ( dzv (ji,jj+1,jk) - dzv (ji,jj,jk) ) & 879 894 & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) & … … 897 912 ENDIF 898 913 ! add to the general momentum trend 899 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)900 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)914 puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) 915 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) 901 916 END DO 902 917 END DO … … 920 935 ENDIF 921 936 ! add to the general momentum trend 922 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)923 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)937 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) 938 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) 924 939 END DO 925 940 END DO … … 931 946 932 947 933 SUBROUTINE hpg_prj( kt )948 SUBROUTINE hpg_prj( kt, Kmm, puu, pvv, Krhs ) 934 949 !!--------------------------------------------------------------------- 935 950 !! *** ROUTINE hpg_prj *** … … 940 955 !! all vertical coordinate systems 941 956 !! 942 !! ** Action : - Update ( ua,va) with the now hydrastatic pressure trend957 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 943 958 !!---------------------------------------------------------------------- 944 959 INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear 945 INTEGER, INTENT(in) :: kt ! ocean time-step index 960 INTEGER , INTENT( in ) :: kt ! ocean time-step index 961 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 962 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 946 963 !! 947 964 INTEGER :: ji, jj, jk, jkk ! dummy loop indices … … 975 992 DO jj = 2, jpjm1 976 993 DO ji = 2, jpim1 977 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji+1,jj) ) > &994 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 978 995 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 979 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji+1,jj) + ht_0(ji+1,jj) ) &996 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) & 980 997 & > rn_wdmin1 + rn_wdmin2 981 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 ) .AND. ( &982 & MAX( ssh n(ji,jj) , sshn(ji+1,jj) ) > &998 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. ( & 999 & MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 983 1000 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 984 1001 … … 986 1003 zcpx(ji,jj) = 1.0_wp 987 1004 ELSE IF(ll_tmp2) THEN 988 ! no worries about ssh n(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here989 zcpx(ji,jj) = ABS( (ssh n(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &990 & / (ssh n(ji+1,jj) - sshn(ji ,jj)) )1005 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 1006 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1007 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 991 1008 992 1009 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) … … 995 1012 END IF 996 1013 997 ll_tmp1 = MIN( ssh n(ji,jj) , sshn(ji,jj+1) ) > &1014 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 998 1015 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 999 & MAX( ssh n(ji,jj) + ht_0(ji,jj), sshn(ji,jj+1) + ht_0(ji,jj+1) ) &1016 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) & 1000 1017 & > rn_wdmin1 + rn_wdmin2 1001 ll_tmp2 = ( ABS( ssh n(ji,jj) - sshn(ji,jj+1) ) > 1.E-12 ) .AND. ( &1002 & MAX( ssh n(ji,jj) , sshn(ji,jj+1) ) > &1018 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. ( & 1019 & MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 1003 1020 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1004 1021 … … 1006 1023 zcpy(ji,jj) = 1.0_wp 1007 1024 ELSE IF(ll_tmp2) THEN 1008 ! no worries about ssh n(ji,jj+1) - sshn(ji,jj) = 0, it won't happen ! here1009 zcpy(ji,jj) = ABS( (ssh n(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &1010 & / (ssh n(ji,jj+1) - sshn(ji,jj)) )1025 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 1026 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1027 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 1011 1028 zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 1012 1029 … … 1031 1048 ELSEIF( jk < jpkm1 ) THEN 1032 1049 DO jkk = jk+1, jpk 1033 zrhh(ji,jj,jkk) = interp1(gde3w _n(ji,jj,jkk ), gde3w_n(ji,jj,jkk-1), &1034 & gde3w _n(ji,jj,jkk-2), rhd (ji,jj,jkk-1), rhd(ji,jj,jkk-2))1050 zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk ), gde3w(ji,jj,jkk-1), & 1051 & gde3w(ji,jj,jkk-2), rhd (ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 1035 1052 END DO 1036 1053 ENDIF … … 1041 1058 DO jj = 1, jpj 1042 1059 DO ji = 1, jpi 1043 zdept(ji,jj,1) = 0.5_wp * e3w _n(ji,jj,1) - sshn(ji,jj) * znad1060 zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) * znad 1044 1061 END DO 1045 1062 END DO … … 1048 1065 DO jj = 1, jpj 1049 1066 DO ji = 1, jpi 1050 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w _n(ji,jj,jk)1067 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w(ji,jj,jk,Kmm) 1051 1068 END DO 1052 1069 END DO … … 1065 1082 DO ji = 2, jpi 1066 1083 zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), & 1067 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w _n(ji,jj,1)1084 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm) 1068 1085 1069 1086 ! assuming linear profile across the top half surface layer 1070 zhpi(ji,jj,1) = 0.5_wp * e3w _n(ji,jj,1) * zrhdt11087 zhpi(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) * zrhdt1 1071 1088 END DO 1072 1089 END DO … … 1090 1107 DO ji = 2, jpim1 1091 1108 !!gm BUG ? if it is ssh at u- & v-point then it should be: 1092 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh n(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj)) * &1109 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * & 1093 1110 ! & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1094 ! zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh n(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1)) * &1111 ! zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji,jj+1) * ssh(ji,jj+1,Kmm)) * & 1095 1112 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1096 1113 !!gm not this: 1097 zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh n(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * &1114 zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * & 1098 1115 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1099 zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh n(ji,jj) + e1e2v(ji+1, jj) * sshn(ji,jj+1)) * &1116 zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * & 1100 1117 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1101 1118 END DO … … 1106 1123 DO jj = 2, jpjm1 1107 1124 DO ji = 2, jpim1 1108 zu(ji,jj,1) = - ( e3u _n(ji,jj,1) - zsshu_n(ji,jj) * znad)1109 zv(ji,jj,1) = - ( e3v _n(ji,jj,1) - zsshv_n(ji,jj) * znad)1125 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad) 1126 zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad) 1110 1127 END DO 1111 1128 END DO … … 1114 1131 DO jj = 2, jpjm1 1115 1132 DO ji = 2, jpim1 1116 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u _n(ji,jj,jk)1117 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v _n(ji,jj,jk)1133 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm) 1134 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm) 1118 1135 END DO 1119 1136 END DO … … 1123 1140 DO jj = 2, jpjm1 1124 1141 DO ji = 2, jpim1 1125 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u _n(ji,jj,jk)1126 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v _n(ji,jj,jk)1142 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm) 1143 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm) 1127 1144 END DO 1128 1145 END DO … … 1176 1193 DO WHILE ( -zdept(jid,jj,jk1) < zuijk ) 1177 1194 IF( jk1 == 1 ) THEN 1178 zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh n(jid,jj)*znad)1195 zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad) 1179 1196 zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), & 1180 1197 bsp(jid,jj,1), csp(jid,jj,1), & … … 1196 1213 IF( .NOT.ln_linssh ) THEN 1197 1214 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 1198 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh n(ji+1,jj)-sshn(ji,jj)) )1215 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) ) 1199 1216 ELSE 1200 1217 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) … … 1204 1221 zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 1205 1222 ENDIF 1206 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)1223 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 1207 1224 ENDIF 1208 1225 … … 1234 1251 DO WHILE ( -zdept(ji,jjd,jk1) < zvijk ) 1235 1252 IF( jk1 == 1 ) THEN 1236 zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh n(ji,jjd)*znad)1253 zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad) 1237 1254 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), & 1238 1255 bsp(ji,jjd,1), csp(ji,jjd,1), & … … 1255 1272 IF( .NOT.ln_linssh ) THEN 1256 1273 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 1257 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh n(ji,jj+1)-sshn(ji,jj)) )1274 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) ) 1258 1275 ELSE 1259 1276 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) … … 1264 1281 ENDIF 1265 1282 1266 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk)1283 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk) 1267 1284 ENDIF 1268 1285 ! -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynkeg.F90
r10893 r10928 187 187 ENDIF 188 188 ! 189 IF(ln_ctl) CALL prt_ctl( tab3d_1= ua, clinfo1=' keg - Ua: ', mask1=umask, &190 & tab3d_2= va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )189 IF(ln_ctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' keg - Ua: ', mask1=umask, & 190 & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 191 191 ! 192 192 IF( ln_timing ) CALL timing_stop('dyn_keg') -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf.F90
r10874 r10928 43 43 CONTAINS 44 44 45 SUBROUTINE dyn_ldf( kt )45 SUBROUTINE dyn_ldf( kt, Kbb, Kmm, puu, pvv, Krhs ) 46 46 !!---------------------------------------------------------------------- 47 47 !! *** ROUTINE dyn_ldf *** … … 49 49 !! ** Purpose : compute the lateral ocean dynamics physics. 50 50 !!---------------------------------------------------------------------- 51 INTEGER, INTENT(in) :: kt ! ocean time-step index 51 INTEGER , INTENT( in ) :: kt ! ocean time-step index 52 INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs ! ocean time level indices 53 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 52 54 ! 53 55 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv … … 58 60 IF( l_trddyn ) THEN ! temporary save of momentum trends 59 61 ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 60 ztrdu(:,:,:) = ua(:,:,:)61 ztrdv(:,:,:) = va(:,:,:)62 ztrdu(:,:,:) = puu(:,:,:,Krhs) 63 ztrdv(:,:,:) = pvv(:,:,:,Krhs) 62 64 ENDIF 63 65 64 66 SELECT CASE ( nldf_dyn ) ! compute lateral mixing trend and add it to the general trend 65 67 ! 66 CASE ( np_lap ) ; CALL dyn_ldf_lap( kt, ub, vb, ua, va, 1 ) ! iso-level laplacian 67 CASE ( np_lap_i ) ; CALL dyn_ldf_iso( kt ) ! rotated laplacian 68 CASE ( np_blp ) ; CALL dyn_ldf_blp( kt, ub, vb, ua, va ) ! iso-level bi-laplacian 68 CASE ( np_lap ) 69 CALL dyn_ldf_lap( kt, Kbb, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs), 1 ) ! iso-level laplacian 70 CASE ( np_lap_i ) 71 CALL dyn_ldf_iso( kt, Kbb, Kmm, puu, pvv, Krhs ) ! rotated laplacian 72 CASE ( np_blp ) 73 CALL dyn_ldf_blp( kt, Kbb, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! iso-level bi-laplacian 69 74 ! 70 75 END SELECT 71 76 72 77 IF( l_trddyn ) THEN ! save the horizontal diffusive trends for further diagnostics 73 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)74 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)78 ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 79 ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) 75 80 CALL trd_dyn( ztrdu, ztrdv, jpdyn_ldf, kt ) 76 81 DEALLOCATE ( ztrdu , ztrdv ) 77 82 ENDIF 78 83 ! ! print sum trends (used for debugging) 79 IF(ln_ctl) CALL prt_ctl( tab3d_1= ua, clinfo1=' ldf - Ua: ', mask1=umask, &80 & tab3d_2= va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )84 IF(ln_ctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ldf - Ua: ', mask1=umask, & 85 & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 81 86 ! 82 87 IF( ln_timing ) CALL timing_stop('dyn_ldf') -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_iso.F90
r10425 r10928 60 60 61 61 62 SUBROUTINE dyn_ldf_iso( kt )62 SUBROUTINE dyn_ldf_iso( kt, Kbb, Kmm, puu, pvv, Krhs ) 63 63 !!---------------------------------------------------------------------- 64 64 !! *** ROUTINE dyn_ldf_iso *** … … 81 81 !! horizontal fluxes associated with the rotated lateral mixing: 82 82 !! u-component: 83 !! ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t di[ u b]84 !! - ahmt e2t * mi-1(uslp) dk[ mi(mk(u b)) ]85 !! zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f dj[ u b]86 !! - ahmf e1f * mi(vslp) dk[ mj(mk(u b)) ]83 !! ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t di[ uu ] 84 !! - ahmt e2t * mi-1(uslp) dk[ mi(mk(uu)) ] 85 !! zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f dj[ uu ] 86 !! - ahmf e1f * mi(vslp) dk[ mj(mk(uu)) ] 87 87 !! v-component: 88 !! zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t di[ v b]89 !! - ahmf e2t * mj(uslp) dk[ mi(mk(v b)) ]90 !! zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f dj[ ub]91 !! - ahmt e1f * mj-1(vslp) dk[ mj(mk(v b)) ]88 !! zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t di[ vv ] 89 !! - ahmf e2t * mj(uslp) dk[ mi(mk(vv)) ] 90 !! zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f dj[ vv ] 91 !! - ahmt e1f * mj-1(vslp) dk[ mj(mk(vv)) ] 92 92 !! take the horizontal divergence of the fluxes: 93 93 !! diffu = 1/(e1u*e2u*e3u) { di [ ziut ] + dj-1[ zjuf ] } 94 94 !! diffv = 1/(e1v*e2v*e3v) { di-1[ zivf ] + dj [ zjvt ] } 95 !! Add this trend to the general trend (u a,va):96 !! u a = ua+ diffu95 !! Add this trend to the general trend (uu(rhs),vv(rhs)): 96 !! uu(rhs) = uu(rhs) + diffu 97 97 !! CAUTION: here the isopycnal part is with a coeff. of aht. This 98 98 !! should be modified for applications others than orca_r2 (!!bug) 99 99 !! 100 100 !! ** Action : 101 !! -( ua,va) updated with the before geopotential harmonic mixing trend101 !! -(puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the before geopotential harmonic mixing trend 102 102 !! -(akzu,akzv) to accompt for the diagonal vertical component 103 103 !! of the rotated operator in dynzdf module 104 104 !!---------------------------------------------------------------------- 105 INTEGER, INTENT( in ) :: kt ! ocean time-step index 105 INTEGER , INTENT( in ) :: kt ! ocean time-step index 106 INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs ! ocean time level indices 107 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 106 108 ! 107 109 INTEGER :: ji, jj, jk ! dummy loop indices … … 128 130 DO jj = 2, jpjm1 129 131 DO ji = 2, jpim1 130 uslp (ji,jj,jk) = - ( gdept _b(ji+1,jj,jk) - gdept_b(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)131 vslp (ji,jj,jk) = - ( gdept _b(ji,jj+1,jk) - gdept_b(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)132 wslpi(ji,jj,jk) = - ( gdepw _b(ji+1,jj,jk) - gdepw_b(ji-1,jj,jk) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5133 wslpj(ji,jj,jk) = - ( gdepw _b(ji,jj+1,jk) - gdepw_b(ji,jj-1,jk) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5132 uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 133 vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 134 wslpi(ji,jj,jk) = - ( gdepw(ji+1,jj,jk,Kbb) - gdepw(ji-1,jj,jk,Kbb) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5 135 wslpj(ji,jj,jk) = - ( gdepw(ji,jj+1,jk,Kbb) - gdepw(ji,jj-1,jk,Kbb) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5 134 136 END DO 135 137 END DO … … 151 153 ! zdkv(jk=1)=zdkv(jk=2) 152 154 153 zdk1u(:,:) = ( ub(:,:,jk) -ub(:,:,jk+1) ) * umask(:,:,jk+1)154 zdk1v(:,:) = ( vb(:,:,jk) -vb(:,:,jk+1) ) * vmask(:,:,jk+1)155 zdk1u(:,:) = ( puu(:,:,jk,Kbb) -puu(:,:,jk+1,Kbb) ) * umask(:,:,jk+1) 156 zdk1v(:,:) = ( pvv(:,:,jk,Kbb) -pvv(:,:,jk+1,Kbb) ) * vmask(:,:,jk+1) 155 157 156 158 IF( jk == 1 ) THEN … … 158 160 zdkv(:,:) = zdk1v(:,:) 159 161 ELSE 160 zdku(:,:) = ( ub(:,:,jk-1) - ub(:,:,jk) ) * umask(:,:,jk)161 zdkv(:,:) = ( vb(:,:,jk-1) - vb(:,:,jk) ) * vmask(:,:,jk)162 zdku(:,:) = ( puu(:,:,jk-1,Kbb) - puu(:,:,jk,Kbb) ) * umask(:,:,jk) 163 zdkv(:,:) = ( pvv(:,:,jk-1,Kbb) - pvv(:,:,jk,Kbb) ) * vmask(:,:,jk) 162 164 ENDIF 163 165 … … 171 173 DO jj = 2, jpjm1 172 174 DO ji = fs_2, jpi ! vector opt. 173 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( e3u _n(ji,jj,jk), e3u_n(ji-1,jj,jk) ) * r1_e1t(ji,jj)175 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( e3u(ji,jj,jk,Kmm), e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj) 174 176 175 177 zmskt = 1._wp / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & … … 178 180 zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 179 181 180 ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) &182 ziut(ji,jj) = ( zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) ) & 181 183 & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & 182 184 & +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk) … … 186 188 DO jj = 2, jpjm1 187 189 DO ji = fs_2, jpi ! vector opt. 188 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t _n(ji,jj,jk) * r1_e1t(ji,jj)190 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj) 189 191 190 192 zmskt = 1._wp / MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk+1) & … … 193 195 zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 194 196 195 ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) &197 ziut(ji,jj) = ( zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) ) & 196 198 & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & 197 199 & +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk) … … 203 205 DO jj = 1, jpjm1 204 206 DO ji = 1, fs_jpim1 ! vector opt. 205 zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f _n(ji,jj,jk) * r1_e2f(ji,jj)207 zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj) 206 208 207 209 zmskf = 1._wp / MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & … … 210 212 zcof2 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 211 213 212 zjuf(ji,jj) = ( zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) &214 zjuf(ji,jj) = ( zabe2 * ( puu(ji,jj+1,jk,Kbb) - puu(ji,jj,jk,Kbb) ) & 213 215 & + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) & 214 216 & +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) * fmask(ji,jj,jk) … … 224 226 DO jj = 2, jpjm1 225 227 DO ji = 1, fs_jpim1 ! vector opt. 226 zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f _n(ji,jj,jk) * r1_e1f(ji,jj)228 zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj) 227 229 228 230 zmskf = 1._wp / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & … … 231 233 zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 232 234 233 zivf(ji,jj) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) &235 zivf(ji,jj) = ( zabe1 * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji,jj,jk,Kbb) ) & 234 236 & + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj) & 235 237 & + zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk) … … 241 243 DO jj = 2, jpj 242 244 DO ji = 1, fs_jpim1 ! vector opt. 243 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( e3v _n(ji,jj,jk), e3v_n(ji,jj-1,jk) ) * r1_e2t(ji,jj)245 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( e3v(ji,jj,jk,Kmm), e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj) 244 246 245 247 zmskt = 1._wp / MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & … … 248 250 zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 249 251 250 zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) &252 zjvt(ji,jj) = ( zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) ) & 251 253 & + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) & 252 254 & +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk) … … 256 258 DO jj = 2, jpj 257 259 DO ji = 1, fs_jpim1 ! vector opt. 258 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t _n(ji,jj,jk) * r1_e2t(ji,jj)260 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj) 259 261 260 262 zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & … … 263 265 zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 264 266 265 zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) &267 zjvt(ji,jj) = ( zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) ) & 266 268 & + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) & 267 269 & +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk) … … 275 277 DO jj = 2, jpjm1 276 278 DO ji = 2, jpim1 !!gm Question vectop possible??? !!bug 277 ua(ji,jj,jk) = ua(ji,jj,jk) + ( ziut(ji+1,jj) - ziut(ji,jj ) &278 & + zjuf(ji ,jj) - zjuf(ji,jj-1) ) * r1_e1e2u(ji,jj) / e3u _n(ji,jj,jk)279 va(ji,jj,jk) = va(ji,jj,jk) + ( zivf(ji,jj ) - zivf(ji-1,jj) &280 & + zjvt(ji,jj+1) - zjvt(ji,jj ) ) * r1_e1e2v(ji,jj) / e3v _n(ji,jj,jk)279 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( ziut(ji+1,jj) - ziut(ji,jj ) & 280 & + zjuf(ji ,jj) - zjuf(ji,jj-1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 281 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zivf(ji,jj ) - zivf(ji-1,jj) & 282 & + zjvt(ji,jj+1) - zjvt(ji,jj ) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 281 283 END DO 282 284 END DO … … 286 288 287 289 ! print sum trends (used for debugging) 288 IF(ln_ctl) CALL prt_ctl( tab3d_1= ua, clinfo1=' ldfh - Ua: ', mask1=umask, &289 & tab3d_2= va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )290 IF(ln_ctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ldfh - Ua: ', mask1=umask, & 291 & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 290 292 291 293 … … 306 308 DO ji = 2, jpi 307 309 ! i-gradient of u at jj 308 zdiu (ji,jk) = tmask(ji,jj ,jk) * ( ub(ji,jj ,jk) - ub(ji-1,jj ,jk) )310 zdiu (ji,jk) = tmask(ji,jj ,jk) * ( puu(ji,jj ,jk,Kbb) - puu(ji-1,jj ,jk,Kbb) ) 309 311 ! j-gradient of u and v at jj 310 zdju (ji,jk) = fmask(ji,jj ,jk) * ( ub(ji,jj+1,jk) - ub(ji ,jj ,jk) )311 zdjv (ji,jk) = tmask(ji,jj ,jk) * ( vb(ji,jj ,jk) - vb(ji ,jj-1,jk) )312 zdju (ji,jk) = fmask(ji,jj ,jk) * ( puu(ji,jj+1,jk,Kbb) - puu(ji ,jj ,jk,Kbb) ) 313 zdjv (ji,jk) = tmask(ji,jj ,jk) * ( pvv(ji,jj ,jk,Kbb) - pvv(ji ,jj-1,jk,Kbb) ) 312 314 ! j-gradient of u and v at jj+1 313 zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj ,jk) - ub(ji ,jj-1,jk) )314 zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji ,jj ,jk) )315 zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( puu(ji,jj ,jk,Kbb) - puu(ji ,jj-1,jk,Kbb) ) 316 zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pvv(ji,jj+1,jk,Kbb) - pvv(ji ,jj ,jk,Kbb) ) 315 317 END DO 316 318 END DO … … 318 320 DO ji = 1, jpim1 319 321 ! i-gradient of v at jj 320 zdiv (ji,jk) = fmask(ji,jj ,jk) * ( vb(ji+1,jj,jk) - vb(ji ,jj ,jk) )322 zdiv (ji,jk) = fmask(ji,jj ,jk) * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji ,jj ,jk,Kbb) ) 321 323 END DO 322 324 END DO … … 391 393 DO jk = 1, jpkm1 392 394 DO ji = 2, jpim1 393 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)394 va(ji,jj,jk) = va(ji,jj,jk) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)395 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 396 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 395 397 END DO 396 398 END DO -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_lap_blp.F90
r10874 r10928 35 35 CONTAINS 36 36 37 SUBROUTINE dyn_ldf_lap( kt, pub, pvb, pua, pva, kpass )37 SUBROUTINE dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs, kpass ) 38 38 !!---------------------------------------------------------------------- 39 39 !! *** ROUTINE dyn_ldf_lap *** … … 45 45 !! writen as : grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) 46 46 !! 47 !! ** Action : - pu a, pva increased by the harmonic operator applied on pub, pvb.47 !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv. 48 48 !!---------------------------------------------------------------------- 49 49 INTEGER , INTENT(in ) :: kt ! ocean time-step index 50 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 50 51 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 51 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu b, pvb! before velocity [m/s]52 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu a, pva! velocity trend [m/s2]52 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu, pv ! before velocity [m/s] 53 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! velocity trend [m/s2] 53 54 ! 54 55 INTEGER :: ji, jj, jk ! dummy loop indices … … 76 77 !!gm open question here : e3f at before or now ? probably now... 77 78 !!gm note that ahmf has already been multiplied by fmask 78 zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f _n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) &79 & * ( e2v(ji ,jj-1) * pv b(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk) &80 & - e1u(ji-1,jj ) * pu b(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk) )79 zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & 80 & * ( e2v(ji ,jj-1) * pv(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk) & 81 & - e1u(ji-1,jj ) * pu(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk) ) 81 82 ! ! ahm * div (computed from 2 to jpi/jpj) 82 83 !!gm note that ahmt has already been multiplied by tmask 83 zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t _b(ji,jj,jk) &84 & * ( e2u(ji,jj)*e3u _b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk) &85 & + e1v(ji,jj)*e3v _b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk) )84 zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb) & 85 & * ( e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk) & 86 & + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk) ) 86 87 END DO 87 88 END DO … … 89 90 DO jj = 2, jpjm1 ! - curl( curl) + grad( div ) 90 91 DO ji = fs_2, fs_jpim1 ! vector opt. 91 pu a(ji,jj,jk) = pua(ji,jj,jk) + zsign * ( &92 & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u _n(ji,jj,jk) &92 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * ( & 93 & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & 93 94 & + ( zdiv(ji+1,jj) - zdiv(ji,jj ) ) * r1_e1u(ji,jj) ) 94 95 ! 95 pv a(ji,jj,jk) = pva(ji,jj,jk) + zsign * ( &96 & ( zcur(ji,jj ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v _n(ji,jj,jk) &96 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * ( & 97 & ( zcur(ji,jj ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm) & 97 98 & + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) ) 98 99 END DO … … 105 106 106 107 107 SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva)108 SUBROUTINE dyn_ldf_blp( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs ) 108 109 !!---------------------------------------------------------------------- 109 110 !! *** ROUTINE dyn_ldf_blp *** … … 116 117 !! It is computed by two successive calls to dyn_ldf_lap routine 117 118 !! 118 !! ** Action : pt aupdated with the before rotated bilaplacian diffusion119 !! ** Action : pt(:,:,:,:,Krhs) updated with the before rotated bilaplacian diffusion 119 120 !!---------------------------------------------------------------------- 120 121 INTEGER , INTENT(in ) :: kt ! ocean time-step index 121 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pub, pvb ! before velocity fields 122 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! momentum trend 122 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 123 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu, pv ! before velocity fields 124 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! momentum trend 123 125 ! 124 126 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zulap, zvlap ! laplacian at u- and v-point … … 134 136 zvlap(:,:,:) = 0._wp 135 137 ! 136 CALL dyn_ldf_lap( kt, pub, pvb, zulap, zvlap, 1 ) ! rotated laplacian applied to ptb (output in zlap)138 CALL dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, zulap, zvlap, 1 ) ! rotated laplacian applied to pt (output in zlap,Kbb) 137 139 ! 138 140 CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1., zvlap, 'V', -1. ) ! Lateral boundary conditions 139 141 ! 140 CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 ) ! rotated laplacian applied to zlap (output in pta)142 CALL dyn_ldf_lap( kt, Kbb, Kmm, zulap, zvlap, pu_rhs, pv_rhs, 2 ) ! rotated laplacian applied to zlap (output in pt(:,:,:,:,Krhs)) 141 143 ! 142 144 END SUBROUTINE dyn_ldf_blp -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynspg.F90
r10919 r10928 175 175 ENDIF 176 176 ! ! print mean trends (used for debugging) 177 IF(ln_ctl) CALL prt_ctl( tab3d_1= ua, clinfo1=' spg - Ua: ', mask1=umask, &178 & tab3d_2= va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )177 IF(ln_ctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' spg - Ua: ', mask1=umask, & 178 & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 179 179 ! 180 180 IF( ln_timing ) CALL timing_stop('dyn_spg') -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90
r10893 r10928 181 181 ! 182 182 ! ! print sum trends (used for debugging) 183 IF(ln_ctl) CALL prt_ctl( tab3d_1= ua, clinfo1=' vor - Ua: ', mask1=umask, &184 & tab3d_2= va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )183 IF(ln_ctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor - Ua: ', mask1=umask, & 184 & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 185 185 ! 186 186 IF( ln_timing ) CALL timing_stop('dyn_vor') -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzad.F90
r10893 r10928 116 116 ENDIF 117 117 ! ! Control print 118 IF(ln_ctl) CALL prt_ctl( tab3d_1= ua, clinfo1=' zad - Ua: ', mask1=umask, &119 & tab3d_2= va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )118 IF(ln_ctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' zad - Ua: ', mask1=umask, & 119 & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 120 120 ! 121 121 IF( ln_timing ) CALL timing_stop('dyn_zad') -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90
r10893 r10928 130 130 IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 131 131 DO jk = 1, jpkm1 ! remove barotropic velocities 132 puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - u a_b(:,:) ) * umask(:,:,jk)133 pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - v a_b(:,:) ) * vmask(:,:,jk)132 puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - uu_b(:,:,Kaa) ) * umask(:,:,jk) 133 pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - vv_b(:,:,Kaa) ) * vmask(:,:,jk) 134 134 END DO 135 135 DO jj = 2, jpjm1 ! Add bottom/top stress due to barotropic component only … … 139 139 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 140 140 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 141 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * u a_b(ji,jj) / ze3ua142 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * v a_b(ji,jj) / ze3va141 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 142 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 143 143 END DO 144 144 END DO … … 150 150 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 151 151 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 152 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * u a_b(ji,jj) / ze3ua153 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * v a_b(ji,jj) / ze3va152 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 153 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 154 154 END DO 155 155 END DO … … 494 494 ENDIF 495 495 ! ! print mean trends (used for debugging) 496 IF(ln_ctl) CALL prt_ctl( tab3d_1= ua, clinfo1=' zdf - Ua: ', mask1=umask, &497 & tab3d_2= va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )496 IF(ln_ctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' zdf - Ua: ', mask1=umask, & 497 & tab3d_2=pvv(:,:,:,Kaa), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 498 498 ! 499 499 IF( ln_timing ) CALL timing_stop('dyn_zdf') -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90
r10927 r10928 428 428 ! 429 429 IF( ln_diurnal_only ) THEN ! diurnal only: a subset of the initialisation routines 430 CALL istate_init( Nbb, Nnn )! ocean initial state (Dynamics and tracers)430 CALL istate_init( Nbb, Nnn, Naa ) ! ocean initial state (Dynamics and tracers) 431 431 CALL sbc_init( Nbb, Nnn ) ! Forcings : surface module 432 432 CALL tra_qsr_init ! penetrative solar radiation qsr … … 440 440 ENDIF 441 441 442 CALL istate_init( Nbb, Nnn )! ocean initial state (Dynamics and tracers)442 CALL istate_init( Nbb, Nnn, Naa ) ! ocean initial state (Dynamics and tracers) 443 443 444 444 ! ! external forcing … … 464 464 465 465 ! ! Dynamics 466 IF( lk_c1d ) CALL dyn_dmp_init ! internal momentum damping467 CALL dyn_adv_init ! advection (vector or flux form)468 CALL dyn_vor_init ! vorticity term including Coriolis469 CALL dyn_ldf_init ! lateral mixing470 CALL dyn_hpg_init 471 CALL dyn_spg_init ! surface pressure gradient466 IF( lk_c1d ) CALL dyn_dmp_init ! internal momentum damping 467 CALL dyn_adv_init ! advection (vector or flux form) 468 CALL dyn_vor_init ! vorticity term including Coriolis 469 CALL dyn_ldf_init ! lateral mixing 470 CALL dyn_hpg_init( Nnn ) ! horizontal gradient of Hydrostatic pressure 471 CALL dyn_spg_init ! surface pressure gradient 472 472 473 473 #if defined key_top -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90
r10927 r10928 195 195 CALL dyn_adv( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS 196 196 CALL dyn_vor( kstp, Nnn , uu, vv, Nrhs ) ! vorticity ==> RHS 197 CALL dyn_ldf ( kstp) ! lateral mixing197 CALL dyn_ldf( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! lateral mixing 198 198 IF( ln_zdfosm ) CALL dyn_osm( kstp, Nnn , uu, vv, Nrhs ) ! OSMOSIS non-local velocity fluxes ==> RHS 199 CALL dyn_hpg ( kstp) ! horizontal gradient of Hydrostatic pressure200 CALL dyn_spg 199 CALL dyn_hpg( kstp, Nnn, uu, vv, Nrhs ) ! horizontal gradient of Hydrostatic pressure 200 CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa ) ! surface pressure gradient 201 201 202 202 ! With split-explicit free surface, since now transports have been updated and ssha as well -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90
r10922 r10928 315 315 IF( ln_ctl ) CALL prt_ctl_init ! Print control 316 316 317 CALL istate_init ! ocean initial state (Dynamics and tracers)317 CALL istate_init( Nnn, Naa ) ! ocean initial state (Dynamics and tracers) 318 318 319 319 CALL sbc_init( Nbb, Nnn ) ! Forcings : surface module … … 524 524 END SUBROUTINE nemo_set_cfctl 525 525 526 SUBROUTINE istate_init 526 SUBROUTINE istate_init( Kmm, Kaa ) 527 527 !!---------------------------------------------------------------------- 528 528 !! *** ROUTINE istate_init *** … … 530 530 !! ** Purpose : Initialization to zero of the dynamics and tracers. 531 531 !!---------------------------------------------------------------------- 532 INTEGER, INTENT(in) :: Kmm, Kaa ! ocean time level indices 532 533 ! 533 534 ! now fields ! after fields ! 534 u n (:,:,:) = 0._wp ; ua(:,:,:) = 0._wp !535 v n (:,:,:) = 0._wp ; va(:,:,:) = 0._wp !536 w n(:,:,:) = 0._wp ! !535 uu (:,:,:,Kmm) = 0._wp ; uu(:,:,:,Kaa) = 0._wp ! 536 vv (:,:,:,Kmm) = 0._wp ; vv(:,:,:,Kaa) = 0._wp ! 537 ww (:,:,:) = 0._wp ! ! 537 538 hdivn(:,:,:) = 0._wp ! ! 538 ts n (:,:,:,:) = 0._wp ! !539 ts (:,:,:,:,Kmm) = 0._wp ! ! 539 540 ! 540 541 rhd (:,:,:) = 0.e0
Note: See TracChangeset
for help on using the changeset viewer.
