- Timestamp:
- 2017-01-02T11:06:49+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90
r5767 r7522 22 22 USE c1d ! 1D configuration: lk_c1d 23 23 USE dom_oce ! ocean domain: variables 24 USE domvvl ! variable volume 24 25 USE zdf_oce ! ocean vertical physics: variables 25 26 USE sbc_oce ! surface module: variables … … 28 29 USE trabbl ! active tracer: bottom boundary layer 29 30 USE ldfslp ! lateral diffusion: iso-neutral slopes 31 USE sbcrnf ! river runoffs 30 32 USE ldfeiv ! eddy induced velocity coef. 31 33 USE ldftra_oce ! ocean tracer lateral physics … … 39 41 USE prtctl ! print control 40 42 USE fldread ! read input fields 43 USE wrk_nemo ! Memory allocation 41 44 USE timing ! Timing 45 USE trc, ONLY : ln_rsttr, numrtr, numrtw, lrst_trc 42 46 43 47 IMPLICIT NONE … … 46 50 PUBLIC dta_dyn_init ! called by opa.F90 47 51 PUBLIC dta_dyn ! called by step.F90 48 49 CHARACTER(len=100) :: cn_dir !: Root directory for location of ssr files 50 LOGICAL :: ln_dynwzv !: vertical velocity read in a file (T) or computed from u/v (F) 51 LOGICAL :: ln_dynbbl !: bbl coef read in a file (T) or computed (F) 52 LOGICAL :: ln_degrad !: degradation option enabled or not 53 LOGICAL :: ln_dynrnf !: read runoff data in file (T) or set to zero (F) 54 55 INTEGER , PARAMETER :: jpfld = 21 ! maximum number of fields to read 52 PUBLIC dta_dyn_swp ! called by step.F90 53 54 CHARACTER(len=100) :: cn_dir !: Root directory for location of ssr files 55 LOGICAL :: ln_dynrnf !: read runoff data in file (T) or set to zero (F) 56 LOGICAL :: ln_dynrnf_depth !: read runoff data in file (T) or set to zero (F) 57 REAL(wp) :: fwbcorr 58 59 60 INTEGER , PARAMETER :: jpfld = 20 ! maximum number of fields to read 56 61 INTEGER , SAVE :: jf_tem ! index of temperature 57 62 INTEGER , SAVE :: jf_sal ! index of salinity 58 INTEGER , SAVE :: jf_uwd ! index of u- wind59 INTEGER , SAVE :: jf_vwd ! index of v- wind60 INTEGER , SAVE :: jf_wwd ! index of w-wind63 INTEGER , SAVE :: jf_uwd ! index of u-transport 64 INTEGER , SAVE :: jf_vwd ! index of v-transport 65 INTEGER , SAVE :: jf_wwd ! index of v-transport 61 66 INTEGER , SAVE :: jf_avt ! index of Kz 62 67 INTEGER , SAVE :: jf_mld ! index of mixed layer deptht 63 68 INTEGER , SAVE :: jf_emp ! index of water flux 69 INTEGER , SAVE :: jf_empb ! index of water flux 64 70 INTEGER , SAVE :: jf_qsr ! index of solar radiation 65 71 INTEGER , SAVE :: jf_wnd ! index of wind speed 66 72 INTEGER , SAVE :: jf_ice ! index of sea ice cover 67 73 INTEGER , SAVE :: jf_rnf ! index of river runoff 74 INTEGER , SAVE :: jf_fmf ! index of downward salt flux 68 75 INTEGER , SAVE :: jf_ubl ! index of u-bbl coef 69 76 INTEGER , SAVE :: jf_vbl ! index of v-bbl coef 70 INTEGER , SAVE :: jf_ahu ! index of u-diffusivity coef 71 INTEGER , SAVE :: jf_ahv ! index of v-diffusivity coef 72 INTEGER , SAVE :: jf_ahw ! index of w-diffusivity coef 73 INTEGER , SAVE :: jf_eiu ! index of u-eiv 74 INTEGER , SAVE :: jf_eiv ! index of v-eiv 75 INTEGER , SAVE :: jf_eiw ! index of w-eiv 76 INTEGER , SAVE :: jf_fmf ! index of downward salt flux 77 78 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_dyn ! structure of input fields (file informations, fields read) 77 INTEGER , SAVE :: jf_div ! index of e3t 78 79 80 TYPE(FLD), ALLOCATABLE, SAVE, DIMENSION(:) :: sf_dyn ! structure of input fields (file informations, fields read) 79 81 ! ! 80 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wdta ! vertical velocity at 2 time step81 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,: ) :: wnow ! vertical velocity at 2 time step82 82 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uslpdta ! zonal isopycnal slopes 83 83 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vslpdta ! meridional isopycnal slopes 84 84 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpidta ! zonal diapycnal slopes 85 85 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpjdta ! meridional diapycnal slopes 86 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: uslpnow ! zonal isopycnal slopes 87 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: vslpnow ! meridional isopycnal slopes 88 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wslpinow ! zonal diapycnal slopes 89 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wslpjnow ! meridional diapycnal slopes 90 91 INTEGER :: nrecprev_tem , nrecprev_uwd 86 87 INTEGER, SAVE :: nprevrec, nsecdyn 92 88 93 89 !! * Substitutions … … 113 109 !!---------------------------------------------------------------------- 114 110 ! 115 USE oce, ONLY: zts => tsa 116 USE oce, ONLY: zuslp => ua , zvslp => va 117 USE oce, ONLY: zwslpi => rotb , zwslpj => rotn 118 USE oce, ONLY: zu => ub , zv => vb, zw => hdivb 119 ! 111 USE oce, ONLY: zhdivtr => ua 120 112 INTEGER, INTENT(in) :: kt ! ocean time-step index 121 ! 122 INTEGER :: ji, jj ! dummy loop indices 123 INTEGER :: isecsbc ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step 124 REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation 125 REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation 126 INTEGER :: iswap_tem, iswap_uwd ! 113 INTEGER :: ji, jj, jk 114 REAL(wp), POINTER, DIMENSION(:,:) :: zemp 115 ! 127 116 !!---------------------------------------------------------------------- 128 117 … … 130 119 IF( nn_timing == 1 ) CALL timing_start( 'dta_dyn') 131 120 ! 132 isecsbc = nsec_year + nsec1jan000 133 ! 134 IF( kt == nit000 ) THEN 135 nrecprev_tem = 0 136 nrecprev_uwd = 0 137 ! 138 CALL fld_read( kt, 1, sf_dyn ) !== read data at kt time step ==! 139 ! 140 IF( lk_ldfslp .AND. .NOT.lk_c1d .AND. sf_dyn(jf_tem)%ln_tint ) THEN ! Computes slopes (here avt is used as workspace) 141 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:) ! temperature 142 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:) ! salinity 143 avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:) ! vertical diffusive coef. 144 CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 145 uslpdta (:,:,:,1) = zuslp (:,:,:) 146 vslpdta (:,:,:,1) = zvslp (:,:,:) 147 wslpidta(:,:,:,1) = zwslpi(:,:,:) 148 wslpjdta(:,:,:,1) = zwslpj(:,:,:) 149 ENDIF 150 IF( ln_dynwzv .AND. sf_dyn(jf_uwd)%ln_tint ) THEN ! compute vertical velocity from u/v 151 zu(:,:,:) = sf_dyn(jf_uwd)%fdta(:,:,:,1) 152 zv(:,:,:) = sf_dyn(jf_vwd)%fdta(:,:,:,1) 153 CALL dta_dyn_wzv( zu, zv, zw ) 154 wdta(:,:,:,1) = zw(:,:,:) * tmask(:,:,:) 155 ENDIF 156 ELSE 157 nrecprev_tem = sf_dyn(jf_tem)%nrec_a(2) 158 nrecprev_uwd = sf_dyn(jf_uwd)%nrec_a(2) 159 ! 160 CALL fld_read( kt, 1, sf_dyn ) !== read data at kt time step ==! 161 ! 162 ENDIF 163 ! 164 IF( lk_ldfslp .AND. .NOT.lk_c1d ) THEN ! Computes slopes (here avt is used as workspace) 165 iswap_tem = 0 166 IF( kt /= nit000 .AND. ( sf_dyn(jf_tem)%nrec_a(2) - nrecprev_tem ) /= 0 ) iswap_tem = 1 167 IF( ( isecsbc > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap_tem == 1 ) .OR. kt == nit000 ) THEN ! read/update the after data 168 IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt 169 IF( sf_dyn(jf_tem)%ln_tint ) THEN ! time interpolation of data 170 IF( kt /= nit000 ) THEN 171 uslpdta (:,:,:,1) = uslpdta (:,:,:,2) ! swap the data 172 vslpdta (:,:,:,1) = vslpdta (:,:,:,2) 173 wslpidta(:,:,:,1) = wslpidta(:,:,:,2) 174 wslpjdta(:,:,:,1) = wslpjdta(:,:,:,2) 175 ENDIF 176 ! 177 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:) ! temperature 178 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity 179 avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef. 180 CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 181 ! 182 uslpdta (:,:,:,2) = zuslp (:,:,:) 183 vslpdta (:,:,:,2) = zvslp (:,:,:) 184 wslpidta(:,:,:,2) = zwslpi(:,:,:) 185 wslpjdta(:,:,:,2) = zwslpj(:,:,:) 186 ELSE 187 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) 188 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) 189 avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) 190 CALL dta_dyn_slp( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 191 uslpnow (:,:,:) = zuslp (:,:,:) 192 vslpnow (:,:,:) = zvslp (:,:,:) 193 wslpinow(:,:,:) = zwslpi(:,:,:) 194 wslpjnow(:,:,:) = zwslpj(:,:,:) 195 ENDIF 196 ENDIF 197 IF( sf_dyn(jf_tem)%ln_tint ) THEN 198 ztinta = REAL( isecsbc - sf_dyn(jf_tem)%nrec_b(2), wp ) & 199 & / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp ) 200 ztintb = 1. - ztinta 201 uslp (:,:,:) = ztintb * uslpdta (:,:,:,1) + ztinta * uslpdta (:,:,:,2) 202 vslp (:,:,:) = ztintb * vslpdta (:,:,:,1) + ztinta * vslpdta (:,:,:,2) 203 wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1) + ztinta * wslpidta(:,:,:,2) 204 wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1) + ztinta * wslpjdta(:,:,:,2) 205 ELSE 206 uslp (:,:,:) = uslpnow (:,:,:) 207 vslp (:,:,:) = vslpnow (:,:,:) 208 wslpi(:,:,:) = wslpinow(:,:,:) 209 wslpj(:,:,:) = wslpjnow(:,:,:) 210 ENDIF 211 ENDIF 212 ! 213 IF( ln_dynwzv ) THEN ! compute vertical velocity from u/v 214 iswap_uwd = 0 215 IF( kt /= nit000 .AND. ( sf_dyn(jf_uwd)%nrec_a(2) - nrecprev_uwd ) /= 0 ) iswap_uwd = 1 216 IF( ( isecsbc > sf_dyn(jf_uwd)%nrec_b(2) .AND. iswap_uwd == 1 ) .OR. kt == nit000 ) THEN ! read/update the after data 217 IF(lwp) WRITE(numout,*) ' Compute new vertical velocity at kt = ', kt 218 IF(lwp) WRITE(numout,*) 219 IF( sf_dyn(jf_uwd)%ln_tint ) THEN ! time interpolation of data 220 IF( kt /= nit000 ) THEN 221 wdta(:,:,:,1) = wdta(:,:,:,2) ! swap the data for initialisation 222 ENDIF 223 zu(:,:,:) = sf_dyn(jf_uwd)%fdta(:,:,:,2) 224 zv(:,:,:) = sf_dyn(jf_vwd)%fdta(:,:,:,2) 225 CALL dta_dyn_wzv( zu, zv, zw ) 226 wdta(:,:,:,2) = zw(:,:,:) * tmask(:,:,:) 227 ELSE 228 zu(:,:,:) = sf_dyn(jf_uwd)%fnow(:,:,:) 229 zv(:,:,:) = sf_dyn(jf_vwd)%fnow(:,:,:) 230 CALL dta_dyn_wzv( zu, zv, zw ) 231 wnow(:,:,:) = zw(:,:,:) * tmask(:,:,:) 232 ENDIF 233 ENDIF 234 IF( sf_dyn(jf_uwd)%ln_tint ) THEN 235 ztinta = REAL( isecsbc - sf_dyn(jf_uwd)%nrec_b(2), wp ) & 236 & / REAL( sf_dyn(jf_uwd)%nrec_a(2) - sf_dyn(jf_uwd)%nrec_b(2), wp ) 237 ztintb = 1. - ztinta 238 wn(:,:,:) = ztintb * wdta(:,:,:,1) + ztinta * wdta(:,:,:,2) 239 ELSE 240 wn(:,:,:) = wnow(:,:,:) 241 ENDIF 242 ENDIF 121 ! 122 nsecdyn = nsec_year + nsec1jan000 ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step 123 ! 124 IF( kt == nit000 ) THEN ; nprevrec = 0 125 ELSE ; nprevrec = sf_dyn(jf_tem)%nrec_a(2) 126 ENDIF 127 ! 128 CALL fld_read( kt, 1, sf_dyn ) != read data at kt time step ==! 129 ! 130 IF( lk_ldfslp .AND. .NOT.lk_c1d ) CALL dta_dyn_slp( kt ) ! Computation of slopes 243 131 ! 244 132 tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature 245 133 tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity 246 ! 134 wndm(:,:) = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1) ! wind speed - needed for gas exchange 135 fmmflx(:,:) = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1) ! downward salt flux (v3.5+) 136 fr_i(:,:) = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1) ! Sea-ice fraction 137 qsr (:,:) = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1) ! solar radiation 138 emp (:,:) = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1) ! E-P 139 IF( ln_dynrnf ) THEN 140 rnf (:,:) = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1) ! E-P 141 IF( ln_dynrnf_depth .AND. lk_vvl ) CALL dta_dyn_hrnf 142 ENDIF 143 ! 144 un(:,:,:) = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:) ! effective u-transport 145 vn(:,:,:) = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:) ! effective v-transport 146 wn(:,:,:) = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:) ! effective v-transport 147 ! 148 IF( lk_vvl ) THEN 149 CALL wrk_alloc(jpi, jpj, zemp ) 150 zhdivtr(:,:,:) = sf_dyn(jf_div)%fnow(:,:,:) * tmask(:,:,:) ! effective u-transport 151 emp_b (:,:) = sf_dyn(jf_empb)%fnow(:,:,1) * tmask(:,:,1) ! E-P 152 zemp(:,:) = 0.5_wp * ( emp(:,:) + emp_b(:,:) ) + rnf(:,:) + fwbcorr * tmask(:,:,1) 153 CALL dta_dyn_ssh( kt, zhdivtr, sshb, zemp, ssha, fse3t_a(:,:,:) ) != ssh, vertical scale factor & vertical transport 154 CALL wrk_dealloc(jpi, jpj, zemp ) 155 ! Write in the tracer restart file 156 ! ******************************* 157 IF( lrst_trc ) THEN 158 IF(lwp) WRITE(numout,*) 159 IF(lwp) WRITE(numout,*) 'dta_dyn_ssh : ssh field written in tracer restart file ', & 160 & 'at it= ', kt,' date= ', ndastp 161 IF(lwp) WRITE(numout,*) '~~~~' 162 CALL iom_rstput( kt, nitrst, numrtw, 'sshn', ssha ) 163 CALL iom_rstput( kt, nitrst, numrtw, 'sshb', sshn ) 164 ENDIF 165 ENDIF 247 166 ! 248 167 CALL eos ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop … … 251 170 252 171 rn2b(:,:,:) = rn2(:,:,:) ! need for zdfmxl 253 CALL zdf_mxl( kt ) ! In any case, we need mxl 254 ! 255 avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coefficient 256 un (:,:,:) = sf_dyn(jf_uwd)%fnow(:,:,:) * umask(:,:,:) ! u-velocity 257 vn (:,:,:) = sf_dyn(jf_vwd)%fnow(:,:,:) * vmask(:,:,:) ! v-velocity 258 IF( .NOT.ln_dynwzv ) & ! w-velocity read in file 259 wn (:,:,:) = sf_dyn(jf_wwd)%fnow(:,:,:) * tmask(:,:,:) 260 hmld(:,:) = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1) ! mixed layer depht 261 wndm(:,:) = sf_dyn(jf_wnd)%fnow(:,:,1) * tmask(:,:,1) ! wind speed - needed for gas exchange 262 emp (:,:) = sf_dyn(jf_emp)%fnow(:,:,1) * tmask(:,:,1) ! E-P 263 fmmflx(:,:) = sf_dyn(jf_fmf)%fnow(:,:,1) * tmask(:,:,1) ! downward salt flux (v3.5+) 264 fr_i(:,:) = sf_dyn(jf_ice)%fnow(:,:,1) * tmask(:,:,1) ! Sea-ice fraction 265 qsr (:,:) = sf_dyn(jf_qsr)%fnow(:,:,1) * tmask(:,:,1) ! solar radiation 266 IF( ln_dynrnf ) & 267 rnf (:,:) = sf_dyn(jf_rnf)%fnow(:,:,1) * tmask(:,:,1) ! river runoffs 268 269 ! ! bbl diffusive coef 172 CALL zdf_mxl( kt ) ! In any case, we need mxl 173 ! 174 hmld(:,:) = sf_dyn(jf_mld)%fnow(:,:,1) * tmask(:,:,1) ! mixed layer depht 175 avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coefficient 176 ! 270 177 #if defined key_trabbl && ! defined key_c1d 271 IF( ln_dynbbl ) THEN ! read in a file 272 ahu_bbl(:,:) = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1) 273 ahv_bbl(:,:) = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1) 274 ELSE ! Compute bbl coefficients if needed 275 tsb(:,:,:,:) = tsn(:,:,:,:) 276 CALL bbl( kt, nit000, 'TRC') 277 END IF 178 ahu_bbl(:,:) = sf_dyn(jf_ubl)%fnow(:,:,1) * umask(:,:,1) ! bbl diffusive coef 179 ahv_bbl(:,:) = sf_dyn(jf_vbl)%fnow(:,:,1) * vmask(:,:,1) 278 180 #endif 279 #if ( ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv ) && ! defined key_c1d 280 aeiw(:,:) = sf_dyn(jf_eiw)%fnow(:,:,1) * tmask(:,:,1) ! w-eiv 281 ! ! Computes the horizontal values from the vertical value 282 DO jj = 2, jpjm1 283 DO ji = fs_2, fs_jpim1 ! vector opt. 284 aeiu(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji+1,jj ) ) ! Average the diffusive coefficient at u- v- points 285 aeiv(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji ,jj+1) ) ! at u- v- points 286 END DO 287 END DO 288 CALL lbc_lnk( aeiu, 'U', 1. ) ; CALL lbc_lnk( aeiv, 'V', 1. ) ! lateral boundary condition 289 #endif 290 291 #if defined key_degrad && ! defined key_c1d 292 ! ! degrad option : diffusive and eiv coef are 3D 293 ahtu(:,:,:) = sf_dyn(jf_ahu)%fnow(:,:,:) * umask(:,:,:) 294 ahtv(:,:,:) = sf_dyn(jf_ahv)%fnow(:,:,:) * vmask(:,:,:) 295 ahtw(:,:,:) = sf_dyn(jf_ahw)%fnow(:,:,:) * tmask(:,:,:) 296 # if defined key_traldf_eiv 297 aeiu(:,:,:) = sf_dyn(jf_eiu)%fnow(:,:,:) * umask(:,:,:) 298 aeiv(:,:,:) = sf_dyn(jf_eiv)%fnow(:,:,:) * vmask(:,:,:) 299 aeiw(:,:,:) = sf_dyn(jf_eiw)%fnow(:,:,:) * tmask(:,:,:) 300 # endif 301 #endif 181 ! 182 ! 183 CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop 302 184 ! 303 185 IF(ln_ctl) THEN ! print control … … 308 190 CALL prt_ctl(tab3d_1=wn , clinfo1=' wn - : ', mask1=tmask, ovlap=1, kdim=jpk ) 309 191 CALL prt_ctl(tab3d_1=avt , clinfo1=' kz - : ', mask1=tmask, ovlap=1, kdim=jpk ) 310 CALL prt_ctl(tab2d_1=fr_i , clinfo1=' fr_i - : ', mask1=tmask, ovlap=1 )311 CALL prt_ctl(tab2d_1=hmld , clinfo1=' hmld - : ', mask1=tmask, ovlap=1 )312 CALL prt_ctl(tab2d_1=fmmflx , clinfo1=' fmmflx - : ', mask1=tmask, ovlap=1 )313 CALL prt_ctl(tab2d_1=emp , clinfo1=' emp - : ', mask1=tmask, ovlap=1 )314 CALL prt_ctl(tab2d_1=wndm , clinfo1=' wspd - : ', mask1=tmask, ovlap=1 )315 CALL prt_ctl(tab2d_1=qsr , clinfo1=' qsr - : ', mask1=tmask, ovlap=1 )192 ! CALL prt_ctl(tab2d_1=fr_i , clinfo1=' fr_i - : ', mask1=tmask, ovlap=1 ) 193 ! CALL prt_ctl(tab2d_1=hmld , clinfo1=' hmld - : ', mask1=tmask, ovlap=1 ) 194 ! CALL prt_ctl(tab2d_1=fmmflx , clinfo1=' fmmflx - : ', mask1=tmask, ovlap=1 ) 195 ! CALL prt_ctl(tab2d_1=emp , clinfo1=' emp - : ', mask1=tmask, ovlap=1 ) 196 ! CALL prt_ctl(tab2d_1=wndm , clinfo1=' wspd - : ', mask1=tmask, ovlap=1 ) 197 ! CALL prt_ctl(tab2d_1=qsr , clinfo1=' qsr - : ', mask1=tmask, ovlap=1 ) 316 198 ENDIF 317 199 ! … … 335 217 INTEGER :: inum, idv, idimv ! local integer 336 218 INTEGER :: ios ! Local integer output status for namelist read 337 !! 338 CHARACTER(len=100) :: cn_dir ! Root directory for location of core files 339 TYPE(FLD_N), DIMENSION(jpfld) :: slf_d ! array of namelist informations on the fields to read 340 TYPE(FLD_N) :: sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd, sn_rnf ! informations about the fields to be read 341 TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl ! " " 342 TYPE(FLD_N) :: sn_ahu, sn_ahv, sn_ahw, sn_eiu, sn_eiv, sn_eiw, sn_fmf ! " " 343 !!---------------------------------------------------------------------- 344 ! 345 NAMELIST/namdta_dyn/cn_dir, ln_dynwzv, ln_dynbbl, ln_degrad, ln_dynrnf, & 346 & sn_tem, sn_sal, sn_mld, sn_emp, sn_ice, sn_qsr, sn_wnd, sn_rnf, & 347 & sn_uwd, sn_vwd, sn_wwd, sn_avt, sn_ubl, sn_vbl, & 348 & sn_ahu, sn_ahv, sn_ahw, sn_eiu, sn_eiv, sn_eiw, sn_fmf 219 INTEGER :: ji, jj, jk 220 REAL(wp) :: zcoef 221 INTEGER :: nkrnf_max 222 REAL(wp) :: hrnf_max 223 !! 224 CHARACTER(len=100) :: cn_dir ! Root directory for location of core files 225 TYPE(FLD_N), DIMENSION(jpfld) :: slf_d ! array of namelist informations on the fields to read 226 TYPE(FLD_N) :: sn_uwd, sn_vwd, sn_wwd, sn_empb, sn_emp ! informations about the fields to be read 227 TYPE(FLD_N) :: sn_tem , sn_sal , sn_avt ! " " 228 TYPE(FLD_N) :: sn_mld, sn_qsr, sn_wnd , sn_ice , sn_fmf ! " " 229 TYPE(FLD_N) :: sn_ubl, sn_vbl, sn_rnf ! " " 230 TYPE(FLD_N) :: sn_div ! informations about the fields to be read 231 !!---------------------------------------------------------------------- 232 233 NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth, fwbcorr, & 234 & sn_uwd, sn_vwd, sn_wwd, sn_emp, & 235 & sn_avt, sn_tem, sn_sal, sn_mld , sn_qsr , & 236 & sn_wnd, sn_ice, sn_fmf, & 237 & sn_ubl, sn_vbl, sn_rnf, & 238 & sn_empb, sn_div 349 239 ! 350 240 REWIND( numnam_ref ) ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data … … 363 253 WRITE(numout,*) '~~~~~~~ ' 364 254 WRITE(numout,*) ' Namelist namdta_dyn' 365 WRITE(numout,*) ' vertical velocity read from file (T) or computed (F) ln_dynwzv = ', ln_dynwzv 366 WRITE(numout,*) ' bbl coef read from file (T) or computed (F) ln_dynbbl = ', ln_dynbbl 367 WRITE(numout,*) ' degradation option enabled (T) or not (F) ln_degrad = ', ln_degrad 368 WRITE(numout,*) ' river runoff option enabled (T) or not (F) ln_dynrnf = ', ln_dynrnf 255 WRITE(numout,*) ' runoffs option enabled (T) or not (F) ln_dynrnf = ', ln_dynrnf 256 WRITE(numout,*) ' runoffs is spread in vertical ln_dynrnf_depth = ', ln_dynrnf_depth 257 WRITE(numout,*) ' annual global mean of empmr for ssh correction fwbcorr = ', fwbcorr 369 258 WRITE(numout,*) 370 259 ENDIF 371 260 ! 372 IF( ln_degrad .AND. .NOT.lk_degrad ) THEN 373 CALL ctl_warn( 'dta_dyn_init: degradation option requires key_degrad activated ; force ln_degrad to false' ) 374 ln_degrad = .FALSE. 375 ENDIF 376 IF( ln_dynbbl .AND. ( .NOT.lk_trabbl .OR. lk_c1d ) ) THEN 377 CALL ctl_warn( 'dta_dyn_init: bbl option requires key_trabbl activated ; force ln_dynbbl to false' ) 378 ln_dynbbl = .FALSE. 379 ENDIF 380 381 jf_tem = 1 ; jf_sal = 2 ; jf_mld = 3 ; jf_emp = 4 ; jf_fmf = 5 ; jf_ice = 6 ; jf_qsr = 7 382 jf_wnd = 8 ; jf_uwd = 9 ; jf_vwd = 10 ; jf_wwd = 11 ; jf_avt = 12 ; jfld = jf_avt 383 ! 384 slf_d(jf_tem) = sn_tem ; slf_d(jf_sal) = sn_sal ; slf_d(jf_mld) = sn_mld 385 slf_d(jf_emp) = sn_emp ; slf_d(jf_fmf ) = sn_fmf ; slf_d(jf_ice) = sn_ice 386 slf_d(jf_qsr) = sn_qsr ; slf_d(jf_wnd) = sn_wnd ; slf_d(jf_avt) = sn_avt 387 slf_d(jf_uwd) = sn_uwd ; slf_d(jf_vwd) = sn_vwd ; slf_d(jf_wwd) = sn_wwd 388 261 262 jf_uwd = 1 ; jf_vwd = 2 ; jf_wwd = 3 ; jf_emp = 4 ; jf_avt = 5 263 jf_tem = 6 ; jf_sal = 7 ; jf_mld = 8 ; jf_qsr = 9 264 jf_wnd = 10 ; jf_ice = 11 ; jf_fmf = 12 ; jfld = jf_fmf 265 266 ! 267 slf_d(jf_uwd) = sn_uwd ; slf_d(jf_vwd) = sn_vwd ; slf_d(jf_wwd) = sn_wwd 268 slf_d(jf_emp) = sn_emp ; slf_d(jf_avt) = sn_avt 269 slf_d(jf_tem) = sn_tem ; slf_d(jf_sal) = sn_sal ; slf_d(jf_mld) = sn_mld 270 slf_d(jf_qsr) = sn_qsr ; slf_d(jf_wnd) = sn_wnd ; slf_d(jf_ice) = sn_ice 271 slf_d(jf_fmf) = sn_fmf 272 273 274 ! 275 IF( lk_vvl ) THEN 276 jf_div = jfld + 1 ; jf_empb = jfld + 2 ; jfld = jf_empb 277 slf_d(jf_div) = sn_div ; slf_d(jf_empb) = sn_empb 278 ENDIF 279 ! 280 IF( lk_trabbl ) THEN 281 jf_ubl = jfld + 1 ; jf_vbl = jfld + 2 ; jfld = jf_vbl 282 slf_d(jf_ubl) = sn_ubl ; slf_d(jf_vbl) = sn_vbl 283 ENDIF 389 284 ! 390 285 IF( ln_dynrnf ) THEN 391 jf_rnf = jfld + 1 ;jfld = jf_rnf392 slf_d(jf_rnf) = sn_rnf286 jf_rnf = jfld + 1 ; jfld = jf_rnf 287 slf_d(jf_rnf) = sn_rnf 393 288 ELSE 394 rnf (:,:) = 0._wp 395 ENDIF 396 397 ! 398 IF( .NOT.ln_degrad ) THEN ! no degrad option 399 IF( lk_traldf_eiv .AND. ln_dynbbl ) THEN ! eiv & bbl 400 jf_ubl = jfld + 1 ; jf_vbl = jfld + 2 ; jf_eiw = jfld + 3 ; jfld = jf_eiw 401 slf_d(jf_ubl) = sn_ubl ; slf_d(jf_vbl) = sn_vbl ; slf_d(jf_eiw) = sn_eiw 402 ENDIF 403 IF( .NOT.lk_traldf_eiv .AND. ln_dynbbl ) THEN ! no eiv & bbl 404 jf_ubl = jfld + 1 ; jf_vbl = jfld + 2 ; jfld = jf_vbl 405 slf_d(jf_ubl) = sn_ubl ; slf_d(jf_vbl) = sn_vbl 406 ENDIF 407 IF( lk_traldf_eiv .AND. .NOT.ln_dynbbl ) THEN ! eiv & no bbl 408 jf_eiw = jfld + 1 ; jfld = jf_eiw ; slf_d(jf_eiw) = sn_eiw 409 ENDIF 410 ELSE 411 jf_ahu = jfld + 1 ; jf_ahv = jfld + 2 ; jf_ahw = jfld + 3 ; jfld = jf_ahw 412 slf_d(jf_ahu) = sn_ahu ; slf_d(jf_ahv) = sn_ahv ; slf_d(jf_ahw) = sn_ahw 413 IF( lk_traldf_eiv .AND. ln_dynbbl ) THEN ! eiv & bbl 414 jf_ubl = jfld + 1 ; jf_vbl = jfld + 2 ; 415 slf_d(jf_ubl) = sn_ubl ; slf_d(jf_vbl) = sn_vbl 416 jf_eiu = jfld + 3 ; jf_eiv = jfld + 4 ; jf_eiw = jfld + 5 ; jfld = jf_eiw 417 slf_d(jf_eiu) = sn_eiu ; slf_d(jf_eiv) = sn_eiv ; slf_d(jf_eiw) = sn_eiw 418 ENDIF 419 IF( .NOT.lk_traldf_eiv .AND. ln_dynbbl ) THEN ! no eiv & bbl 420 jf_ubl = jfld + 1 ; jf_vbl = jfld + 2 ; jfld = jf_vbl 421 slf_d(jf_ubl) = sn_ubl ; slf_d(jf_vbl) = sn_vbl 422 ENDIF 423 IF( lk_traldf_eiv .AND. .NOT.ln_dynbbl ) THEN ! eiv & no bbl 424 jf_eiu = jfld + 1 ; jf_eiv = jfld + 2 ; jf_eiw = jfld + 3 ; jfld = jf_eiw 425 slf_d(jf_eiu) = sn_eiu ; slf_d(jf_eiv) = sn_eiv ; slf_d(jf_eiw) = sn_eiw 426 ENDIF 427 ENDIF 289 rnf(:,:) = 0._wp 290 ENDIF 291 428 292 429 293 ALLOCATE( sf_dyn(jfld), STAT=ierr ) ! set sf structure 430 IF( ierr > 0 ) THEN294 IF( ierr > 0 ) THEN 431 295 CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' ) ; RETURN 432 296 ENDIF 433 297 ! ! fill sf with slf_i and control print 434 298 CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' ) 299 ! 435 300 ! Open file for each variable to get his number of dimension 436 301 DO ifpr = 1, jfld … … 456 321 ALLOCATE( uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2), & 457 322 & wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2), STAT=ierr2 ) 458 ELSE 459 ALLOCATE( uslpnow (jpi,jpj,jpk) , vslpnow (jpi,jpj,jpk) , & 460 & wslpinow(jpi,jpj,jpk) , wslpjnow(jpi,jpj,jpk) , STAT=ierr2 ) 461 ENDIF 462 IF( ierr2 > 0 ) THEN 463 CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' ) ; RETURN 323 ! 324 IF( ierr2 > 0 ) THEN 325 CALL ctl_stop( 'dta_dyn_init : unable to allocate slope arrays' ) ; RETURN 326 ENDIF 464 327 ENDIF 465 328 ENDIF 466 IF( ln_dynwzv ) THEN ! slopes 467 IF( sf_dyn(jf_uwd)%ln_tint ) THEN ! time interpolation 468 ALLOCATE( wdta(jpi,jpj,jpk,2), STAT=ierr3 ) 469 ELSE 470 ALLOCATE( wnow(jpi,jpj,jpk) , STAT=ierr3 ) 471 ENDIF 472 IF( ierr3 > 0 ) THEN 473 CALL ctl_stop( 'dta_dyn_init : unable to allocate wdta arrays' ) ; RETURN 474 ENDIF 475 ENDIF 476 ! 477 CALL dta_dyn( nit000 ) 478 ! 479 END SUBROUTINE dta_dyn_init 480 481 SUBROUTINE dta_dyn_wzv( pu, pv, pw ) 482 !!---------------------------------------------------------------------- 483 !! *** ROUTINE wzv *** 484 !! 485 !! ** Purpose : Compute the now vertical velocity after the array swap 486 !! 487 !! ** Method : - compute the now divergence given by : 488 !! * z-coordinate ONLY !!!! 489 !! hdiv = 1/(e1t*e2t) [ di(e2u u) + dj(e1v v) ] 490 !! - Using the incompressibility hypothesis, the vertical 491 !! velocity is computed by integrating the horizontal divergence 492 !! from the bottom to the surface. 493 !! The boundary conditions are w=0 at the bottom (no flux). 494 !!---------------------------------------------------------------------- 495 USE oce, ONLY: zhdiv => hdivn 496 ! 497 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu, pv !: horizontal velocities 498 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: pw !: vertical velocity 499 !! 500 INTEGER :: ji, jj, jk 501 REAL(wp) :: zu, zu1, zv, zv1, zet 502 !!---------------------------------------------------------------------- 503 ! 504 ! Computation of vertical velocity using horizontal divergence 505 zhdiv(:,:,:) = 0._wp 506 DO jk = 1, jpkm1 507 DO jj = 2, jpjm1 508 DO ji = fs_2, fs_jpim1 ! vector opt. 509 zu = pu(ji ,jj ,jk) * umask(ji ,jj ,jk) * e2u(ji ,jj ) * fse3u(ji ,jj ,jk) 510 zu1 = pu(ji-1,jj ,jk) * umask(ji-1,jj ,jk) * e2u(ji-1,jj ) * fse3u(ji-1,jj ,jk) 511 zv = pv(ji ,jj ,jk) * vmask(ji ,jj ,jk) * e1v(ji ,jj ) * fse3v(ji ,jj ,jk) 512 zv1 = pv(ji ,jj-1,jk) * vmask(ji ,jj-1,jk) * e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) 513 zet = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 514 zhdiv(ji,jj,jk) = ( zu - zu1 + zv - zv1 ) * zet 329 ! 330 IF( lk_vvl ) THEN 331 IF( .NOT. sf_dyn(jf_uwd)%ln_clim .AND. ln_rsttr .AND. & ! Restart: read in restart file 332 iom_varid( numrtr, 'sshn', ldstop = .FALSE. ) > 0 ) THEN 333 IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation' 334 CALL iom_get( numrtr, jpdom_autoglo, 'sshn', sshn(:,:) ) 335 CALL iom_get( numrtr, jpdom_autoglo, 'sshb', sshb(:,:) ) 336 ELSE 337 IF(lwp) WRITE(numout,*) ' sshn forcing fields read in the restart file for initialisation' 338 CALL iom_open( 'restart', inum ) 339 CALL iom_get( inum, jpdom_autoglo, 'sshn', sshn(:,:) ) 340 CALL iom_get( inum, jpdom_autoglo, 'sshb', sshb(:,:) ) 341 CALL iom_close( inum ) ! close file 342 ENDIF 343 ! 344 DO jk = 1, jpkm1 345 fse3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + sshn(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 346 ENDDO 347 fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 348 349 ! Horizontal scale factor interpolations 350 ! -------------------------------------- 351 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 352 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 353 354 ! Vertical scale factor interpolations 355 ! ------------------------------------ 356 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n(:,:,:), 'W' ) 357 358 fse3t_b(:,:,:) = fse3t_n(:,:,:) 359 fse3u_b(:,:,:) = fse3u_n(:,:,:) 360 fse3v_b(:,:,:) = fse3v_n(:,:,:) 361 362 ! t- and w- points depth 363 ! ---------------------- 364 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 365 fsdepw_n(:,:,1) = 0.0_wp 366 367 DO jk = 2, jpk 368 DO jj = 1,jpj 369 DO ji = 1,jpi 370 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere 371 ! tmask = wmask, ie everywhere expect at jk = mikt 372 ! 1 for jk = 373 ! mikt 374 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 375 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 376 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 377 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 378 END DO 379 END DO 380 END DO 381 382 fsdept_b(:,:,:) = fsdept_n(:,:,:) 383 fsdepw_b(:,:,:) = fsdepw_n(:,:,:) 384 ! 385 ENDIF 386 ! 387 IF( ln_dynrnf .AND. ln_dynrnf_depth ) THEN ! read depht over which runoffs are distributed 388 IF(lwp) WRITE(numout,*) 389 IF(lwp) WRITE(numout,*) ' read in the file depht over which runoffs are distributed' 390 CALL iom_open ( "runoffs", inum ) ! open file 391 CALL iom_get ( inum, jpdom_data, 'rodepth', h_rnf ) ! read the river mouth array 392 CALL iom_close( inum ) ! close file 393 ! 394 nk_rnf(:,:) = 0 ! set the number of level over which river runoffs are applied 395 DO jj = 1, jpj 396 DO ji = 1, jpi 397 IF( h_rnf(ji,jj) > 0._wp ) THEN 398 jk = 2 399 DO WHILE ( jk /= mbkt(ji,jj) .AND. gdept_0(ji,jj,jk) < h_rnf(ji,jj) ) ; jk = jk + 1 400 END DO 401 nk_rnf(ji,jj) = jk 402 ELSEIF( h_rnf(ji,jj) == -1._wp ) THEN ; nk_rnf(ji,jj) = 1 403 ELSEIF( h_rnf(ji,jj) == -999._wp ) THEN ; nk_rnf(ji,jj) = mbkt(ji,jj) 404 ELSE 405 CALL ctl_stop( 'sbc_rnf_init: runoff depth not positive, and not -999 or -1, rnf value in file fort.999' ) 406 WRITE(999,*) 'ji, jj, h_rnf(ji,jj) :', ji, jj, h_rnf(ji,jj) 407 ENDIF 515 408 END DO 516 409 END DO 410 DO jj = 1, jpj ! set the associated depth 411 DO ji = 1, jpi 412 h_rnf(ji,jj) = 0._wp 413 DO jk = 1, nk_rnf(ji,jj) 414 h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk) 415 END DO 416 END DO 417 END DO 418 ELSE ! runoffs applied at the surface 419 nk_rnf(:,:) = 1 420 h_rnf (:,:) = fse3t(:,:,1) 421 ENDIF 422 nkrnf_max = MAXVAL( nk_rnf(:,:) ) 423 hrnf_max = MAXVAL( h_rnf(:,:) ) 424 IF( lk_mpp ) THEN 425 CALL mpp_max( nkrnf_max ) ! max over the global domain 426 CALL mpp_max( hrnf_max ) ! max over the global domain 427 ENDIF 428 IF(lwp) WRITE(numout,*) ' ' 429 IF(lwp) WRITE(numout,*) ' max depht of runoff : ', hrnf_max,' max level : ', nkrnf_max 430 IF(lwp) WRITE(numout,*) ' ' 431 ! 432 CALL dta_dyn( nit000 ) 433 ! 434 END SUBROUTINE dta_dyn_init 435 436 SUBROUTINE dta_dyn_swp( kt ) 437 !!--------------------------------------------------------------------- 438 !! *** ROUTINE dta_dyn_swp *** 439 !! 440 !! ** Purpose : Swap and the data and compute the vertical scale factor at U/V/W point 441 !! and the depht 442 !! 443 !!--------------------------------------------------------------------- 444 INTEGER, INTENT(in) :: kt ! time step 445 INTEGER :: ji, jj, jk 446 REAL(wp) :: zcoef 447 ! 448 !!--------------------------------------------------------------------- 449 450 IF( kt == nit000 ) THEN 451 IF(lwp) WRITE(numout,*) 452 IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height' 453 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 454 ENDIF 455 456 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:)) ! before <-- now filtered 457 sshn(:,:) = ssha(:,:) 458 459 fse3t_n(:,:,:) = fse3t_a(:,:,:) 460 461 ! Reconstruction of all vertical scale factors at now and before time steps 462 ! ============================================================================= 463 464 ! Horizontal scale factor interpolations 465 ! -------------------------------------- 466 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 467 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 468 469 ! Vertical scale factor interpolations 470 ! ------------------------------------ 471 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 472 473 fse3t_b(:,:,:) = fse3t_n(:,:,:) 474 fse3u_b(:,:,:) = fse3u_n(:,:,:) 475 fse3v_b(:,:,:) = fse3v_n(:,:,:) 476 477 ! t- and w- points depth 478 ! ---------------------- 479 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 480 fsdepw_n(:,:,1) = 0.0_wp 481 482 DO jk = 2, jpk 483 DO jj = 1,jpj 484 DO ji = 1,jpi 485 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 486 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 487 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 488 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 489 END DO 490 END DO 491 END DO 492 493 fsdept_b(:,:,:) = fsdept_n(:,:,:) 494 fsdepw_b(:,:,:) = fsdepw_n(:,:,:) 495 496 ! 497 END SUBROUTINE dta_dyn_swp 498 499 SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb, pemp, pssha, pe3ta ) 500 !!---------------------------------------------------------------------- 501 !! *** ROUTINE dta_dyn_wzv *** 502 !! 503 !! ** Purpose : compute the after ssh (ssha) and the now vertical velocity 504 !! 505 !! ** Method : Using the incompressibility hypothesis, 506 !! - the ssh increment is computed by integrating the horizontal divergence 507 !! and multiply by the time step. 508 !! 509 !! - compute the after scale factor : repartition of ssh INCREMENT proportionnaly 510 !! to the level thickness ( z-star case ) 511 !! 512 !! - the vertical velocity is computed by integrating the horizontal divergence 513 !! from the bottom to the surface minus the scale factor evolution. 514 !! The boundary conditions are w=0 at the bottom (no flux) 515 !! 516 !! ** action : ssha / e3t_a / wn 517 !! 518 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 519 !!---------------------------------------------------------------------- 520 !! * Arguments 521 INTEGER, INTENT(in ) :: kt ! time-step 522 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(in ) :: phdivtr ! horizontal divergence transport 523 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: psshb ! now ssh 524 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(in ) :: pemp ! evaporation minus precipitation 525 REAL(wp), DIMENSION(jpi,jpj) , OPTIONAL, INTENT(inout) :: pssha ! after ssh 526 REAL(wp), DIMENSION(jpi,jpj,jpk), OPTIONAL, INTENT(out) :: pe3ta ! after vertical scale factor 527 !! * Local declarations 528 INTEGER :: jk 529 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv 530 REAL(wp) :: z2dt 531 !!---------------------------------------------------------------------- 532 533 ! 534 z2dt = 2._wp * rdt 535 ! 536 zhdiv(:,:) = 0._wp 537 DO jk = 1, jpkm1 538 zhdiv(:,:) = zhdiv(:,:) + phdivtr(:,:,jk) * tmask(:,:,jk) 517 539 END DO 518 CALL lbc_lnk( zhdiv, 'T', 1. ) ! Lateral boundary conditions on zhdiv519 !520 ! computation of vertical velocity from the bottom521 pw(:,:,jpk) = 0._wp522 DO jk = jpkm1, 1, -1523 pw(:,:,jk) = pw(:,:,jk+1) - fse3t(:,:,jk) * zhdiv(:,:,jk)540 ! ! Sea surface elevation time-stepping 541 pssha(:,:) = ( psshb(:,:) - z2dt * ( r1_rau0 * pemp(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 542 ! ! 543 ! ! After acale factors at t-points ( z_star coordinate ) 544 DO jk = 1, jpkm1 545 pe3ta(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + pssha(:,:) * tmask(:,:,1) / ( ht_0(:,:) + 1.0 - tmask(:,:,1) ) ) 524 546 END DO 525 547 ! 526 END SUBROUTINE dta_dyn_wzv 527 528 SUBROUTINE dta_dyn_slp( kt, pts, puslp, pvslp, pwslpi, pwslpj ) 548 END SUBROUTINE dta_dyn_ssh 549 550 551 SUBROUTINE dta_dyn_hrnf 552 !!---------------------------------------------------------------------- 553 !! *** ROUTINE sbc_rnf *** 554 !! 555 !! ** Purpose : update the horizontal divergence with the runoff inflow 556 !! 557 !! ** Method : 558 !! CAUTION : rnf is positive (inflow) decreasing the 559 !! divergence and expressed in m/s 560 !! 561 !! ** Action : phdivn decreased by the runoff inflow 562 !!---------------------------------------------------------------------- 563 !! 564 INTEGER :: ji, jj, jk ! dummy loop indices 565 !!---------------------------------------------------------------------- 566 ! 567 DO jj = 1, jpj ! update the depth over which runoffs are distributed 568 DO ji = 1, jpi 569 h_rnf(ji,jj) = 0._wp 570 DO jk = 1, nk_rnf(ji,jj) ! recalculates h_rnf to be the depth in metres 571 h_rnf(ji,jj) = h_rnf(ji,jj) + fse3t(ji,jj,jk) ! to the bottom of the relevant grid box 572 END DO 573 END DO 574 END DO 575 ! 576 END SUBROUTINE dta_dyn_hrnf 577 578 579 580 SUBROUTINE dta_dyn_slp( kt ) 581 !!--------------------------------------------------------------------- 582 !! *** ROUTINE dta_dyn_slp *** 583 !! 584 !! ** Purpose : Computation of slope 585 !! 586 !!--------------------------------------------------------------------- 587 USE oce, ONLY: zts => tsa 588 ! 589 INTEGER, INTENT(in) :: kt ! time step 590 ! 591 INTEGER :: ji, jj ! dummy loop indices 592 REAL(wp) :: ztinta ! ratio applied to after records when doing time interpolation 593 REAL(wp) :: ztintb ! ratio applied to before records when doing time interpolation 594 INTEGER :: iswap 595 REAL(wp), POINTER, DIMENSION(:,:,:) :: zuslp, zvslp, zwslpi, zwslpj 596 !!--------------------------------------------------------------------- 597 ! 598 CALL wrk_alloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj ) 599 ! 600 IF( sf_dyn(jf_tem)%ln_tint ) THEN ! Computes slopes (here avt is used as workspace) 601 IF( kt == nit000 ) THEN 602 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,1) * tmask(:,:,:) ! temperature 603 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,1) * tmask(:,:,:) ! salinity 604 avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,1) * tmask(:,:,:) ! vertical diffusive coef. 605 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 606 uslpdta (:,:,:,1) = zuslp (:,:,:) 607 vslpdta (:,:,:,1) = zvslp (:,:,:) 608 wslpidta(:,:,:,1) = zwslpi(:,:,:) 609 wslpjdta(:,:,:,1) = zwslpj(:,:,:) 610 ! 611 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:) ! temperature 612 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity 613 avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef. 614 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 615 uslpdta (:,:,:,2) = zuslp (:,:,:) 616 vslpdta (:,:,:,2) = zvslp (:,:,:) 617 wslpidta(:,:,:,2) = zwslpi(:,:,:) 618 wslpjdta(:,:,:,2) = zwslpj(:,:,:) 619 ELSE 620 ! 621 iswap = 0 622 IF( sf_dyn(jf_tem)%nrec_a(2) - nprevrec /= 0 ) iswap = 1 623 IF( nsecdyn > sf_dyn(jf_tem)%nrec_b(2) .AND. iswap == 1 ) THEN ! read/update the after data 624 IF(lwp) WRITE(numout,*) ' Compute new slopes at kt = ', kt 625 uslpdta (:,:,:,1) = uslpdta (:,:,:,2) ! swap the data 626 vslpdta (:,:,:,1) = vslpdta (:,:,:,2) 627 wslpidta(:,:,:,1) = wslpidta(:,:,:,2) 628 wslpjdta(:,:,:,1) = wslpjdta(:,:,:,2) 629 ! 630 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fdta(:,:,:,2) * tmask(:,:,:) ! temperature 631 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fdta(:,:,:,2) * tmask(:,:,:) ! salinity 632 avt(:,:,:) = sf_dyn(jf_avt)%fdta(:,:,:,2) * tmask(:,:,:) ! vertical diffusive coef. 633 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 634 ! 635 uslpdta (:,:,:,2) = zuslp (:,:,:) 636 vslpdta (:,:,:,2) = zvslp (:,:,:) 637 wslpidta(:,:,:,2) = zwslpi(:,:,:) 638 wslpjdta(:,:,:,2) = zwslpj(:,:,:) 639 ENDIF 640 ENDIF 641 ENDIF 642 ! 643 IF( sf_dyn(jf_tem)%ln_tint ) THEN 644 ztinta = REAL( nsecdyn - sf_dyn(jf_tem)%nrec_b(2), wp ) & 645 & / REAL( sf_dyn(jf_tem)%nrec_a(2) - sf_dyn(jf_tem)%nrec_b(2), wp ) 646 ztintb = 1. - ztinta 647 #if defined key_ldfslp && ! defined key_c1d 648 uslp (:,:,:) = ztintb * uslpdta (:,:,:,1) + ztinta * uslpdta (:,:,:,2) 649 vslp (:,:,:) = ztintb * vslpdta (:,:,:,1) + ztinta * vslpdta (:,:,:,2) 650 wslpi(:,:,:) = ztintb * wslpidta(:,:,:,1) + ztinta * wslpidta(:,:,:,2) 651 wslpj(:,:,:) = ztintb * wslpjdta(:,:,:,1) + ztinta * wslpjdta(:,:,:,2) 652 #endif 653 ELSE 654 zts(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:) * tmask(:,:,:) ! temperature 655 zts(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:) * tmask(:,:,:) ! salinity 656 avt(:,:,:) = sf_dyn(jf_avt)%fnow(:,:,:) * tmask(:,:,:) ! vertical diffusive coef. 657 CALL compute_slopes( kt, zts, zuslp, zvslp, zwslpi, zwslpj ) 658 ! 659 #if defined key_ldfslp && ! defined key_c1d 660 uslp (:,:,:) = zuslp (:,:,:) 661 vslp (:,:,:) = zvslp (:,:,:) 662 wslpi(:,:,:) = zwslpi(:,:,:) 663 wslpj(:,:,:) = zwslpj(:,:,:) 664 #endif 665 ENDIF 666 ! 667 CALL wrk_dealloc(jpi, jpj, jpk, zuslp, zvslp, zwslpi, zwslpj ) 668 ! 669 END SUBROUTINE dta_dyn_slp 670 671 SUBROUTINE compute_slopes( kt, pts, puslp, pvslp, pwslpi, pwslpj ) 529 672 !!--------------------------------------------------------------------- 530 673 !! *** ROUTINE dta_dyn_slp *** … … 568 711 #endif 569 712 ! 570 END SUBROUTINE dta_dyn_slp 713 END SUBROUTINE compute_slopes 714 571 715 !!====================================================================== 572 716 END MODULE dtadyn
Note: See TracChangeset
for help on using the changeset viewer.
