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