Changeset 6579 for branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM
- Timestamp:
- 2016-05-19T16:12:37+02:00 (8 years ago)
- Location:
- branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 5 edited
-
DOM/domhgr.F90 (modified) (3 diffs)
-
DOM/domwri.F90 (modified) (7 diffs)
-
DOM/dtatsd.F90 (modified) (3 diffs)
-
DOM/istate.F90 (modified) (8 diffs)
-
usrdef.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domhgr.F90
r6484 r6579 290 290 ! 291 291 IF( ie1e2u_v == 0 ) THEN ! e1e2u and e1e2v have not been read: compute them 292 ! ! e2u and e1v does not include a 293 ! reduction in some strait: apply reduction 292 ! ! e2u and e1v does not include a reduction in some strait: apply reduction 294 293 e1e2u (:,:) = e1u(:,:) * e2u(:,:) 295 294 e1e2v (:,:) = e1v(:,:) * e2v(:,:) … … 301 300 & gphit , gphiu , gphiv , gphif , & 302 301 & e1t , e1u , e1v , e1f , & 303 & e2t , e2u , e2v , e2f , & 302 & e2t , e2u , e2v , e2f , & 304 303 & ie1e2u_v ) 305 304 ! … … 331 330 e2_e1u(:,:) = e2u(:,:) / e1u(:,:) 332 331 e1_e2v(:,:) = e1v(:,:) / e2v(:,:) 333 334 IF( lwp .AND. nn_print >=1 .AND. .NOT.ln_rstart ) THEN ! Control print : Grid informations (if not restart)335 WRITE(numout,*)336 WRITE(numout,*) ' longitude and e1 scale factors'337 WRITE(numout,*) ' ------------------------------'338 WRITE(numout,9300) ( ji, glamt(ji,1), glamu(ji,1), &339 glamv(ji,1), glamf(ji,1), &340 e1t(ji,1), e1u(ji,1), &341 e1v(ji,1), e1f(ji,1), ji = 1, jpi,10)342 9300 FORMAT( 1x, i4, f8.2,1x, f8.2,1x, f8.2,1x, f8.2, 1x, &343 f19.10, 1x, f19.10, 1x, f19.10, 1x, f19.10 )344 !345 WRITE(numout,*)346 WRITE(numout,*) ' latitude and e2 scale factors'347 WRITE(numout,*) ' -----------------------------'348 WRITE(numout,9300) ( jj, gphit(1,jj), gphiu(1,jj), &349 & gphiv(1,jj), gphif(1,jj), &350 & e2t (1,jj), e2u (1,jj), &351 & e2v (1,jj), e2f (1,jj), jj = 1, jpj, 10 )352 ENDIF353 332 354 333 -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
r5836 r6579 71 71 ! 72 72 ! ! horizontal mesh (inum3) 73 CALL iom_rstput( 0, 0, inum0, 'glamt', glamt, ktype = jp_r 4) ! ! latitude74 CALL iom_rstput( 0, 0, inum0, 'glamu', glamu, ktype = jp_r 4)75 CALL iom_rstput( 0, 0, inum0, 'glamv', glamv, ktype = jp_r 4)76 CALL iom_rstput( 0, 0, inum0, 'glamf', glamf, ktype = jp_r 4)77 78 CALL iom_rstput( 0, 0, inum0, 'gphit', gphit, ktype = jp_r 4) ! ! longitude79 CALL iom_rstput( 0, 0, inum0, 'gphiu', gphiu, ktype = jp_r 4)80 CALL iom_rstput( 0, 0, inum0, 'gphiv', gphiv, ktype = jp_r 4)81 CALL iom_rstput( 0, 0, inum0, 'gphif', gphif, ktype = jp_r 4)73 CALL iom_rstput( 0, 0, inum0, 'glamt', glamt, ktype = jp_r8 ) ! ! latitude 74 CALL iom_rstput( 0, 0, inum0, 'glamu', glamu, ktype = jp_r8 ) 75 CALL iom_rstput( 0, 0, inum0, 'glamv', glamv, ktype = jp_r8 ) 76 CALL iom_rstput( 0, 0, inum0, 'glamf', glamf, ktype = jp_r8 ) 77 78 CALL iom_rstput( 0, 0, inum0, 'gphit', gphit, ktype = jp_r8 ) ! ! longitude 79 CALL iom_rstput( 0, 0, inum0, 'gphiu', gphiu, ktype = jp_r8 ) 80 CALL iom_rstput( 0, 0, inum0, 'gphiv', gphiv, ktype = jp_r8 ) 81 CALL iom_rstput( 0, 0, inum0, 'gphif', gphif, ktype = jp_r8 ) 82 82 83 83 CALL iom_rstput( 0, 0, inum0, 'e1t', e1t, ktype = jp_r8 ) ! ! e1 scale factors … … 229 229 230 230 ! ! horizontal mesh (inum3) 231 CALL iom_rstput( 0, 0, inum3, 'glamt', glamt, ktype = jp_r 4) ! ! latitude232 CALL iom_rstput( 0, 0, inum3, 'glamu', glamu, ktype = jp_r 4)233 CALL iom_rstput( 0, 0, inum3, 'glamv', glamv, ktype = jp_r 4)234 CALL iom_rstput( 0, 0, inum3, 'glamf', glamf, ktype = jp_r 4)235 236 CALL iom_rstput( 0, 0, inum3, 'gphit', gphit, ktype = jp_r 4) ! ! longitude237 CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu, ktype = jp_r 4)238 CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv, ktype = jp_r 4)239 CALL iom_rstput( 0, 0, inum3, 'gphif', gphif, ktype = jp_r 4)231 CALL iom_rstput( 0, 0, inum3, 'glamt', glamt, ktype = jp_r8 ) ! ! latitude 232 CALL iom_rstput( 0, 0, inum3, 'glamu', glamu, ktype = jp_r8 ) 233 CALL iom_rstput( 0, 0, inum3, 'glamv', glamv, ktype = jp_r8 ) 234 CALL iom_rstput( 0, 0, inum3, 'glamf', glamf, ktype = jp_r8 ) 235 236 CALL iom_rstput( 0, 0, inum3, 'gphit', gphit, ktype = jp_r8 ) ! ! longitude 237 CALL iom_rstput( 0, 0, inum3, 'gphiu', gphiu, ktype = jp_r8 ) 238 CALL iom_rstput( 0, 0, inum3, 'gphiv', gphiv, ktype = jp_r8 ) 239 CALL iom_rstput( 0, 0, inum3, 'gphif', gphif, ktype = jp_r8 ) 240 240 241 241 CALL iom_rstput( 0, 0, inum3, 'e1t', e1t, ktype = jp_r8 ) ! ! e1 scale factors … … 253 253 ! note that mbkt is set to 1 over land ==> use surface tmask 254 254 zprt(:,:) = ssmask(:,:) * REAL( mbkt(:,:) , wp ) 255 CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i 2) ! ! nb of ocean T-points255 CALL iom_rstput( 0, 0, inum4, 'mbathy', zprt, ktype = jp_i4 ) ! ! nb of ocean T-points 256 256 zprt(:,:) = ssmask(:,:) * REAL( mikt(:,:) , wp ) 257 CALL iom_rstput( 0, 0, inum4, 'misf', zprt, ktype = jp_i 2) ! ! nb of ocean T-points257 CALL iom_rstput( 0, 0, inum4, 'misf', zprt, ktype = jp_i4 ) ! ! nb of ocean T-points 258 258 zprt(:,:) = ssmask(:,:) * REAL( risfdep(:,:) , wp ) 259 CALL iom_rstput( 0, 0, inum4, 'isfdraft', zprt, ktype = jp_r 4) ! ! nb of ocean T-points259 CALL iom_rstput( 0, 0, inum4, 'isfdraft', zprt, ktype = jp_r8 ) ! ! nb of ocean T-points 260 260 261 261 IF( ln_sco ) THEN ! s-coordinate … … 279 279 CALL iom_rstput( 0, 0, inum4, 'gdept_1d' , gdept_1d ) ! ! stretched system 280 280 CALL iom_rstput( 0, 0, inum4, 'gdepw_1d' , gdepw_1d ) 281 CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r 4 )282 CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r 4)281 CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r8 ) 282 CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r8 ) 283 283 ENDIF 284 284 … … 302 302 ! 303 303 IF( nmsh <= 3 ) THEN ! ! 3D depth 304 CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r 4 )304 CALL iom_rstput( 0, 0, inum4, 'gdept_0', gdept_0, ktype = jp_r8 ) 305 305 DO jk = 1,jpk 306 306 DO jj = 1, jpjm1 … … 308 308 zdepu(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji+1,jj ,jk) ) 309 309 zdepv(ji,jj,jk) = MIN( gdept_0(ji,jj,jk) , gdept_0(ji ,jj+1,jk) ) 310 END DO 311 END DO 310 END DO 311 END DO 312 312 END DO 313 CALL lbc_lnk( zdepu, 'U', 1. ) ; CALL lbc_lnk( zdepv, 'V', 1. ) 314 CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r 4)315 CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r 4)316 CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r 4)313 CALL lbc_lnk( zdepu, 'U', 1. ) ; CALL lbc_lnk( zdepv, 'V', 1. ) 314 CALL iom_rstput( 0, 0, inum4, 'gdepu', zdepu, ktype = jp_r8 ) 315 CALL iom_rstput( 0, 0, inum4, 'gdepv', zdepv, ktype = jp_r8 ) 316 CALL iom_rstput( 0, 0, inum4, 'gdepw_0', gdepw_0, ktype = jp_r8 ) 317 317 ELSE ! ! 2D bottom depth 318 318 DO jj = 1,jpj … … 322 322 END DO 323 323 END DO 324 CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r 4 )325 CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r 4 )324 CALL iom_rstput( 0, 0, inum4, 'hdept', zprt, ktype = jp_r8 ) 325 CALL iom_rstput( 0, 0, inum4, 'hdepw', zprw, ktype = jp_r8 ) 326 326 ENDIF 327 327 ! -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
r6140 r6579 102 102 IF( ln_tsd_init .OR. ln_tsd_tradmp ) THEN 103 103 ! 104 IF(lwp) WRITE(numout,*) ' sto dentro if ln_tsd_init prima di allocate if necessary riga 104 ' 104 105 ALLOCATE( sf_tsd(jpts), STAT=ierr0 ) 105 106 IF( ierr0 > 0 ) THEN … … 107 108 ENDIF 108 109 ! 110 IF(lwp) WRITE(numout,*) ' riga 110 : allocate jp_temp jpk ' 109 111 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk) , STAT=ierr0 ) 110 112 IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 ) 113 IF(lwp) WRITE(numout,*) ' riga 113 : dopo allocate jp_temp jpk, 2 ' 111 114 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk) , STAT=ierr2 ) 112 115 IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 ) … … 117 120 ! ! fill sf_tsd with sn_tem & sn_sal and control print 118 121 slf_i(jp_tem) = sn_tem ; slf_i(jp_sal) = sn_sal 122 IF(lwp) WRITE(numout,*) ' riga 122 : prima di fld_fill' 119 123 CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd' ) 124 IF(lwp) WRITE(numout,*) ' riga 124 : dopo di fld_fill' 120 125 ! 121 126 ENDIF -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r6140 r6579 14 14 !! 3.3 ! 2010-10 (C. Ethe) merge TRC-TRA 15 15 !! 3.4 ! 2011-04 (G. Madec) Merge of dtatem and dtasal & suppression of tb,tn/sb,sn 16 !! 3.7 ! 2016-04 (S. Flavoni) change configuration's interface: 17 !! read file or CALL usr_def module to compute initial state (example given for GYRE) 16 18 !!---------------------------------------------------------------------- 17 19 18 20 !!---------------------------------------------------------------------- 19 21 !! istate_init : initial state setting 20 !! istate_tem : analytical profile for initial Temperature21 !! istate_sal : analytical profile for initial Salinity22 !! istate_eel : initial state setting of EEL R5 configuration23 !! istate_gyre : initial state setting of GYRE configuration24 22 !! istate_uvg : initial velocity in geostropic balance 25 23 !!---------------------------------------------------------------------- … … 43 41 USE wrk_nemo ! Memory allocation 44 42 USE timing ! Timing 43 USE usrdef ! User defined routine 44 45 45 46 46 47 IMPLICIT NONE … … 48 49 49 50 PUBLIC istate_init ! routine called by step.F90 51 !SF PUBLIC ini_read ! subroutine ini_read 50 52 51 53 !! * Substitutions … … 75 77 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 76 78 79 !SF initialisation of T & S with data file 77 80 CALL dta_tsd_init ! Initialisation of T & S input data 78 81 IF( lk_c1d ) CALL dta_uvd_init ! Initialization of U & V input data … … 91 94 ! ! Start from rest 92 95 ! ! --------------- 93 numror = 0 94 neuler = 0 95 CALL day_init 96 ! 96 numror = 0 ! define numror = 0 -> no restart file to read 97 neuler = 0 ! Set time-step indicator at nit000 (euler forward) 98 CALL day_init ! model calendar (using both namelist and restart infos) 99 ! ! Initialization of ocean to zero 97 100 ! before fields ! now fields 98 101 sshb (:,:) = 0._wp ; sshn (:,:) = 0._wp … … 101 104 hdivn(:,:,:) = 0._wp 102 105 ! 103 IF( cp_cfg == 'eel' ) THEN 104 CALL istate_eel ! EEL configuration : start from pre-defined U,V T-S fields 105 ELSEIF( cp_cfg == 'gyre' ) THEN 106 CALL istate_gyre ! GYRE configuration : start from pre-defined T-S fields 107 ELSE ! Initial T-S, U-V fields read in files 108 IF ( ln_tsd_init ) THEN ! read 3D T and S data at nit000 109 CALL dta_tsd( nit000, tsb ) 110 tsn(:,:,:,:) = tsb(:,:,:,:) 111 ! 112 ELSE ! Initial T-S fields defined analytically 113 CALL istate_t_s 114 ENDIF 115 IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000 116 CALL wrk_alloc( jpi,jpj,jpk,2, zuvd ) 117 CALL dta_uvd( nit000, zuvd ) 118 ub(:,:,:) = zuvd(:,:,:,1) ; un(:,:,:) = ub(:,:,:) 119 vb(:,:,:) = zuvd(:,:,:,2) ; vn(:,:,:) = vb(:,:,:) 120 CALL wrk_dealloc( jpi,jpj,jpk,2, zuvd ) 121 ENDIF 106 IF( ln_tsd_init ) THEN ! read 3D T and S data at nit000 107 CALL dta_tsd( nit000, tsb ) 108 ELSE ! user defined initial T and S 109 CALL usr_def_ini( tsb ) 110 ENDIF 111 tsn(:,:,:,:) = tsb(:,:,:,:) ! set now to before values 112 ! 113 IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000 114 CALL wrk_alloc( jpi,jpj,jpk,2, zuvd ) 115 CALL dta_uvd( nit000, zuvd ) 116 ub(:,:,:) = zuvd(:,:,:,1) ; un(:,:,:) = ub(:,:,:) 117 vb(:,:,:) = zuvd(:,:,:,2) ; vn(:,:,:) = vb(:,:,:) 118 CALL wrk_dealloc( jpi,jpj,jpk,2, zuvd ) 122 119 ENDIF 123 120 ! … … 131 128 !!gm 132 129 ! 133 ENDIF 130 ENDIF 134 131 ! 135 132 ! Initialize "now" and "before" barotropic velocities: … … 162 159 END SUBROUTINE istate_init 163 160 164 165 SUBROUTINE istate_t_s166 !!---------------------------------------------------------------------167 !! *** ROUTINE istate_t_s ***168 !!169 !! ** Purpose : Intialization of the temperature field with an170 !! analytical profile or a file (i.e. in EEL configuration)171 !!172 !! ** Method : - temperature: use Philander analytic profile173 !! - salinity : use to a constant value 35.5174 !!175 !! References : Philander ???176 !!----------------------------------------------------------------------177 INTEGER :: ji, jj, jk178 REAL(wp) :: zsal = 35.50_wp179 !!----------------------------------------------------------------------180 !181 IF(lwp) WRITE(numout,*)182 IF(lwp) WRITE(numout,*) 'istate_t_s : Philander s initial temperature profile'183 IF(lwp) WRITE(numout,*) '~~~~~~~~~~ and constant salinity (',zsal,' psu)'184 !185 DO jk = 1, jpk186 tsn(:,:,jk,jp_tem) = ( ( ( 7.5 - 0. * ABS( gphit(:,:) )/30. ) * ( 1.-TANH((gdept_n(:,:,jk)-80.)/30.) ) &187 & + 10. * ( 5000. - gdept_n(:,:,jk) ) /5000.) ) * tmask(:,:,jk)188 tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)189 END DO190 tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)191 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)192 !193 END SUBROUTINE istate_t_s194 195 196 SUBROUTINE istate_eel197 !!----------------------------------------------------------------------198 !! *** ROUTINE istate_eel ***199 !!200 !! ** Purpose : Initialization of the dynamics and tracers for EEL R5201 !! configuration (channel with or without a topographic bump)202 !!203 !! ** Method : - set temprature field204 !! - set salinity field205 !! - set velocity field including horizontal divergence206 !! and relative vorticity fields207 !!----------------------------------------------------------------------208 USE divhor ! hor. divergence (div_hor routine)209 USE iom210 !211 INTEGER :: inum ! temporary logical unit212 INTEGER :: ji, jj, jk ! dummy loop indices213 INTEGER :: ijloc214 REAL(wp) :: zh1, zh2, zslope, zcst, zfcor ! temporary scalars215 REAL(wp) :: zt1 = 15._wp ! surface temperature value (EEL R5)216 REAL(wp) :: zt2 = 5._wp ! bottom temperature value (EEL R5)217 REAL(wp) :: zsal = 35.0_wp ! constant salinity (EEL R2, R5 and R6)218 REAL(wp) :: zueel = 0.1_wp ! constant uniform zonal velocity (EEL R5)219 REAL(wp), DIMENSION(jpiglo,jpjglo) :: zssh ! initial ssh over the global domain220 !!----------------------------------------------------------------------221 !222 SELECT CASE ( jp_cfg )223 ! ! ====================224 CASE ( 5 ) ! EEL R5 configuration225 ! ! ====================226 !227 ! set temperature field with a linear profile228 ! -------------------------------------------229 IF(lwp) WRITE(numout,*)230 IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: linear temperature profile'231 IF(lwp) WRITE(numout,*) '~~~~~~~~~~'232 !233 zh1 = gdept_1d( 1 )234 zh2 = gdept_1d(jpkm1)235 !236 zslope = ( zt1 - zt2 ) / ( zh1 - zh2 )237 zcst = ( zt1 * ( zh1 - zh2) - ( zt1 - zt2 ) * zh1 ) / ( zh1 - zh2 )238 !239 DO jk = 1, jpk240 tsn(:,:,jk,jp_tem) = ( zt2 + zt1 * exp( - gdept_n(:,:,jk) / 1000 ) ) * tmask(:,:,jk)241 tsb(:,:,jk,jp_tem) = tsn(:,:,jk,jp_tem)242 END DO243 !244 ! set salinity field to a constant value245 ! --------------------------------------246 IF(lwp) WRITE(numout,*)247 IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal248 IF(lwp) WRITE(numout,*) '~~~~~~~~~~'249 !250 tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)251 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)252 !253 ! set the dynamics: U,V, hdiv (and ssh if necessary)254 ! ----------------255 ! Start EEL5 configuration with barotropic geostrophic velocities256 ! according the sshb and sshn SSH imposed.257 ! we assume a uniform grid (hence the use of e1t(1,1) for delta_y)258 ! we use the Coriolis frequency at mid-channel.259 ub(:,:,:) = zueel * umask(:,:,:)260 un(:,:,:) = ub(:,:,:)261 ijloc = mj0(INT(jpjglo-1)/2)262 zfcor = ff(1,ijloc)263 !264 DO jj = 1, jpjglo265 zssh(:,jj) = - (FLOAT(jj)- FLOAT(jpjglo-1)/2.)*zueel*e1t(1,1)*zfcor/grav266 END DO267 !268 IF(lwp) THEN269 WRITE(numout,*) ' Uniform zonal velocity for EEL R5:',zueel270 WRITE(numout,*) ' Geostrophic SSH profile as a function of y:'271 WRITE(numout,'(12(1x,f6.2))') zssh(1,:)272 ENDIF273 !274 DO jj = 1, nlcj275 DO ji = 1, nlci276 sshb(ji,jj) = zssh( mig(ji) , mjg(jj) ) * tmask(ji,jj,1)277 END DO278 END DO279 sshb(nlci+1:jpi, : ) = 0.e0 ! set to zero extra mpp columns280 sshb( : ,nlcj+1:jpj) = 0.e0 ! set to zero extra mpp rows281 !282 sshn(:,:) = sshb(:,:) ! set now ssh to the before value283 !284 IF( nn_rstssh /= 0 ) THEN285 nn_rstssh = 0 ! hand-made initilization of ssh286 CALL ctl_warn( 'istate_eel: force nn_rstssh = 0' )287 ENDIF288 !289 !!gm Check here call to div_hor should not be necessary290 !!gm div_hor call runoffs not sure they are defined at that level291 CALL div_hor( nit000 ) ! horizontal divergence and relative vorticity (curl)292 ! N.B. the vertical velocity will be computed from the horizontal divergence field293 ! in istate by a call to wzv routine294 295 296 ! ! ==========================297 CASE ( 2 , 6 ) ! EEL R2 or R6 configuration298 ! ! ==========================299 !300 ! set temperature field with a NetCDF file301 ! ----------------------------------------302 IF(lwp) WRITE(numout,*)303 IF(lwp) WRITE(numout,*) 'istate_eel : EEL R2 or R6: read initial temperature in a NetCDF file'304 IF(lwp) WRITE(numout,*) '~~~~~~~~~~'305 !306 CALL iom_open ( 'eel.initemp', inum )307 CALL iom_get ( inum, jpdom_data, 'initemp', tsb(:,:,:,jp_tem) ) ! read before temprature (tb)308 CALL iom_close( inum )309 !310 tsn(:,:,:,jp_tem) = tsb(:,:,:,jp_tem) ! set nox temperature to tb311 !312 ! set salinity field to a constant value313 ! --------------------------------------314 IF(lwp) WRITE(numout,*)315 IF(lwp) WRITE(numout,*) 'istate_eel : EEL R5: constant salinity field, S = ', zsal316 IF(lwp) WRITE(numout,*) '~~~~~~~~~~'317 !318 tsn(:,:,:,jp_sal) = zsal * tmask(:,:,:)319 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)320 !321 ! ! ===========================322 CASE DEFAULT ! NONE existing configuration323 ! ! ===========================324 WRITE(ctmp1,*) 'EEL with a ', jp_cfg,' km resolution is not coded'325 CALL ctl_stop( ctmp1 )326 !327 END SELECT328 !329 END SUBROUTINE istate_eel330 331 332 SUBROUTINE istate_gyre333 !!----------------------------------------------------------------------334 !! *** ROUTINE istate_gyre ***335 !!336 !! ** Purpose : Initialization of the dynamics and tracers for GYRE337 !! configuration (double gyre with rotated domain)338 !!339 !! ** Method : - set temprature field340 !! - set salinity field341 !!----------------------------------------------------------------------342 INTEGER :: ji, jj, jk ! dummy loop indices343 INTEGER :: inum ! temporary logical unit344 INTEGER, PARAMETER :: ntsinit = 0 ! (0/1) (analytical/input data files) T&S initialization345 !!----------------------------------------------------------------------346 !347 SELECT CASE ( ntsinit)348 !349 CASE ( 0 ) ! analytical T/S profil deduced from LEVITUS350 IF(lwp) WRITE(numout,*)351 IF(lwp) WRITE(numout,*) 'istate_gyre : initial analytical T and S profil deduced from LEVITUS '352 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'353 !354 DO jk = 1, jpk355 DO jj = 1, jpj356 DO ji = 1, jpi357 tsn(ji,jj,jk,jp_tem) = ( 16. - 12. * TANH( (gdept_n(ji,jj,jk) - 400) / 700 ) ) &358 & * (-TANH( (500-gdept_n(ji,jj,jk)) / 150 ) + 1) / 2 &359 & + ( 15. * ( 1. - TANH( (gdept_n(ji,jj,jk)-50.) / 1500.) ) &360 & - 1.4 * TANH((gdept_n(ji,jj,jk)-100.) / 100.) &361 & + 7. * (1500. - gdept_n(ji,jj,jk)) / 1500. ) &362 & * (-TANH( (gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2363 tsn(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)364 tsb(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem)365 366 tsn(ji,jj,jk,jp_sal) = ( 36.25 - 1.13 * TANH( (gdept_n(ji,jj,jk) - 305) / 460 ) ) &367 & * (-TANH((500 - gdept_n(ji,jj,jk)) / 150) + 1) / 2 &368 & + ( 35.55 + 1.25 * (5000. - gdept_n(ji,jj,jk)) / 5000. &369 & - 1.62 * TANH( (gdept_n(ji,jj,jk) - 60. ) / 650. ) &370 & + 0.2 * TANH( (gdept_n(ji,jj,jk) - 35. ) / 100. ) &371 & + 0.2 * TANH( (gdept_n(ji,jj,jk) - 1000.) / 5000.) ) &372 & * (-TANH((gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2373 tsn(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)374 tsb(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal)375 END DO376 END DO377 END DO378 !379 CASE ( 1 ) ! T/S data fields read in dta_tem.nc/data_sal.nc files380 IF(lwp) WRITE(numout,*)381 IF(lwp) WRITE(numout,*) 'istate_gyre : initial T and S read from dta_tem.nc/data_sal.nc files'382 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'383 IF(lwp) WRITE(numout,*) ' NetCDF FORMAT'384 385 ! Read temperature field386 ! ----------------------387 CALL iom_open ( 'data_tem', inum )388 CALL iom_get ( inum, jpdom_data, 'votemper', tsn(:,:,:,jp_tem) )389 CALL iom_close( inum )390 391 tsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem) * tmask(:,:,:)392 tsb(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)393 394 ! Read salinity field395 ! -------------------396 CALL iom_open ( 'data_sal', inum )397 CALL iom_get ( inum, jpdom_data, 'vosaline', tsn(:,:,:,jp_sal) )398 CALL iom_close( inum )399 400 tsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) * tmask(:,:,:)401 tsb(:,:,:,jp_sal) = tsn(:,:,:,jp_sal)402 !403 END SELECT404 !405 IF(lwp) THEN406 WRITE(numout,*)407 WRITE(numout,*) ' Initial temperature and salinity profiles:'408 WRITE(numout, "(9x,' level gdept_1d temperature salinity ')" )409 WRITE(numout, "(10x, i4, 3f10.2)" ) ( jk, gdept_1d(jk), tsn(2,2,jk,jp_tem), tsn(2,2,jk,jp_sal), jk = 1, jpk )410 ENDIF411 !412 END SUBROUTINE istate_gyre413 414 415 SUBROUTINE istate_uvg416 !!----------------------------------------------------------------------417 !! *** ROUTINE istate_uvg ***418 !!419 !! ** Purpose : Compute the geostrophic velocities from (tn,sn) fields420 !!421 !! ** Method : Using the hydrostatic hypothesis the now hydrostatic422 !! pressure is computed by integrating the in-situ density from the423 !! surface to the bottom.424 !! p=integral [ rau*g dz ]425 !!----------------------------------------------------------------------426 USE divhor ! hor. divergence (div_hor routine)427 USE lbclnk ! ocean lateral boundary condition (or mpp link)428 !429 INTEGER :: ji, jj, jk ! dummy loop indices430 REAL(wp) :: zmsv, zphv, zmsu, zphu, zalfg ! temporary scalars431 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprn432 !!----------------------------------------------------------------------433 !434 CALL wrk_alloc( jpi,jpj,jpk, zprn)435 !436 IF(lwp) WRITE(numout,*)437 IF(lwp) WRITE(numout,*) 'istate_uvg : Start from Geostrophy'438 IF(lwp) WRITE(numout,*) '~~~~~~~~~~'439 440 ! Compute the now hydrostatic pressure441 ! ------------------------------------442 443 zalfg = 0.5 * grav * rau0444 445 zprn(:,:,1) = zalfg * e3w_n(:,:,1) * ( 1 + rhd(:,:,1) ) ! Surface value446 447 DO jk = 2, jpkm1 ! Vertical integration from the surface448 zprn(:,:,jk) = zprn(:,:,jk-1) &449 & + zalfg * e3w_n(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )450 END DO451 452 ! Compute geostrophic balance453 ! ---------------------------454 DO jk = 1, jpkm1455 DO jj = 2, jpjm1456 DO ji = fs_2, fs_jpim1 ! vertor opt.457 zmsv = 1. / MAX( umask(ji-1,jj+1,jk) + umask(ji ,jj+1,jk) &458 + umask(ji-1,jj ,jk) + umask(ji ,jj ,jk) , 1. )459 zphv = ( zprn(ji ,jj+1,jk) - zprn(ji-1,jj+1,jk) ) * umask(ji-1,jj+1,jk) / e1u(ji-1,jj+1) &460 + ( zprn(ji+1,jj+1,jk) - zprn(ji ,jj+1,jk) ) * umask(ji ,jj+1,jk) / e1u(ji ,jj+1) &461 + ( zprn(ji ,jj ,jk) - zprn(ji-1,jj ,jk) ) * umask(ji-1,jj ,jk) / e1u(ji-1,jj ) &462 + ( zprn(ji+1,jj ,jk) - zprn(ji ,jj ,jk) ) * umask(ji ,jj ,jk) / e1u(ji ,jj )463 zphv = 1. / rau0 * zphv * zmsv * vmask(ji,jj,jk)464 465 zmsu = 1. / MAX( vmask(ji+1,jj ,jk) + vmask(ji ,jj ,jk) &466 + vmask(ji+1,jj-1,jk) + vmask(ji ,jj-1,jk) , 1. )467 zphu = ( zprn(ji+1,jj+1,jk) - zprn(ji+1,jj ,jk) ) * vmask(ji+1,jj ,jk) / e2v(ji+1,jj ) &468 + ( zprn(ji ,jj+1,jk) - zprn(ji ,jj ,jk) ) * vmask(ji ,jj ,jk) / e2v(ji ,jj ) &469 + ( zprn(ji+1,jj ,jk) - zprn(ji+1,jj-1,jk) ) * vmask(ji+1,jj-1,jk) / e2v(ji+1,jj-1) &470 + ( zprn(ji ,jj ,jk) - zprn(ji ,jj-1,jk) ) * vmask(ji ,jj-1,jk) / e2v(ji ,jj-1)471 zphu = 1. / rau0 * zphu * zmsu * umask(ji,jj,jk)472 473 ! Compute the geostrophic velocities474 un(ji,jj,jk) = -2. * zphu / ( ff(ji,jj) + ff(ji ,jj-1) )475 vn(ji,jj,jk) = 2. * zphv / ( ff(ji,jj) + ff(ji-1,jj ) )476 END DO477 END DO478 END DO479 480 IF(lwp) WRITE(numout,*) ' we force to zero bottom velocity'481 482 ! Susbtract the bottom velocity (level jpk-1 for flat bottom case)483 ! to have a zero bottom velocity484 485 DO jk = 1, jpkm1486 un(:,:,jk) = ( un(:,:,jk) - un(:,:,jpkm1) ) * umask(:,:,jk)487 vn(:,:,jk) = ( vn(:,:,jk) - vn(:,:,jpkm1) ) * vmask(:,:,jk)488 END DO489 490 CALL lbc_lnk( un, 'U', -1. )491 CALL lbc_lnk( vn, 'V', -1. )492 493 ub(:,:,:) = un(:,:,:)494 vb(:,:,:) = vn(:,:,:)495 496 !497 !!gm Check here call to div_hor should not be necessary498 !!gm div_hor call runoffs not sure they are defined at that level499 CALL div_hor( nit000 ) ! now horizontal divergence500 !501 CALL wrk_dealloc( jpi,jpj,jpk, zprn)502 !503 END SUBROUTINE istate_uvg504 505 !!=====================================================================506 161 END MODULE istate -
branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/usrdef.F90
r6484 r6579 2 2 !!============================================================================== 3 3 !! *** MODULE usrdef *** 4 !! User defined module: set analytically the minimum information needed to 5 !! define a configuration (domain, initial state, forcing) 6 !! 7 !! === GYRE configuration given as an example === 8 !! 9 !! This module is given as an example, it can be modified 10 !! following the user's desiderata. 11 !! second order accuracy schemes. 4 !! User defined module: used like example to define init, sbc, bathy, ... 12 5 !!============================================================================== 13 !! History : NEMO ! 2016-03 (S. Flavoni) Original code6 !! History : NEMO ! 2016-03 (S. Flavoni) 14 7 !!---------------------------------------------------------------------- 15 8 16 9 !!---------------------------------------------------------------------- 17 !! usr_def_hgr : compute the horizontal mesh 18 !! usr_def_zgr : compute the vertical mesh 19 !! to be defined: 10 !! usr_def_hgr : initialize the horizontal mesh 20 11 !! usr_def_ini : initial state 21 !! usr_def_sbc : compute the surface bounday conditions 22 !! usr_def_xxx : initialize the xxx 12 !! usr_def_sbc : initialize the surface bounday conditions 23 13 !!---------------------------------------------------------------------- 24 14 USE step_oce ! module used in the ocean time stepping module (step.F90) 25 USE mppini ! shared/distributed memory setting (mpp_init routine) 26 USE phycst ! physical constants 15 USE mppini ! shared/distributed memory setting (mpp_init routine) 16 USE phycst ! physical constants 27 17 IMPLICIT NONE 28 18 PRIVATE … … 30 20 31 21 PUBLIC usr_def_hgr 32 PUBLIC usr_def_ zgr22 PUBLIC usr_def_ini 33 23 34 24 !!---------------------------------------------------------------------- 35 !! NEMO/OPA 3.7 , NEMO Consortium (201 6)36 !! $Id: 25 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 26 !! $Id:$ 37 27 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 38 28 !!---------------------------------------------------------------------- 39 29 CONTAINS 40 30 41 42 SUBROUTINE usr_def_hgr( nbench, jp_cfg, kff, pff, & ! Coriolis parameter (if domain not on the sphere) 31 SUBROUTINE usr_def_hgr( nbench, jp_cfg, kff , pff , & ! Coriolis parameter (if domain not on the sphere) 43 32 & pglamt, pglamu, pglamv , pglamf, & ! gridpoints position (required) 44 33 & pgphit, pgphiu, pgphiv , pgphif, & ! 45 34 & pe1t , pe1u , pe1v , pe1f , & ! scale factors (required) 46 35 & pe2t , pe2u , pe2v , pe2f , & ! 47 & ke1e2u_v ) ! u- & v-surfaces (if gridsize reduction is used in strait(s))36 & ke1e2u_v ) ! u- & v-surfaces (if gridsize reduction is used in strait(s)) 48 37 !!------------------------------------------------------------------------------------------------- 49 38 !! *** ROUTINE usr_def_hgr *** … … 76 65 !!------------------------------------------------------------------------------- 77 66 INTEGER :: ji, jj ! dummy loop indices 78 REAL(wp) :: zlam1, zlam0, zcos_alpha, zim1 , zjm1 , ze1 , ze1deg, zf0 79 REAL(wp) :: zphi1, zphi0, zsin_alpha, zim05, zjm05, zbeta, znorme 67 REAL(wp) :: zlam1, zlam0, zcos_alpha, zim1 , zjm1 , ze1 , ze1deg, zf0 ! local scalars 68 REAL(wp) :: zphi1, zphi0, zsin_alpha, zim05, zjm05, zbeta, znorme ! local scalars 80 69 !!------------------------------------------------------------------------------- 81 70 ! … … 91 80 zphi1 = 29._wp 92 81 ! resolution in meters 93 ze1 = 106000. / REAL( jp_cfg , wp ) 82 ze1 = 106000. / REAL( jp_cfg , wp ) 94 83 ! benchmark: forced the resolution to be about 100 km 95 IF( nbench /= 0 ) ze1 = 106000._wp96 zsin_alpha = - SQRT( 2._wp ) * 0.5_wp97 zcos_alpha = SQRT( 2._wp ) * 0.5_wp98 ze1deg = ze1 / (ra * rad)99 IF( nbench /= 0 ) ze1deg = ze1deg / REAL( jp_cfg , wp ) ! benchmark: keep the lat/+lon100 ! ! at the right jp_cfg resolution101 zlam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp )102 zphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp )103 !104 IF( nprint==1 .AND. lwp ) THEN105 WRITE(numout,*) ' ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha106 WRITE(numout,*) ' ze1deg', ze1deg, 'zlam0', zlam0, 'zphi0', zphi0107 ENDIF108 !109 DO jj = 1, jpj110 DO ji = 1, jpi111 zim1 = REAL( ji + nimpp - 1 ) - 1. ; zim05 = REAL( ji + nimpp - 1 ) - 1.5112 zjm1 = REAL( jj + njmpp - 1 ) - 1. ; zjm05 = REAL( jj + njmpp - 1 ) - 1.5113 !114 !glamt(i,j) longitude at T-point115 !gphit(i,j) latitude at T-point116 pglamt(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha117 pgphit(ji,jj) = zphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha118 !119 !glamu(i,j) longitude at U-point120 !gphiu(i,j) latitude at U-point121 pglamu(ji,jj) = zlam0 + zim1 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha122 pgphiu(ji,jj) = zphi0 - zim1 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha123 !124 !glamv(i,j) longitude at V-point125 !gphiv(i,j) latitude at V-point126 pglamv(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm1 * ze1deg * zsin_alpha127 pgphiv(ji,jj) = zphi0 - zim05 * ze1deg * zsin_alpha + zjm1 * ze1deg * zcos_alpha128 !129 !glamf(i,j) longitude at F-point130 !gphif(i,j) latitude at F-point131 pglamf(ji,jj) = zlam0 + zim1 * ze1deg * zcos_alpha + zjm1 * ze1deg * zsin_alpha132 pgphif(ji,jj) = zphi0 - zim1 * ze1deg * zsin_alpha + zjm1 * ze1deg * zcos_alpha133 END DO134 END DO135 !136 ! !== Horizontal scale factors ==! (in meters)137 !138 ! ! constant grid spacing139 pe1t(:,:) = ze1 ; pe2t(:,:) = ze1140 pe1u(:,:) = ze1 ; pe2u(:,:) = ze1141 pe1v(:,:) = ze1 ; pe2v(:,:) = ze1142 pe1f(:,:) = ze1 ; pe2f(:,:) = ze1143 ! ! NO reduction of grid size in some straits144 ke1e2u_v = 0 ! ==>> u_ & v_surfaces will be computed in dom_ghr routine145 !146 !147 148 ! !== Coriolis parameter ==!149 kff = 1 ! indicate not to compute ff afterward150 !151 zbeta = 2. * omega * COS( rad * zphi1 ) / ra ! beta at latitude zphi1152 !SF we overwrite zphi0 (south point in latitude) used just above to define pgphif (value of zphi0=15.5190567531966)153 !SF for computation of coriolis we keep the parameter of XXXX154 zphi0 = 15._wp ! latitude of the most southern grid point155 zf0 = 2. * omega * SIN( rad * zphi0 ) ! compute f0 1st point south156 !157 pff(:,:) = ( zf0 + zbeta * ABS( gphif(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south)158 !159 IF(lwp) WRITE(numout,*) ' beta-plane used. beta = ', zbeta, ' 1/(s.m)'160 !161 IF( nn_timing == 1 ) CALL timing_stop('usr_def_hgr')162 84 ! 163 85 END SUBROUTINE usr_def_hgr 164 86 87 165 88 SUBROUTINE usr_def_zgr() ! Coriolis parameter (if domain not on the sphere) 166 89 ! subroutine for vertical grid 167 90 END SUBROUTINE usr_def_zgr 168 91 92 93 SUBROUTINE usr_def_ini( pts ) 94 !!---------------------------------------------------------------------- 95 !! *** ROUTINE usr_def_ini *** 96 !! 97 !! ** Purpose : Initialization of the dynamics and tracers 98 !! Here GYRE configuration example : (double gyre with rotated domain) 99 !! 100 !! ** Method : - set temprature field 101 !! - set salinity field 102 !!---------------------------------------------------------------------- 103 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: ptsd ! T & S data 104 ! 105 INTEGER :: ji, jj, jk ! dummy loop indices 106 !!---------------------------------------------------------------------- 107 ! 108 IF(lwp) WRITE(numout,*) 109 IF(lwp) WRITE(numout,*) 'usr_def_ini : analytical definition of initial state ' 110 IF(lwp) WRITE(numout,*) ' T and S profiles deduced from LEVITUS ' 111 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 112 ! 113 DO jk = 1, jpk 114 DO jj = 1, jpj 115 DO ji = 1, jpi 116 pts(ji,jj,jk,jp_tem) = ( ( 16. - 12. * TANH( (gdept_n(ji,jj,jk) - 400) / 700 ) ) & 117 & * (-TANH( (500-gdept_n(ji,jj,jk)) / 150 ) + 1) / 2 & 118 & + ( 15. * ( 1. - TANH( (gdept_n(ji,jj,jk)-50.) / 1500.) ) & 119 & - 1.4 * TANH((gdept_n(ji,jj,jk)-100.) / 100.) & 120 & + 7. * (1500. - gdept_n(ji,jj,jk)) / 1500.) & 121 & * (-TANH( (gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2 ) * tmask(ji,jj,jk) 122 123 pts(ji,jj,jk,jp_sal) = ( ( 36.25 - 1.13 * TANH( (gdept_n(ji,jj,jk) - 305) / 460 ) ) & 124 & * (-TANH((500 - gdept_n(ji,jj,jk)) / 150) + 1) / 2 & 125 & + ( 35.55 + 1.25 * (5000. - gdept_n(ji,jj,jk)) / 5000. & 126 & - 1.62 * TANH( (gdept_n(ji,jj,jk) - 60. ) / 650. ) & 127 & + 0.2 * TANH( (gdept_n(ji,jj,jk) - 35. ) / 100. ) & 128 & + 0.2 * TANH( (gdept_n(ji,jj,jk) - 1000.) / 5000.) ) & 129 & * (-TANH((gdept_n(ji,jj,jk) - 500) / 150) + 1) / 2 ) * tmask(ji,jj,jk) 130 END DO 131 END DO 132 END DO 133 ! 134 END SUBROUTINE usr_def_ini 135 169 136 !!====================================================================== 170 137 END MODULE usrdef
Note: See TracChangeset
for help on using the changeset viewer.
