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