Changeset 5338
- Timestamp:
- 2015-06-03T14:36:20+02:00 (9 years ago)
- Location:
- branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM
- Files:
-
- 1 deleted
- 13 edited
- 3 copied
-
CONFIG/SHARED/namelist_ref (modified) (25 diffs)
-
CONFIG/cfg.txt (modified) (1 diff)
-
CONFIG/uspcfg.txt (modified) (1 diff)
-
NEMO/OPA_SRC/DOM/domzgr.F90 (modified) (4 diffs)
-
NEMO/OPA_SRC/DOM/istate.F90 (modified) (1 diff)
-
NEMO/OPA_SRC/DYN/dynadv.F90 (modified) (9 diffs)
-
NEMO/OPA_SRC/DYN/dynkeg.F90 (modified) (9 diffs)
-
NEMO/OPA_SRC/TRA/eosbn2.F90 (modified) (3 diffs)
-
NEMO/OPA_SRC/ZDF/zdfbfr.F90 (modified) (10 diffs)
-
NEMO/OPA_SRC/nemogcm.F90 (modified) (2 diffs)
-
NEMO/OPA_SRC/step.F90 (modified) (5 diffs)
-
NEMO/OPA_SRC/step_oce.F90 (modified) (1 diff)
-
NEMO/OPA_SRC/stopar.F90 (copied) (copied from trunk/NEMOGCM/NEMO/OPA_SRC/stopar.F90)
-
NEMO/OPA_SRC/stopts.F90 (copied) (copied from trunk/NEMOGCM/NEMO/OPA_SRC/stopts.F90)
-
NEMO/OPA_SRC/storng.F90 (copied) (copied from trunk/NEMOGCM/NEMO/OPA_SRC/storng.F90)
-
SETTE/sette.sh (modified) (2 diffs)
-
fcm-make/inc/keys-isomip.cfg (deleted)
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/CONFIG/SHARED/namelist_ref
r5317 r5338 10 10 !! 7 - dynamics (namdyn_adv, namdyn_vor, namdyn_hpg, namdyn_spg, namdyn_ldf) 11 11 !! 8 - Verical physics (namzdf, namzdf_ric, namzdf_tke, namzdf_kpp, namzdf_ddm, namzdf_tmx) 12 !! 9 - diagnostics (namnc4, namtrd, namspr, namflo, namhsb )12 !! 9 - diagnostics (namnc4, namtrd, namspr, namflo, namhsb, namsto) 13 13 !! 10 - miscellaneous (namsol, nammpp, namctl) 14 14 !! 11 - Obs & Assim (namobs, nam_asminc) … … 54 54 !! *** Domain namelists *** 55 55 !!====================================================================== 56 !! namcfg parameters of the configuration 56 !! namcfg parameters of the configuration 57 57 !! namzgr vertical coordinate 58 58 !! namzgr_sco s-coordinate or hybrid z-s-coordinate … … 62 62 ! 63 63 !----------------------------------------------------------------------- 64 &namcfg ! parameters of the configuration 64 &namcfg ! parameters of the configuration 65 65 !----------------------------------------------------------------------- 66 66 cp_cfg = "default" ! name of the configuration … … 76 76 jperio = 0 ! lateral cond. type (between 0 and 6) 77 77 ! = 0 closed ; = 1 cyclic East-West 78 ! = 2 equatorial symmetric ; = 3 North fold T-point pivot 78 ! = 2 equatorial symmetric ; = 3 North fold T-point pivot 79 79 ! = 4 cyclic East-West AND North fold T-point pivot 80 80 ! = 5 North fold F-point pivot 81 81 ! = 6 cyclic East-West AND North fold F-point pivot 82 ln_use_jattr = .false. ! use (T) the file attribute: open_ocean_jstart, if present 82 ln_use_jattr = .false. ! use (T) the file attribute: open_ocean_jstart, if present 83 83 ! in netcdf input files, as the start j-row for reading 84 84 / … … 105 105 !!!!!!! SH94 stretching coefficients (ln_s_sh94 = .true.) 106 106 rn_theta = 6.0 ! surface control parameter (0<=theta<=20) 107 rn_bb = 0.8 ! stretching with SH94 s-sigma 107 rn_bb = 0.8 ! stretching with SH94 s-sigma 108 108 !!!!!!! SF12 stretching coefficient (ln_s_sf12 = .true.) 109 109 rn_alpha = 4.4 ! stretching with SF12 s-sigma … … 114 114 rn_zb_b = -0.2 ! offset for calculating Zb 115 115 !!!!!!!! Other stretching (not SH94 or SF12) [also uses rn_theta above] 116 rn_thetb = 1.0 ! bottom control parameter (0<=thetb<= 1) 116 rn_thetb = 1.0 ! bottom control parameter (0<=thetb<= 1) 117 117 / 118 118 !----------------------------------------------------------------------- … … 122 122 rn_bathy = 0. ! value of the bathymetry. if (=0) bottom flat at jpkm1 123 123 nn_closea = 0 ! remove (=0) or keep (=1) closed seas and lakes (ORCA) 124 nn_msh = 0! create (=1) a mesh file or not (=0)124 nn_msh = 1 ! create (=1) a mesh file or not (=0) 125 125 rn_hmin = -3. ! min depth of the ocean (>0) or min number of ocean level (<0) 126 126 rn_e3zps_min= 20. ! partial step thickness is set larger than the minimum of … … 168 168 nn_baro = 30 ! Number of iterations of barotropic mode 169 169 ! during rn_rdt seconds. Only used if ln_bt_nn_auto=F 170 rn_bt_cmax = 0.8 ! Maximum courant number allowed if ln_bt_nn_auto=T 170 rn_bt_cmax = 0.8 ! Maximum courant number allowed if ln_bt_nn_auto=T 171 171 nn_bt_flt = 1 ! Time filter choice 172 172 ! = 0 None 173 173 ! = 1 Boxcar over nn_baro barotropic steps 174 ! = 2 Boxcar over 2*nn_baro " " 174 ! = 2 Boxcar over 2*nn_baro " " 175 175 / 176 176 !----------------------------------------------------------------------- … … 249 249 ln_rnf = .true. ! runoffs (T => fill namsbc_rnf) 250 250 nn_isf = 0 ! ice shelf melting/freezing (/=0 => fill namsbc_isf) 251 ! 0 =no isf 1 = presence of ISF 252 ! 2 = bg03 parametrisation 3 = rnf file for isf 251 ! 0 =no isf 1 = presence of ISF 252 ! 2 = bg03 parametrisation 3 = rnf file for isf 253 253 ! 4 = ISF fwf specified 254 254 ! option 1 and 4 need ln_isfcav = .true. (domzgr) … … 281 281 &namsbc_flx ! surface boundary condition : flux formulation 282 282 !----------------------------------------------------------------------- 283 ! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! 283 ! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! 284 284 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename ! 285 285 sn_utau = 'utau' , 24 , 'utau' , .false. , .false., 'yearly' , '' , '' , '' … … 324 324 ln_taudif = .false. ! HF tau contribution: use "mean of stress module - module of the mean stress" data 325 325 rn_zqt = 10. ! Air temperature and humidity reference height (m) 326 rn_zu = 10. ! Wind vector reference height (m) 326 rn_zu = 10. ! Wind vector reference height (m) 327 327 rn_pfac = 1. ! multiplicative factor for precipitation (total & snow) 328 328 rn_efac = 1. ! multiplicative factor for evaporation (0. or 1.) 329 rn_vfac = 0. ! multiplicative factor for ocean/ice velocity 329 rn_vfac = 0. ! multiplicative factor for ocean/ice velocity 330 330 ! in the calculation of the wind stress (0.=absolute winds or 1.=relative winds) 331 331 / … … 377 377 ! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! 378 378 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename ! 379 sn_usp = 'sas_grid_U' , 120 , 'vozocrtx' , .true. , .true. , 'yearly' , '' , '' , '' 379 sn_usp = 'sas_grid_U' , 120 , 'vozocrtx' , .true. , .true. , 'yearly' , '' , '' , '' 380 380 sn_vsp = 'sas_grid_V' , 120 , 'vomecrty' , .true. , .true. , 'yearly' , '' , '' , '' 381 381 sn_tem = 'sas_grid_T' , 120 , 'sosstsst' , .true. , .true. , 'yearly' , '' , '' , '' … … 426 426 / 427 427 !----------------------------------------------------------------------- 428 &namsbc_isf ! Top boundary layer (ISF) 428 &namsbc_isf ! Top boundary layer (ISF) 429 429 !----------------------------------------------------------------------- 430 430 ! ! file name ! frequency (hours) ! variable ! time interpol. ! clim ! 'yearly'/ ! weights ! rotation ! … … 503 503 ! Initial mass required for an iceberg of each class 504 504 rn_initial_mass = 8.8e7, 4.1e8, 3.3e9, 1.8e10, 3.8e10, 7.5e10, 1.2e11, 2.2e11, 3.9e11, 7.4e11 505 ! Proportion of calving mass to apportion to each class 505 ! Proportion of calving mass to apportion to each class 506 506 rn_distribution = 0.24, 0.12, 0.15, 0.18, 0.12, 0.07, 0.03, 0.03, 0.03, 0.02 507 507 ! Ratio between effective and real iceberg mass (non-dim) 508 ! i.e. number of icebergs represented at a point 508 ! i.e. number of icebergs represented at a point 509 509 rn_mass_scaling = 2000, 200, 50, 20, 10, 5, 2, 1, 1, 1 510 510 ! thickness of newly calved bergs (m) … … 515 515 rn_bits_erosion_fraction = 0. ! Fraction of erosion melt flux to divert to bergy bits 516 516 rn_sicn_shift = 0. ! Shift of sea-ice concn in erosion flux (0<sicn_shift<1) 517 ln_passive_mode = .false. ! iceberg - ocean decoupling 517 ln_passive_mode = .false. ! iceberg - ocean decoupling 518 518 nn_test_icebergs = 10 ! Create test icebergs of this class (-1 = no) 519 519 ! Put a test iceberg at each gridpoint in box (lon1,lon2,lat1,lat2) 520 520 rn_test_box = 108.0, 116.0, -66.0, -58.0 521 rn_speed_limit = 0. ! CFL speed limit for a berg 521 rn_speed_limit = 0. ! CFL speed limit for a berg 522 522 523 523 ! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! 524 524 ! ! ! (if <0 months) ! name ! (logical) ! (T/F ) ! 'monthly' ! filename ! pairing ! filename ! 525 525 sn_icb = 'calving' , -1 , 'calvingmask', .true. , .true. , 'yearly' , '' , '' , '' 526 527 cn_dir = './' 526 527 cn_dir = './' 528 528 / 529 529 … … 601 601 ! = 2, use tidal harmonic forcing data from files 602 602 ! = 3, use external data AND tidal harmonic forcing 603 cn_dyn3d = 'none' ! 603 cn_dyn3d = 'none' ! 604 604 nn_dyn3d_dta = 0 ! = 0, bdy data are equal to the initial state 605 605 ! = 1, bdy data are read in 'bdydata .nc' files 606 cn_tra = 'none' ! 606 cn_tra = 'none' ! 607 607 nn_tra_dta = 0 ! = 0, bdy data are equal to the initial state 608 608 ! = 1, bdy data are read in 'bdydata .nc' files 609 cn_ice_lim = 'none' ! 609 cn_ice_lim = 'none' ! 610 610 nn_ice_lim_dta = 0 ! = 0, bdy data are equal to the initial state 611 611 ! = 1, bdy data are read in 'bdydata .nc' files … … 616 616 ln_tra_dmp =.false. ! open boudaries conditions for tracers 617 617 ln_dyn3d_dmp =.false. ! open boundary condition for baroclinic velocities 618 rn_time_dmp = 1. ! Damping time scale in days 618 rn_time_dmp = 1. ! Damping time scale in days 619 619 rn_time_dmp_out = 1. ! Outflow damping time scale 620 620 nn_rimwidth = 10 ! width of the relaxation zone … … 669 669 rn_bfri2_max = 1.e-1 ! max. bottom drag coefficient (non linear case and ln_loglayer=T) 670 670 rn_bfeb2 = 2.5e-3 ! bottom turbulent kinetic energy background (m2/s2) 671 rn_bfrz0 = 3.e-3 ! bottom roughness [m] if ln_loglayer=T 671 rn_bfrz0 = 3.e-3 ! bottom roughness [m] if ln_loglayer=T 672 672 ln_bfr2d = .false. ! horizontal variation of the bottom friction coef (read a 2D mask file ) 673 673 rn_bfrien = 50. ! local multiplying factor of bfr (ln_bfr2d=T) … … 715 715 !----------------------------------------------------------------------- 716 716 nn_eos = -1 ! type of equation of state and Brunt-Vaisala frequency 717 ! =-1, TEOS-10 718 ! = 0, EOS-80 717 ! =-1, TEOS-10 718 ! = 0, EOS-80 719 719 ! = 1, S-EOS (simplified eos) 720 720 ln_useCT = .true. ! use of Conservative Temp. ==> surface CT converted in Pot. Temp. in sbcssm … … 807 807 !----------------------------------------------------------------------- 808 808 ln_dynadv_vec = .true. ! vector form (T) or flux form (F) 809 nn_dynkeg = 0 ! scheme for grad(KE): =0 C2 ; =1 Hollingsworth correction 809 810 ln_dynadv_cen2= .false. ! flux form - 2nd order centered scheme 810 811 ln_dynadv_ubs = .false. ! flux form - 3rd order UBS scheme … … 814 815 &nam_vvl ! vertical coordinate options 815 816 !----------------------------------------------------------------------- 816 ln_vvl_zstar = .true. ! zstar vertical coordinate 817 ln_vvl_zstar = .true. ! zstar vertical coordinate 817 818 ln_vvl_ztilde = .false. ! ztilde vertical coordinate: only high frequency variations 818 819 ln_vvl_layer = .false. ! full layer vertical coordinate … … 999 1000 !! namc1d_uvd data: U & V currents ("key_c1d") 1000 1001 !! namc1d_dyndmp U & V newtonian damping ("key_c1d") 1002 !! namsto Stochastic parametrization of EOS 1001 1003 !!====================================================================== 1002 1004 ! … … 1057 1059 ln_dyndmp = .false. ! add a damping term (T) or not (F) 1058 1060 / 1061 !----------------------------------------------------------------------- 1062 &namsto ! Stochastic parametrization of EOS 1063 !----------------------------------------------------------------------- 1064 ln_rststo = .false. ! start from mean parameter (F) or from restart file (T) 1065 ln_rstseed = .true. ! read seed of RNG from restart file 1066 cn_storst_in = "restart_sto" ! suffix of stochastic parameter restart file (input) 1067 cn_storst_out = "restart_sto" ! suffix of stochastic parameter restart file (output) 1068 1069 ln_sto_eos = .false. ! stochastic equation of state 1070 nn_sto_eos = 1 ! number of independent random walks 1071 rn_eos_stdxy = 1.4 ! random walk horz. standard deviation (in grid points) 1072 rn_eos_stdz = 0.7 ! random walk vert. standard deviation (in grid points) 1073 rn_eos_tcor = 1440.0 ! random walk time correlation (in timesteps) 1074 nn_eos_ord = 1 ! order of autoregressive processes 1075 nn_eos_flt = 0 ! passes of Laplacian filter 1076 rn_eos_lim = 2.0 ! limitation factor (default = 3.0) 1077 / 1059 1078 1060 1079 !!====================================================================== … … 1168 1187 ln_sst = .false. ! Logical switch for SST observations 1169 1188 ln_reysst = .false. ! ln_reysst Logical switch for Reynolds observations 1170 ln_ghrsst = .false. ! ln_ghrsst Logical switch for GHRSST observations 1189 ln_ghrsst = .false. ! ln_ghrsst Logical switch for GHRSST observations 1171 1190 1172 1191 ln_sstfb = .false. ! Logical switch for feedback SST data … … 1195 1214 sstfbfiles = 'sst_01.nc' 1196 1215 ! seaicefiles Sea Ice input observation file names 1197 seaicefiles = 'seaice_01.nc' 1216 seaicefiles = 'seaice_01.nc' 1198 1217 ! velavcurfiles Vel. cur. daily av. input file name 1199 1218 ! velhvcurfiles Vel. cur. high freq. input file name -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/CONFIG/cfg.txt
r5313 r5338 8 8 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 9 9 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC 10 ISOMIP OPA_SRC11 10 GYRE OPA_SRC 12 11 ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/CONFIG/uspcfg.txt
r5312 r5338 1 1 ORCA1_CICE # ORCA2_LIM # OPA_SRC TOP_SRC # http://gws-access.ceda.ac.uk/public/nemo/uspconfigs/ORCA1_CICE/v3.6.0/ORCA1_CICE_ctl.txt 2 ISOMIP # GYRE # OPA_SRC # http://gws-access.ceda.ac.uk/public/nemo/uspconfigs/ISOMIP/v3.6.0/ISOMIP_ctl.txt -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5313 r5338 472 472 risfdep(:,:)=0.e0 473 473 misfdep(:,:)=1 474 !475 ! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code476 IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN477 risfdep(:,:)=200.e0478 misfdep(:,:)=1479 ij0 = 1 ; ij1 = 40480 DO jj = mj0(ij0), mj1(ij1)481 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp482 END DO483 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp484 !485 ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN486 !487 risfdep(:,:)=0.e0488 misfdep(:,:)=1489 ij0 = 1 ; ij1 = 40490 DO jj = mj0(ij0), mj1(ij1)491 risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp492 END DO493 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp494 END IF495 474 ! 496 475 DEALLOCATE( idta, zdta ) … … 969 948 !! 970 949 INTEGER :: ji, jj, jk ! dummy loop indices 971 INTEGER :: ik, it 950 INTEGER :: ik, it, ikb, ikt ! temporary integers 972 951 LOGICAL :: ll_print ! Allow control print for debugging 973 952 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points … … 1152 1131 IF ( ln_isfcav ) THEN 1153 1132 ! (ISF) define e3uw (adapted for 2 cells in the water column) 1154 ! Need to test if the modification of only mikt and mbkt levels is enough 1155 DO jk = 2,jpk 1156 DO jj = 1, jpjm1 1157 DO ji = 1, fs_jpim1 ! vector opt. 1158 e3uw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji+1,jj ,jk) ) & 1159 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji+1,jj ,jk-1) ) 1160 e3vw_0(ji,jj,jk) = MIN( gdept_0(ji,jj,jk), gdept_0(ji ,jj+1,jk) ) & 1161 & - MAX( gdept_0(ji,jj,jk-1), gdept_0(ji ,jj+1,jk-1) ) 1162 END DO 1163 END DO 1133 DO jj = 2, jpjm1 1134 DO ji = 2, fs_jpim1 ! vector opt. 1135 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 1136 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 1137 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) & 1138 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) ) 1139 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 1140 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 1141 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) & 1142 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) ) 1143 END DO 1164 1144 END DO 1165 1145 END IF 1166 1146 1167 1147 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 1168 1148 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) … … 1538 1518 1539 1519 ! remove single point "bay" on isf coast line in the ice shelf draft' 1540 DO jk = 1, jpk1520 DO jk = 2, jpk 1541 1521 WHERE (misfdep==0) misfdep=jpk 1542 1522 zmask=0 -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r5313 r5338 110 110 ELSEIF( cp_cfg == 'gyre' ) THEN 111 111 CALL istate_gyre ! GYRE configuration : start from pre-defined T-S fields 112 ELSEIF( cp_cfg == 'isomip' .OR. cp_cfg == 'isomip2') THEN113 IF(lwp) WRITE(numout,*) 'Initialization of T+S for ISOMIP domain'114 tsn(:,:,:,jp_tem)=-1.9*tmask(:,:,:) ! ISOMIP configuration : start from constant T+S fields115 tsn(:,:,:,jp_sal)=34.4*tmask(:,:,:)116 tsb(:,:,:,:)=tsn(:,:,:,:)117 112 ELSE ! Initial T-S, U-V fields read in files 118 113 IF ( ln_tsd_init ) THEN ! read 3D T and S data at nit000 -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
r5313 r5338 5 5 !!============================================================================== 6 6 !! History : 1.0 ! 2006-11 (G. Madec) Original code 7 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 7 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 8 !! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option 8 9 !!---------------------------------------------------------------------- 9 10 … … 17 18 USE dynkeg ! kinetic energy gradient (dyn_keg routine) 18 19 USE dynzad ! vertical advection (dyn_zad routine) 20 ! 19 21 USE in_out_manager ! I/O manager 20 22 USE lib_mpp ! MPP library … … 25 27 26 28 PUBLIC dyn_adv ! routine called by step module 27 PUBLIC dyn_adv_init ! routine called by opa module29 PUBLIC dyn_adv_init ! routine called by opa module 28 30 31 ! !* namdyn_adv namelist * 29 32 LOGICAL, PUBLIC :: ln_dynadv_vec !: vector form flag 33 INTEGER, PUBLIC :: nn_dynkeg !: scheme of kinetic energy gradient: =0 C2 ; =1 Hollingsworth 30 34 LOGICAL, PUBLIC :: ln_dynadv_cen2 !: flux form - 2nd order centered scheme flag 31 35 LOGICAL, PUBLIC :: ln_dynadv_ubs !: flux form - 3rd order UBS scheme flag … … 38 42 # include "vectopt_loop_substitute.h90" 39 43 !!---------------------------------------------------------------------- 40 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)44 !! NEMO/OPA 3.6 , NEMO Consortium (2015) 41 45 !! $Id$ 42 46 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 63 67 SELECT CASE ( nadv ) ! compute advection trend and add it to general trend 64 68 CASE ( 0 ) 65 CALL dyn_keg ( kt ) ! vector form : horizontal gradient of kinetic energy66 CALL dyn_zad ( kt ) ! vector form : vertical advection69 CALL dyn_keg ( kt, nn_dynkeg ) ! vector form : horizontal gradient of kinetic energy 70 CALL dyn_zad ( kt ) ! vector form : vertical advection 67 71 CASE ( 1 ) 68 CALL dyn_keg ( kt ) ! vector form : horizontal gradient of kinetic energy69 CALL dyn_zad_zts ( kt ) ! vector form : vertical advection with sub-timestepping72 CALL dyn_keg ( kt, nn_dynkeg ) ! vector form : horizontal gradient of kinetic energy 73 CALL dyn_zad_zts ( kt ) ! vector form : vertical advection with sub-timestepping 70 74 CASE ( 2 ) 71 CALL dyn_adv_cen2( kt ) ! 2nd order centered scheme75 CALL dyn_adv_cen2( kt ) ! 2nd order centered scheme 72 76 CASE ( 3 ) 73 CALL dyn_adv_ubs ( kt ) ! 3rd order UBS scheme77 CALL dyn_adv_ubs ( kt ) ! 3rd order UBS scheme 74 78 ! 75 CASE (-1 ) ! esopa: test all possibility with control print76 CALL dyn_keg ( kt )79 CASE (-1 ) ! esopa: test all possibility with control print 80 CALL dyn_keg ( kt, nn_dynkeg ) 77 81 CALL dyn_zad ( kt ) 78 82 CALL dyn_adv_cen2( kt ) … … 92 96 !! momentum advection formulation & scheme and set nadv 93 97 !!---------------------------------------------------------------------- 94 INTEGER :: ioptio 95 INTEGER :: ios ! Local integer output status for namelist read 96 !! 97 NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts 98 INTEGER :: ioptio, ios ! Local integer 99 ! 100 NAMELIST/namdyn_adv/ ln_dynadv_vec, nn_dynkeg, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts 98 101 !!---------------------------------------------------------------------- 99 102 ! 100 103 REWIND( numnam_ref ) ! Namelist namdyn_adv in reference namelist : Momentum advection scheme 101 104 READ ( numnam_ref, namdyn_adv, IOSTAT = ios, ERR = 901) … … 112 115 WRITE(numout,*) '~~~~~~~~~~~' 113 116 WRITE(numout,*) ' Namelist namdyn_adv : chose a advection formulation & scheme for momentum' 114 WRITE(numout,*) ' Vector/flux form (T/F) ln_dynadv_vec = ', ln_dynadv_vec 115 WRITE(numout,*) ' 2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2 116 WRITE(numout,*) ' 3rd order UBS advection scheme ln_dynadv_ubs = ', ln_dynadv_ubs 117 WRITE(numout,*) ' Sub timestepping of vertical advection ln_dynzad_zts = ', ln_dynzad_zts 117 WRITE(numout,*) ' Vector/flux form (T/F) ln_dynadv_vec = ', ln_dynadv_vec 118 WRITE(numout,*) ' = 0 standard scheme ; =1 Hollingsworth scheme nn_dynkeg = ', nn_dynkeg 119 WRITE(numout,*) ' 2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2 120 WRITE(numout,*) ' 3rd order UBS advection scheme ln_dynadv_ubs = ', ln_dynadv_ubs 121 WRITE(numout,*) ' Sub timestepping of vertical advection ln_dynzad_zts = ', ln_dynzad_zts 118 122 ENDIF 119 123 … … 126 130 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' ) 127 131 IF( ln_dynzad_zts .AND. .NOT. ln_dynadv_vec ) & 128 129 IF( ln_dynzad_zts .AND. ln_isfcav ) &130 CALL ctl_stop( 'Sub timestepping of vertical advection does not work with ln_isfcav = .TRUE.' )132 CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' ) 133 IF( nn_dynkeg /= nkeg_C2 .AND. nn_dynkeg /= nkeg_HW ) & 134 CALL ctl_stop( 'KEG scheme wrong value of nn_dynkeg' ) 131 135 132 136 ! ! Set nadv … … 139 143 IF(lwp) THEN ! Print the choice 140 144 WRITE(numout,*) 141 IF( nadv == 0 ) WRITE(numout,*) ' vector form : keg + zad + vor is used' 145 IF( nadv == 0 ) WRITE(numout,*) ' vector form : keg + zad + vor is used' 142 146 IF( nadv == 1 ) WRITE(numout,*) ' vector form : keg + zad_zts + vor is used' 147 IF( nadv == 0 .OR. nadv == 1 ) THEN 148 IF( nn_dynkeg == nkeg_C2 ) WRITE(numout,*) 'with Centered standard keg scheme' 149 IF( nn_dynkeg == nkeg_HW ) WRITE(numout,*) 'with Hollingsworth keg scheme' 150 ENDIF 143 151 IF( nadv == 2 ) WRITE(numout,*) ' flux form : 2nd order scheme is used' 144 152 IF( nadv == 3 ) WRITE(numout,*) ' flux form : UBS scheme is used' -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
r5312 r5338 4 4 !! Ocean dynamics: kinetic energy gradient trend 5 5 !!====================================================================== 6 !! History : 1.0 ! 87-09 (P. Andrich, m.-a. Foujols) Original code 7 !! 7.0 ! 97-05 (G. Madec) Split dynber into dynkeg and dynhpg 8 !! 9.0 ! 02-07 (G. Madec) F90: Free form and module 6 !! History : 1.0 ! 1987-09 (P. Andrich, M.-A. Foujols) Original code 7 !! 7.0 ! 1997-05 (G. Madec) Split dynber into dynkeg and dynhpg 8 !! NEMO 1.0 ! 2002-07 (G. Madec) F90: Free form and module 9 !! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option 9 10 !!---------------------------------------------------------------------- 10 11 … … 18 19 ! 19 20 USE in_out_manager ! I/O manager 21 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 20 22 USE lib_mpp ! MPP library 21 23 USE prtctl ! Print control … … 28 30 PUBLIC dyn_keg ! routine called by step module 29 31 32 INTEGER, PARAMETER, PUBLIC :: nkeg_C2 = 0 !: 2nd order centered scheme (standard scheme) 33 INTEGER, PARAMETER, PUBLIC :: nkeg_HW = 1 !: Hollingsworth et al., QJRMS, 1983 34 ! 35 REAL(wp) :: r1_48 = 1._wp / 48._wp !: =1/(4*2*6) 36 30 37 !! * Substitutions 31 38 # include "vectopt_loop_substitute.h90" 32 39 !!---------------------------------------------------------------------- 33 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)40 !! NEMO/OPA 3.6 , NEMO Consortium (2015) 34 41 !! $Id$ 35 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 37 44 CONTAINS 38 45 39 SUBROUTINE dyn_keg( kt )46 SUBROUTINE dyn_keg( kt, kscheme ) 40 47 !!---------------------------------------------------------------------- 41 48 !! *** ROUTINE dyn_keg *** … … 45 52 !! general momentum trend. 46 53 !! 47 !! ** Method : Compute the now horizontal kinetic energy 54 !! ** Method : * kscheme = nkeg_C2 : 2nd order centered scheme that 55 !! conserve kinetic energy. Compute the now horizontal kinetic energy 48 56 !! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ] 57 !! * kscheme = nkeg_HW : Hollingsworth correction following 58 !! Arakawa (2001). The now horizontal kinetic energy is given by: 59 !! zhke = 1/6 [ mi-1( 2 * un^2 + ((un(j+1)+un(j-1))/2)^2 ) 60 !! + mj-1( 2 * vn^2 + ((vn(i+1)+vn(i-1))/2)^2 ) ] 61 !! 49 62 !! Take its horizontal gradient and add it to the general momentum 50 63 !! trend (ua,va). … … 54 67 !! ** Action : - Update the (ua, va) with the hor. ke gradient trend 55 68 !! - send this trends to trd_dyn (l_trddyn=T) for post-processing 69 !! 70 !! ** References : Arakawa, A., International Geophysics 2001. 71 !! Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983. 56 72 !!---------------------------------------------------------------------- 57 INTEGER, INTENT( in ) :: kt ! ocean time-step index 73 INTEGER, INTENT( in ) :: kt ! ocean time-step index 74 INTEGER, INTENT( in ) :: kscheme ! =0/1 type of KEG scheme 58 75 ! 59 76 INTEGER :: ji, jj, jk ! dummy loop indices … … 63 80 !!---------------------------------------------------------------------- 64 81 ! 65 IF( nn_timing == 1 ) CALL timing_start('dyn_keg')82 IF( nn_timing == 1 ) CALL timing_start('dyn_keg') 66 83 ! 67 CALL wrk_alloc( jpi, jpj, jpk,zhke )84 CALL wrk_alloc( jpi,jpj,jpk, zhke ) 68 85 ! 69 86 IF( kt == nit000 ) THEN 70 87 IF(lwp) WRITE(numout,*) 71 IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend '88 IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme 72 89 IF(lwp) WRITE(numout,*) '~~~~~~~' 73 90 ENDIF 74 91 75 92 IF( l_trddyn ) THEN ! Save ua and va trends 76 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )93 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 77 94 ztrdu(:,:,:) = ua(:,:,:) 78 95 ztrdv(:,:,:) = va(:,:,:) 79 96 ENDIF 80 97 81 ! ! =============== 82 DO jk = 1, jpkm1 ! Horizontal slab 83 ! ! =============== 84 DO jj = 2, jpj ! Horizontal kinetic energy at T-point 85 DO ji = fs_2, jpi ! vector opt. 86 zu = 0.25 * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 87 & + un(ji ,jj ,jk) * un(ji ,jj ,jk) ) 88 zv = 0.25 * ( vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 89 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) 90 zhke(ji,jj,jk) = zv + zu 91 !!gm simplier coding ==>> ~ faster 92 ! don't forget to suppress local zu zv scalars 93 ! zhke(ji,jj,jk) = 0.25 * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 94 ! & + un(ji ,jj ,jk) * un(ji ,jj ,jk) & 95 ! & + vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 96 ! & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) 97 !!gm end <<== 98 END DO 99 END DO 100 DO jj = 2, jpjm1 ! add the gradient of kinetic energy to the general momentum trends 98 zhke(:,:,jpk) = 0._wp 99 100 SELECT CASE ( kscheme ) !== Horizontal kinetic energy at T-point ==! 101 ! 102 CASE ( nkeg_C2 ) !-- Standard scheme --! 103 DO jk = 1, jpkm1 104 DO jj = 2, jpj 105 DO ji = fs_2, jpi ! vector opt. 106 zu = un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 107 & + un(ji ,jj ,jk) * un(ji ,jj ,jk) 108 zv = vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 109 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) 110 zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 111 END DO 112 END DO 113 END DO 114 ! 115 CASE ( nkeg_HW ) !-- Hollingsworth scheme --! 116 DO jk = 1, jpkm1 117 DO jj = 2, jpjm1 118 DO ji = fs_2, jpim1 ! vector opt. 119 zu = 8._wp * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 120 & + un(ji ,jj ,jk) * un(ji ,jj ,jk) ) & 121 & + ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) & 122 & + ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) * ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) 123 ! 124 zv = 8._wp * ( vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 125 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) & 126 & + ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) & 127 & + ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) * ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) 128 zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 129 END DO 130 END DO 131 END DO 132 CALL lbc_lnk( zhke, 'T', 1. ) 133 ! 134 END SELECT 135 ! 136 DO jk = 1, jpkm1 !== grad( KE ) added to the general momentum trends ==! 137 DO jj = 2, jpjm1 101 138 DO ji = fs_2, fs_jpim1 ! vector opt. 102 139 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) … … 104 141 END DO 105 142 END DO 106 !!gm idea to be tested ==>> is it faster on scalar computers ? 107 ! DO jj = 2, jpjm1 ! add the gradient of kinetic energy to the general momentum trends 108 ! DO ji = fs_2, fs_jpim1 ! vector opt. 109 ! ua(ji,jj,jk) = ua(ji,jj,jk) - 0.25 * ( + un(ji+1,jj ,jk) * un(ji+1,jj ,jk) & 110 ! & + vn(ji+1,jj-1,jk) * vn(ji+1,jj-1,jk) & 111 ! & + vn(ji+1,jj ,jk) * vn(ji+1,jj ,jk) & 112 ! ! 113 ! & - un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 114 ! & - vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 115 ! & - vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) / e1u(ji,jj) 116 ! ! 117 ! va(ji,jj,jk) = va(ji,jj,jk) - 0.25 * ( un(ji-1,jj+1,jk) * un(ji-1,jj+1,jk) & 118 ! & + un(ji ,jj+1,jk) * un(ji ,jj+1,jk) & 119 ! & + vn(ji ,jj+1,jk) * vn(ji ,jj+1,jk) & 120 ! ! 121 ! & - un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 122 ! & - un(ji ,jj ,jk) * un(ji ,jj ,jk) & 123 ! & - vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) / e2v(ji,jj) 124 ! END DO 125 ! END DO 126 !!gm en idea <<== 127 ! ! =============== 128 END DO ! End of slab 129 ! ! =============== 130 131 IF( l_trddyn ) THEN ! save the Kinetic Energy trends for diagnostic 143 END DO 144 ! 145 IF( l_trddyn ) THEN ! save the Kinetic Energy trends for diagnostic 132 146 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 133 147 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 134 148 CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt ) 135 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )149 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 136 150 ENDIF 137 151 ! … … 139 153 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 140 154 ! 141 CALL wrk_dealloc( jpi, jpj, jpk,zhke )155 CALL wrk_dealloc( jpi,jpj,jpk, zhke ) 142 156 ! 143 IF( nn_timing == 1 ) CALL timing_stop('dyn_keg')157 IF( nn_timing == 1 ) CALL timing_stop('dyn_keg') 144 158 ! 145 159 END SUBROUTINE dyn_keg -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r5313 r5338 47 47 USE lbclnk ! ocean lateral boundary conditions 48 48 USE timing ! Timing 49 USE stopar ! Stochastic T/S fluctuations 50 USE stopts ! Stochastic T/S fluctuations 49 51 50 52 IMPLICIT NONE … … 313 315 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pdep ! depth [m] 314 316 ! 315 INTEGER :: ji, jj, jk ! dummy loop indices 316 REAL(wp) :: zt , zh , zs , ztm ! local scalars 317 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 317 INTEGER :: ji, jj, jk, jsmp ! dummy loop indices 318 INTEGER :: jdof 319 REAL(wp) :: zt , zh , zstemp, zs , ztm ! local scalars 320 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 321 REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign ! local vectors 318 322 !!---------------------------------------------------------------------- 319 323 ! … … 324 328 CASE( -1, 0 ) !== polynomial TEOS-10 / EOS-80 ==! 325 329 ! 326 DO jk = 1, jpkm1 327 DO jj = 1, jpj 328 DO ji = 1, jpi 329 ! 330 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 331 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 332 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 333 ztm = tmask(ji,jj,jk) ! tmask 334 ! 335 zn3 = EOS013*zt & 336 & + EOS103*zs+EOS003 337 ! 338 zn2 = (EOS022*zt & 339 & + EOS112*zs+EOS012)*zt & 340 & + (EOS202*zs+EOS102)*zs+EOS002 341 ! 342 zn1 = (((EOS041*zt & 343 & + EOS131*zs+EOS031)*zt & 344 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 345 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 346 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 347 ! 348 zn0 = (((((EOS060*zt & 349 & + EOS150*zs+EOS050)*zt & 350 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 351 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 352 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 353 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 354 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 355 ! 356 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 357 ! 358 prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface 359 ! 360 prd(ji,jj,jk) = ( zn * r1_rau0 - 1._wp ) * ztm ! density anomaly (masked) 330 ! Stochastic equation of state 331 IF ( ln_sto_eos ) THEN 332 ALLOCATE(zn0_sto(1:2*nn_sto_eos)) 333 ALLOCATE(zn_sto(1:2*nn_sto_eos)) 334 ALLOCATE(zsign(1:2*nn_sto_eos)) 335 DO jsmp = 1, 2*nn_sto_eos, 2 336 zsign(jsmp) = 1._wp 337 zsign(jsmp+1) = -1._wp 338 END DO 339 ! 340 DO jk = 1, jpkm1 341 DO jj = 1, jpj 342 DO ji = 1, jpi 343 ! 344 ! compute density (2*nn_sto_eos) times: 345 ! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts) 346 ! (2) for t-dt, s-ds (with the opposite fluctuation) 347 DO jsmp = 1, nn_sto_eos*2 348 jdof = (jsmp + 1) / 2 349 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 350 zt = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0 ! temperature 351 zstemp = pts (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp) 352 zs = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 ) ! square root salinity 353 ztm = tmask(ji,jj,jk) ! tmask 354 ! 355 zn3 = EOS013*zt & 356 & + EOS103*zs+EOS003 357 ! 358 zn2 = (EOS022*zt & 359 & + EOS112*zs+EOS012)*zt & 360 & + (EOS202*zs+EOS102)*zs+EOS002 361 ! 362 zn1 = (((EOS041*zt & 363 & + EOS131*zs+EOS031)*zt & 364 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 365 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 366 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 367 ! 368 zn0_sto(jsmp) = (((((EOS060*zt & 369 & + EOS150*zs+EOS050)*zt & 370 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 371 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 372 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 373 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 374 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 375 ! 376 zn_sto(jsmp) = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp) 377 END DO 378 ! 379 ! compute stochastic density as the mean of the (2*nn_sto_eos) densities 380 prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp 381 DO jsmp = 1, nn_sto_eos*2 382 prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp) ! potential density referenced at the surface 383 ! 384 prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_rau0 - 1._wp ) ! density anomaly (masked) 385 END DO 386 prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos 387 prd (ji,jj,jk) = 0.5_wp * prd (ji,jj,jk) * ztm / nn_sto_eos 388 END DO 361 389 END DO 362 390 END DO 363 END DO 364 ! 391 DEALLOCATE(zn0_sto,zn_sto,zsign) 392 ! Non-stochastic equation of state 393 ELSE 394 DO jk = 1, jpkm1 395 DO jj = 1, jpj 396 DO ji = 1, jpi 397 ! 398 zh = pdep(ji,jj,jk) * r1_Z0 ! depth 399 zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature 400 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 401 ztm = tmask(ji,jj,jk) ! tmask 402 ! 403 zn3 = EOS013*zt & 404 & + EOS103*zs+EOS003 405 ! 406 zn2 = (EOS022*zt & 407 & + EOS112*zs+EOS012)*zt & 408 & + (EOS202*zs+EOS102)*zs+EOS002 409 ! 410 zn1 = (((EOS041*zt & 411 & + EOS131*zs+EOS031)*zt & 412 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 413 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 414 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 415 ! 416 zn0 = (((((EOS060*zt & 417 & + EOS150*zs+EOS050)*zt & 418 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 419 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 420 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 421 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 422 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 423 ! 424 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 425 ! 426 prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface 427 ! 428 prd(ji,jj,jk) = ( zn * r1_rau0 - 1._wp ) * ztm ! density anomaly (masked) 429 END DO 430 END DO 431 END DO 432 ENDIF 433 365 434 CASE( 1 ) !== simplified EOS ==! 366 435 ! -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
r5313 r5338 171 171 END DO 172 172 END DO 173 CALL lbc_lnk( bfrua, 'U', 1. ) ; CALL lbc_lnk( bfrva, 'V', 1. ) ! Lateral boundary condition 174 173 175 IF ( ln_isfcav ) THEN 174 176 DO jj = 2, jpjm1 175 177 DO ji = 2, jpim1 176 178 ! (ISF) ======================================================================== 177 ikbu = miku(ji,jj) ! ocean bottomlevel at u- and v-points178 ikbv = mikv(ji,jj) ! ( deepest ocean u- and v-points)179 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 180 ikbv = mikv(ji,jj) ! (1st wet ocean u- and v-points) 179 181 ! 180 182 zvu = 0.25 * ( vn(ji,jj ,ikbu) + vn(ji+1,jj ,ikbu) & … … 183 185 & + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv) ) 184 186 ! 185 zecu = SQRT( un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_ bfeb2 )186 zecv = SQRT( vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_ bfeb2 )187 zecu = SQRT( un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_tfeb2 ) 188 zecv = SQRT( vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_tfeb2 ) 187 189 ! 188 190 tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj ) ) * zecu * (1._wp - umask(ji,jj,1)) … … 202 204 END DO 203 205 END DO 204 END IF205 !206 CALL lbc_lnk( bfrua, 'U', 1. ) ; CALL lbc_lnk( bfrva, 'V', 1. ) ! Lateral boundary condition206 CALL lbc_lnk( tfrua, 'U', 1. ) ; CALL lbc_lnk( tfrva, 'V', 1. ) ! Lateral boundary condition 207 END IF 208 ! 207 209 ! 208 210 IF(ln_ctl) CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr - u: ', mask1=umask, & … … 277 279 IF(lwp) WRITE(numout,*) ' coef rn_bfri2 enhancement factor rn_bfrien = ',rn_bfrien 278 280 ENDIF 279 IF(lwp) WRITE(numout,*) ' top friction coef. rn_bfri1 = ', rn_bfri1 280 IF( ln_tfr2d ) THEN 281 IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_tfr2d = ', ln_tfr2d 282 IF(lwp) WRITE(numout,*) ' coef rn_tfri2 enhancement factor rn_tfrien = ',rn_tfrien 283 ENDIF 281 IF ( ln_isfcav ) THEN 282 IF(lwp) WRITE(numout,*) ' top friction coef. rn_bfri1 = ', rn_tfri1 283 IF( ln_tfr2d ) THEN 284 IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_tfr2d = ', ln_tfr2d 285 IF(lwp) WRITE(numout,*) ' coef rn_tfri2 enhancement factor rn_tfrien = ',rn_tfrien 286 ENDIF 287 END IF 284 288 ! 285 289 IF(ln_bfr2d) THEN … … 295 299 bfrua(:,:) = - bfrcoef2d(:,:) 296 300 bfrva(:,:) = - bfrcoef2d(:,:) 301 ! 302 IF ( ln_isfcav ) THEN 303 IF(ln_tfr2d) THEN 304 ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 305 CALL iom_open('tfr_coef.nc',inum) 306 CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array 307 CALL iom_close(inum) 308 tfrcoef2d(:,:) = rn_tfri1 * ( 1 + rn_tfrien * tfrcoef2d(:,:) ) 309 ELSE 310 tfrcoef2d(:,:) = rn_tfri1 ! initialize tfrcoef2d to the namelist variable 311 ENDIF 312 ! 313 tfrua(:,:) = - tfrcoef2d(:,:) 314 tfrva(:,:) = - tfrcoef2d(:,:) 315 END IF 297 316 ! 298 317 CASE( 2 ) … … 311 330 IF(lwp) WRITE(numout,*) ' coef rn_bfri2 enhancement factor rn_bfrien = ',rn_bfrien 312 331 ENDIF 313 IF(lwp) WRITE(numout,*) ' quadratic top friction' 314 IF(lwp) WRITE(numout,*) ' friction coef. rn_bfri2 = ', rn_tfri2 315 IF(lwp) WRITE(numout,*) ' Max. coef. (log case) rn_tfri2_max = ', rn_tfri2_max 316 IF(lwp) WRITE(numout,*) ' background tke rn_tfeb2 = ', rn_tfeb2 317 IF(lwp) WRITE(numout,*) ' log formulation ln_tfr2d = ', ln_loglayer 318 IF(lwp) WRITE(numout,*) ' bottom roughness rn_tfrz0 [m] = ', rn_tfrz0 319 IF( rn_tfrz0<=0.e0 ) THEN 320 WRITE(ctmp1,*) ' bottom roughness must be strictly positive' 321 CALL ctl_stop( ctmp1 ) 322 ENDIF 323 IF( ln_tfr2d ) THEN 324 IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_tfr2d = ', ln_tfr2d 325 IF(lwp) WRITE(numout,*) ' coef rn_tfri2 enhancement factor rn_tfrien = ',rn_tfrien 326 ENDIF 332 IF ( ln_isfcav ) THEN 333 IF(lwp) WRITE(numout,*) ' quadratic top friction' 334 IF(lwp) WRITE(numout,*) ' friction coef. rn_tfri2 = ', rn_tfri2 335 IF(lwp) WRITE(numout,*) ' Max. coef. (log case) rn_tfri2_max = ', rn_tfri2_max 336 IF(lwp) WRITE(numout,*) ' background tke rn_tfeb2 = ', rn_tfeb2 337 IF(lwp) WRITE(numout,*) ' log formulation ln_tfr2d = ', ln_loglayer 338 IF(lwp) WRITE(numout,*) ' top roughness rn_tfrz0 [m] = ', rn_tfrz0 339 IF( rn_tfrz0<=0.e0 ) THEN 340 WRITE(ctmp1,*) ' top roughness must be strictly positive' 341 CALL ctl_stop( ctmp1 ) 342 ENDIF 343 IF( ln_tfr2d ) THEN 344 IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_tfr2d = ', ln_tfr2d 345 IF(lwp) WRITE(numout,*) ' coef rn_tfri2 enhancement factor rn_tfrien = ',rn_tfrien 346 ENDIF 347 END IF 327 348 ! 328 349 IF(ln_bfr2d) THEN … … 336 357 bfrcoef2d(:,:) = rn_bfri2 ! initialize bfrcoef2d to the namelist variable 337 358 ENDIF 359 360 IF ( ln_isfcav ) THEN 361 IF(ln_tfr2d) THEN 362 ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 363 CALL iom_open('tfr_coef.nc',inum) 364 CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array 365 CALL iom_close(inum) 366 ! 367 tfrcoef2d(:,:) = rn_tfri2 * ( 1 + rn_tfrien * tfrcoef2d(:,:) ) 368 ELSE 369 tfrcoef2d(:,:) = rn_tfri2 ! initialize tfrcoef2d to the namelist variable 370 ENDIF 371 END IF 338 372 ! 339 373 IF ( ln_loglayer.AND.(.NOT.lk_vvl) ) THEN ! set "log layer" bottom friction once for all … … 346 380 END DO 347 381 END DO 382 IF ( ln_isfcav ) THEN 383 DO jj = 1, jpj 384 DO ji = 1, jpi 385 ikbt = mikt(ji,jj) 386 ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_tfrz0 ))**2._wp 387 tfrcoef2d(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 388 tfrcoef2d(ji,jj) = MIN(tfrcoef2d(ji,jj), rn_tfri2_max) 389 END DO 390 END DO 391 END IF 348 392 ENDIF 349 393 ! … … 398 442 zminbfr = MIN( zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) ) ) 399 443 zmaxbfr = MAX( zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) ) ) 444 ! (ISF) 445 IF ( ln_isfcav ) THEN 446 ikbu = miku(ji,jj) ! 1st wet ocean level at u- and v-points 447 ikbv = mikv(ji,jj) 448 zfru = 0.5 * fse3u(ji,jj,ikbu) / rdt 449 zfrv = 0.5 * fse3v(ji,jj,ikbv) / rdt 450 IF( ABS( tfrcoef2d(ji,jj) ) > zfru ) THEN 451 IF( ln_ctl ) THEN 452 WRITE(numout,*) 'TFR ', narea, nimpp+ji, njmpp+jj, ikbu 453 WRITE(numout,*) 'TFR ', ABS( tfrcoef2d(ji,jj) ), zfru 454 ENDIF 455 ictu = ictu + 1 456 ENDIF 457 IF( ABS( tfrcoef2d(ji,jj) ) > zfrv ) THEN 458 IF( ln_ctl ) THEN 459 WRITE(numout,*) 'TFR ', narea, nimpp+ji, njmpp+jj, ikbv 460 WRITE(numout,*) 'TFR ', tfrcoef2d(ji,jj), zfrv 461 ENDIF 462 ictv = ictv + 1 463 ENDIF 464 zmintfr = MIN( zmintfr, MIN( zfru, ABS( tfrcoef2d(ji,jj) ) ) ) 465 zmaxtfr = MAX( zmaxtfr, MIN( zfrv, ABS( tfrcoef2d(ji,jj) ) ) ) 466 END IF 467 ! END ISF 400 468 END DO 401 469 END DO … … 405 473 CALL mpp_min( zminbfr ) 406 474 CALL mpp_max( zmaxbfr ) 475 IF ( ln_isfcav) CALL mpp_min( zmintfr ) 476 IF ( ln_isfcav) CALL mpp_max( zmaxtfr ) 407 477 ENDIF 408 478 IF( .NOT.ln_bfrimp) THEN 409 479 IF( lwp .AND. ictu + ictv > 0 ) THEN 410 WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points '411 WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points '480 WRITE(numout,*) ' Bottom/Top friction stability check failed at ', ictu, ' U-points ' 481 WRITE(numout,*) ' Bottom/Top friction stability check failed at ', ictv, ' V-points ' 412 482 WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr 413 WRITE(numout,*) ' Bottomfriction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr414 WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary'483 IF ( ln_isfcav ) WRITE(numout,*) ' Top friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr 484 WRITE(numout,*) ' Bottom/Top friction coefficient will be reduced where necessary' 415 485 ENDIF 416 486 ENDIF -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r5313 r5338 82 82 USE crsini ! initialise grid coarsening utility 83 83 USE lbcnfd, ONLY: isendto, nsndto, nfsloop, nfeloop ! Setup of north fold exchanges 84 USE stopar 85 USE stopts 84 86 85 87 IMPLICIT NONE … … 432 434 IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 ) CALL cla_init ! Cross Land Advection 433 435 CALL icb_init( rdt, nit000) ! initialise icebergs instance 436 CALL sto_par_init ! Stochastic parametrization 437 IF( ln_sto_eos ) CALL sto_pts_init ! RRandom T/S fluctuations 434 438 435 439 #if defined key_top -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/step.F90
r5313 r5338 106 106 107 107 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 108 ! Update stochastic parameters and random T/S fluctuations 109 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 110 CALL sto_par( kstp ) ! Stochastic parameters 111 112 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 108 113 ! Ocean physics update (ua, va, tsa used as workspace) 109 114 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 145 150 ! 146 151 IF( lk_ldfslp ) THEN ! slope of lateral mixing 152 IF(ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 147 153 CALL eos( tsb, rhd, gdept_0(:,:,:) ) ! before in situ density 148 154 IF( ln_zps .AND. .NOT. ln_isfcav) & … … 180 186 ! Note that the computation of vertical velocity above, hence "after" sea level 181 187 ! is necessary to compute momentum advection for the rhs of barotropic loop: 188 IF(ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 182 189 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 183 190 IF( ln_zps .AND. .NOT. ln_isfcav) & … … 261 268 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 262 269 CALL tra_nxt( kstp ) ! tracer fields at next time step 270 IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 263 271 CALL eos ( tsa, rhd, rhop, fsdept_n(:,:,:) ) ! Time-filtered in situ density for hpg computation 264 272 IF( ln_zps .AND. .NOT. ln_isfcav) & … … 271 279 ELSE ! centered hpg (eos then time stepping) 272 280 IF ( .NOT. lk_dynspg_ts ) THEN ! eos already called in time-split case 281 IF( ln_sto_eos ) CALL sto_pts( tsn ) ! Random T/S fluctuations 273 282 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 274 283 IF( ln_zps .AND. .NOT. ln_isfcav) & -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
r5312 r5338 53 53 54 54 USE dynnxt ! time-stepping (dyn_nxt routine) 55 56 USE stopar ! Stochastic parametrization (sto_par routine) 57 USE stopts 55 58 56 59 USE bdy_par ! for lk_bdy -
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/SETTE/sette.sh
r5312 r5338 1000 1000 export TEST_NAME="LONG" 1001 1001 cd ${CONFIG_DIR} 1002 . ./makenemo -m ${CMP_NAM} -n ISOMIP_LONG - rISOMIP -j 8 del_key ${DEL_KEYS}1002 . ./makenemo -m ${CMP_NAM} -n ISOMIP_LONG -u ISOMIP -j 8 del_key ${DEL_KEYS} 1003 1003 cd ${SETTE_DIR} 1004 1004 . ./param.cfg … … 1068 1068 export TEST_NAME="REPRO_1_4" 1069 1069 cd ${CONFIG_DIR} 1070 . ./makenemo -m ${CMP_NAM} -n ISOMIP_4 - rISOMIP -j 8 add_key "key_mpp_rep" del_key ${DEL_KEYS}1070 . ./makenemo -m ${CMP_NAM} -n ISOMIP_4 -u ISOMIP -j 8 add_key "key_mpp_rep" del_key ${DEL_KEYS} 1071 1071 cd ${SETTE_DIR} 1072 1072 . ./param.cfg
Note: See TracChangeset
for help on using the changeset viewer.
