Changeset 14058 for NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/domqco.F90
- Timestamp:
- 2020-12-03T16:53:41+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/domqco.F90
r14046 r14058 8 8 !! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates 9 9 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 10 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping11 !! 4.x ! 2020-02 (G. Madec, S. Techene) pure z* (quasi-eulerian) coordinate12 !!---------------------------------------------------------------------- 13 14 !!---------------------------------------------------------------------- 15 !! dom_q e_init: define initial vertical scale factors, depths and column thickness16 !! dom_q e_r3c : Compute ssh/h_0 ratioat t-, u-, v-, and optionally f-points17 !! qe_rst_read : read/write restart file18 !! dom_qe_ctl: Check the vvl options10 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) add time level indices for prognostic variables 11 !! - ! 2020-02 (S. Techene, G. Madec) quasi-eulerian coordinate (z* or s*) 12 !!---------------------------------------------------------------------- 13 14 !!---------------------------------------------------------------------- 15 !! dom_qco_init : define initial vertical scale factors, depths and column thickness 16 !! dom_qco_zgr : Set ssh/h_0 ratio at t 17 !! dom_qco_r3c : Compute ssh/h_0 ratio at t-, u-, v-, and optionally f-points 18 !! qco_ctl : Check the vvl options 19 19 !!---------------------------------------------------------------------- 20 20 USE oce ! ocean dynamics and tracers … … 55 55 LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints 56 56 57 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport58 59 57 !! * Substitutions 60 58 # include "do_loop_substitute.h90" … … 79 77 !! 80 78 !!---------------------------------------------------------------------- 81 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 79 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! time level indices 80 !!---------------------------------------------------------------------- 82 81 ! 83 82 IF(lwp) WRITE(numout,*) … … 85 84 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 86 85 ! 87 CALL dom_qco_ctl ! choose vertical coordinate (z_star, z_tilde or layer) 88 ! 89 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 90 CALL qe_rst_read( nit000, Kbb, Kmm ) 91 ! 92 CALL dom_qco_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column 86 CALL qco_ctl ! choose vertical coordinate (z_star, z_tilde or layer) 87 ! 88 CALL dom_qco_zgr( Kbb, Kmm ) ! interpolation scale factor, depth and water column 89 ! 90 #if defined key_agrif 91 ! We need to define r3[tuv](Kaa) for AGRIF initialisation (should not be a 92 ! problem for the restartability...) 93 r3t(:,:,Kaa) = r3t(:,:,Kmm) 94 r3u(:,:,Kaa) = r3u(:,:,Kmm) 95 r3v(:,:,Kaa) = r3v(:,:,Kmm) 96 #endif 93 97 ! 94 98 END SUBROUTINE dom_qco_init 95 99 96 100 97 SUBROUTINE dom_qco_zgr( Kbb, Kmm, Kaa)101 SUBROUTINE dom_qco_zgr( Kbb, Kmm ) 98 102 !!---------------------------------------------------------------------- 99 103 !! *** ROUTINE dom_qco_init *** 100 104 !! 101 !! ** Purpose : Initialization of all ssh. to h._0 ratio 102 !! 103 !! ** Method : - interpolate scale factors 104 !! 105 !! ** Action : - r3(t/u/v)_b 106 !! - r3(t/u/v/f)_n 107 !! 108 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 109 !!---------------------------------------------------------------------- 110 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 105 !! ** Purpose : Initialization of all r3. = ssh./h._0 ratios 106 !! 107 !! ** Method : Call domqco using Kbb and Kmm 108 !! NB: dom_qco_zgr is called by dom_qco_init it uses ssh from ssh_init 109 !! 110 !! ** Action : - r3(t/u/v)(Kbb) 111 !! - r3(t/u/v/f)(Kmm) 112 !!---------------------------------------------------------------------- 113 INTEGER, INTENT(in) :: Kbb, Kmm ! time level indices 111 114 !!---------------------------------------------------------------------- 112 115 ! 113 116 ! !== Set of all other vertical scale factors ==! (now and before) 114 117 ! ! Horizontal interpolation of e3t 115 CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb) )118 CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb) ) 116 119 CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) ) 117 120 ! … … 143 146 ! !== ratio at u-,v-point ==! 144 147 ! 145 IF( ln_dynadv_vec ) THEN !- Vector Form (thickness weighted averaging) 148 !!st IF( ln_dynadv_vec ) THEN !- Vector Form (thickness weighted averaging) 149 #if ! defined key_qcoTest_FluxForm 150 ! ! no 'key_qcoTest_FluxForm' : surface weighted ssh average 146 151 DO_2D( 0, 0, 0, 0 ) 147 152 pr3u(ji,jj) = 0.5_wp * ( e1e2t(ji ,jj) * pssh(ji ,jj) & … … 150 155 & + e1e2t(ji,jj+1) * pssh(ji,jj+1) ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj) 151 156 END_2D 152 ELSE !- Flux Form (simple averaging) 157 !!st ELSE !- Flux Form (simple averaging) 158 #else 153 159 DO_2D( 0, 0, 0, 0 ) 154 pr3u(ji,jj) = 0.5_wp * ( pssh(ji ,jj) + pssh(ji+1,jj) ) * r1_hu_0(ji,jj)155 pr3v(ji,jj) = 0.5_wp * ( pssh(ji,jj ) + pssh(ji,jj+1) ) * r1_hv_0(ji,jj)160 pr3u(ji,jj) = 0.5_wp * ( pssh(ji,jj) + pssh(ji+1,jj ) ) * r1_hu_0(ji,jj) 161 pr3v(ji,jj) = 0.5_wp * ( pssh(ji,jj) + pssh(ji ,jj+1) ) * r1_hv_0(ji,jj) 156 162 END_2D 157 ENDIF 163 !!st ENDIF 164 #endif 158 165 ! 159 166 IF( .NOT.PRESENT( pr3f ) ) THEN !- lbc on ratio at u-, v-points only … … 163 170 ELSE !== ratio at f-point ==! 164 171 ! 165 IF( ln_dynadv_vec ) THEN !- Vector Form (thickness weighted averaging) 166 DO_2D( 1, 0, 1, 0 ) ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line 172 !!st IF( ln_dynadv_vec ) THEN !- Vector Form (thickness weighted averaging) 173 #if ! defined key_qcoTest_FluxForm 174 ! ! no 'key_qcoTest_FluxForm' : surface weighted ssh average 175 176 DO_2D( 0, 0, 0, 0 ) ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line 167 177 pr3f(ji,jj) = 0.25_wp * ( e1e2t(ji ,jj ) * pssh(ji ,jj ) & 168 178 & + e1e2t(ji+1,jj ) * pssh(ji+1,jj ) & … … 170 180 & + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1) ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj) 171 181 END_2D 172 ELSE !- Flux Form (simple averaging) 173 DO_2D( 1, 0, 1, 0 ) ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line 174 pr3f(ji,jj) = 0.25_wp * ( pssh(ji ,jj ) + pssh(ji+1,jj ) & 175 & + pssh(ji ,jj+1) + pssh(ji+1,jj+1) ) * r1_hf_0(ji,jj) 182 !!st ELSE !- Flux Form (simple averaging) 183 #else 184 DO_2D( 0, 0, 0, 0 ) ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line 185 pr3f(ji,jj) = 0.25_wp * ( pssh(ji,jj ) + pssh(ji+1,jj ) & 186 & + pssh(ji,jj+1) + pssh(ji+1,jj+1) ) * r1_hf_0(ji,jj) 176 187 END_2D 177 ENDIF 188 !!st ENDIF 189 #endif 178 190 ! ! lbc on ratio at u-,v-,f-points 179 191 CALL lbc_lnk_multi( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp, pr3f, 'F', 1._wp ) … … 184 196 185 197 186 SUBROUTINE q e_rst_read( kt, Kbb, Kmm )198 SUBROUTINE qco_ctl 187 199 !!--------------------------------------------------------------------- 188 !! *** ROUTINE qe_rst_read *** 189 !! 190 !! ** Purpose : Read ssh in restart file 191 !! 192 !! ** Method : use of IOM library 193 !! if the restart does not contain ssh, 194 !! it is set to the _0 values. 195 !!---------------------------------------------------------------------- 196 INTEGER , INTENT(in) :: kt ! ocean time-step 197 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 198 ! 199 INTEGER :: ji, jj, jk 200 INTEGER :: id1, id2 ! local integers 201 !!---------------------------------------------------------------------- 202 ! 203 IF( ln_rstart ) THEN !* Read the restart file 204 CALL rst_read_open ! open the restart file if necessary 205 ! 206 id1 = iom_varid( numror, 'sshb', ldstop = .FALSE. ) 207 id2 = iom_varid( numror, 'sshn', ldstop = .FALSE. ) 208 ! 209 ! ! --------- ! 210 ! ! all cases ! 211 ! ! --------- ! 212 ! 213 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 214 CALL iom_get( numror, jpdom_auto, 'sshb' , ssh(:,:,Kbb) ) 215 CALL iom_get( numror, jpdom_auto, 'sshn' , ssh(:,:,Kmm) ) 216 ! needed to restart if land processor not computed 217 IF(lwp) write(numout,*) 'qe_rst_read : ssh(:,:,Kbb) and ssh(:,:,Kmm) found in restart files' 218 WHERE ( ssmask(:,:) == 0.0_wp ) !!gm/st ==> sm should not be necessary on ssh when it was required on e3 219 ssh(:,:,Kmm) = 0._wp 220 ssh(:,:,Kbb) = 0._wp 221 END WHERE 222 IF( l_1st_euler ) THEN 223 ssh(:,:,Kbb) = ssh(:,:,Kmm) 224 ENDIF 225 ELSE IF( id1 > 0 ) THEN 226 IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kmm) not found in restart files' 227 IF(lwp) write(numout,*) 'sshn set equal to sshb.' 228 IF(lwp) write(numout,*) 'neuler is forced to 0' 229 CALL iom_get( numror, jpdom_auto, 'sshb', ssh(:,:,Kbb) ) 230 ssh(:,:,Kmm) = ssh(:,:,Kbb) 231 l_1st_euler = .TRUE. 232 ELSE IF( id2 > 0 ) THEN 233 IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kbb) not found in restart files' 234 IF(lwp) write(numout,*) 'sshb set equal to sshn.' 235 IF(lwp) write(numout,*) 'neuler is forced to 0' 236 CALL iom_get( numror, jpdom_auto, 'sshn', ssh(:,:,Kmm) ) 237 ssh(:,:,Kbb) = ssh(:,:,Kmm) 238 l_1st_euler = .TRUE. 239 ELSE 240 IF(lwp) write(numout,*) 'qe_rst_read WARNING : ssh(:,:,Kmm) not found in restart file' 241 IF(lwp) write(numout,*) 'ssh_b and ssh_n set to zero' 242 IF(lwp) write(numout,*) 'neuler is forced to 0' 243 ssh(:,:,:) = 0._wp 244 l_1st_euler = .TRUE. 245 ENDIF 246 ! 247 ELSE !* Initialize at "rest" 248 ! 249 IF( ll_wd ) THEN ! MJB ll_wd edits start here - these are essential 250 ! 251 IF( cn_cfg == 'wad' ) THEN ! Wetting and drying test case 252 CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 253 ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones 254 ssh(:,: ,Kmm) = ssh(:,: ,Kbb) 255 uu (:,:,: ,Kmm) = uu (:,:,: ,Kbb) 256 vv (:,:,: ,Kmm) = vv (:,:,: ,Kbb) 257 ELSE ! if not test case 258 ssh(:,:,Kmm) = -ssh_ref 259 ssh(:,:,Kbb) = -ssh_ref 260 ! 261 DO_2D( 1, 1, 1, 1 ) 262 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 263 ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 264 ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 265 ENDIF 266 END_2D 267 ENDIF 268 269 DO ji = 1, jpi 270 DO jj = 1, jpj 271 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 272 CALL ctl_stop( 'qe_rst_read: ht_0 must be positive at potentially wet points' ) 273 ENDIF 274 END DO 275 END DO 276 ! 277 ELSE 278 ! 279 ! Just to read set ssh in fact, called latter once vertical grid 280 ! is set up: 281 ! CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 282 ! ! 283 ssh(:,:,:) = 0._wp 284 ! 285 ENDIF ! end of ll_wd edits 286 ! 287 ENDIF 288 ! 289 END SUBROUTINE qe_rst_read 290 291 292 SUBROUTINE dom_qco_ctl 293 !!--------------------------------------------------------------------- 294 !! *** ROUTINE dom_qco_ctl *** 200 !! *** ROUTINE qco_ctl *** 295 201 !! 296 202 !! ** Purpose : Control the consistency between namelist options … … 312 218 IF(lwp) THEN ! Namelist print 313 219 WRITE(numout,*) 314 WRITE(numout,*) ' dom_qco_ctl : choice/control of the variable vertical coordinate'315 WRITE(numout,*) '~~~~~~~~ ~~~'220 WRITE(numout,*) 'qco_ctl : choice/control of the variable vertical coordinate' 221 WRITE(numout,*) '~~~~~~~~' 316 222 WRITE(numout,*) ' Namelist nam_vvl : chose a vertical coordinate' 317 223 WRITE(numout,*) ' zstar ln_vvl_zstar = ', ln_vvl_zstar … … 357 263 #endif 358 264 ! 359 END SUBROUTINE dom_qco_ctl265 END SUBROUTINE qco_ctl 360 266 361 267 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
