Changeset 9490 for branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC
- Timestamp:
- 2018-04-23T10:44:07+02:00 (6 years ago)
- Location:
- branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 21 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90
r9125 r9490 53 53 WRITE(numout,*) 'dia_tmb_init : Output Top, Middle, Bottom Diagnostics' 54 54 WRITE(numout,*) '~~~~~~~~~~~~' 55 WRITE(numout,*) ' Namelist nam_diatmb : set tmb outputs '56 WRITE(numout,*) ' Switch for TMB diagnostics (T) or not (F) ln_diatmb = ', ln_diatmb55 WRITE(numout,*) ' Namelist nam_diatmb : set tmb outputs ' 56 WRITE(numout,*) ' Switch for TMB diagnostics (T) or not (F) ln_diatmb = ', ln_diatmb 57 57 ENDIF 58 58 ! -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r9436 r9490 148 148 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hu_b , r1_hu_n , r1_hu_a !: inverse of u-depth [1/m] 149 149 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hv_b , r1_hv_n , r1_hv_a !: inverse of v-depth [1/m] 150 151 150 152 151 INTEGER, PUBLIC :: nla10 !: deepest W level Above ~10m (nlb10 - 1) -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r9405 r9490 316 316 IF(lwp) THEN ! control print 317 317 WRITE(numout,*) ' Namelist : namrun --- run parameters' 318 WRITE(numout,*) ' job numbernn_no = ', nn_no318 WRITE(numout,*) ' Assimilation cycle nn_no = ', nn_no 319 319 WRITE(numout,*) ' experiment name for output cn_exp = ', TRIM( cn_exp ) 320 320 WRITE(numout,*) ' file prefix restart input cn_ocerst_in = ', TRIM( cn_ocerst_in ) … … 351 351 ENDIF 352 352 353 no = nn_no ! conversion DOCTOR names into model names (this should disappear soon) 354 cexper = cn_exp 353 cexper = cn_exp ! conversion DOCTOR names into model names (this should disappear soon) 355 354 nrstdt = nn_rstctl 356 355 nit000 = nn_it000 -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
r9168 r9490 9 9 !! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module 10 10 !! 3.3 ! 2010-10 (C. Bricaud, S. Masson) use of fldread 11 !! 3.4 ! 2010-11 (G. Madec, C. Ethe) Merge of dtatem and dtasal + suppression ofCPP keys11 !! 3.4 ! 2010-11 (G. Madec, C. Ethe) Merge of dtatem and dtasal + remove CPP keys 12 12 !!---------------------------------------------------------------------- 13 13 … … 29 29 PUBLIC dta_tsd ! called by istate.F90 and tradmp.90 30 30 31 LOGICAL , PUBLIC :: ln_tsd_init !: T & S data flag 32 LOGICAL , PUBLIC :: ln_tsd_tradmp !: internal damping toward input data flag 31 ! !!* namtsd namelist : Temperature & Salinity Data * 32 LOGICAL , PUBLIC :: ln_tsd_init !: T & S data flag 33 LOGICAL , PUBLIC :: ln_tsd_dmp !: internal damping toward input data flag 33 34 34 35 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tsd ! structure of input SST (file informations, fields read) … … 58 59 TYPE(FLD_N) :: sn_tem, sn_sal 59 60 !! 60 NAMELIST/namtsd/ ln_tsd_init, ln_tsd_ tradmp, cn_dir, sn_tem, sn_sal61 NAMELIST/namtsd/ ln_tsd_init, ln_tsd_dmp, cn_dir, sn_tem, sn_sal 61 62 !!---------------------------------------------------------------------- 62 63 ! … … 72 73 IF(lwm) WRITE ( numond, namtsd ) 73 74 74 IF( PRESENT( ld_tradmp ) ) ln_tsd_ tradmp = .TRUE. ! forces the initialization when tradmp is used75 IF( PRESENT( ld_tradmp ) ) ln_tsd_dmp = .TRUE. ! forces the initialization when tradmp is used 75 76 76 77 IF(lwp) THEN ! control print … … 79 80 WRITE(numout,*) '~~~~~~~~~~~~ ' 80 81 WRITE(numout,*) ' Namelist namtsd' 81 WRITE(numout,*) ' Initialisation of ocean T & S with T &S input data ln_tsd_init 82 WRITE(numout,*) ' damping of ocean T & S toward T &S input data ln_tsd_ tradmp = ', ln_tsd_tradmp82 WRITE(numout,*) ' Initialisation of ocean T & S with T &S input data ln_tsd_init = ', ln_tsd_init 83 WRITE(numout,*) ' damping of ocean T & S toward T &S input data ln_tsd_dmp = ', ln_tsd_dmp 83 84 WRITE(numout,*) 84 IF( .NOT.ln_tsd_init .AND. .NOT.ln_tsd_ tradmp ) THEN85 IF( .NOT.ln_tsd_init .AND. .NOT.ln_tsd_dmp ) THEN 85 86 WRITE(numout,*) 86 WRITE(numout,*) ' T & S data not used'87 WRITE(numout,*) ' ===>> T & S data not used' 87 88 ENDIF 88 89 ENDIF … … 95 96 ! 96 97 ! ! allocate the arrays (if necessary) 97 IF( ln_tsd_init .OR. ln_tsd_tradmp) THEN98 IF( ln_tsd_init .OR. ln_tsd_dmp ) THEN 98 99 ! 99 100 ALLOCATE( sf_tsd(jpts), STAT=ierr0 ) … … 129 130 !! - 'key_orca_lev10' interpolates on 10 times more levels 130 131 !! - s- or mixed z-s coordinate: vertical interpolation on model mesh 131 !! - ln_tsd_ tradmp=F: deallocates the T-S data structure132 !! - ln_tsd_dmp=F: deallocates the T-S data structure 132 133 !! as T-S data are no are used 133 134 !! … … 149 150 ! 150 151 ! !== ORCA_R2 configuration and T & S damping ==! 151 IF( cn_cfg == "orca" .AND. nn_cfg == 2 .AND. ln_tsd_ tradmp ) THEN ! some hand made alterations152 IF( cn_cfg == "orca" .AND. nn_cfg == 2 .AND. ln_tsd_dmp ) THEN ! some hand made alterations 152 153 ! 153 154 ij0 = 101 ; ij1 = 109 ! Reduced T & S in the Alboran Sea … … 238 239 ENDIF 239 240 ! 240 IF( .NOT.ln_tsd_ tradmp ) THEN !== deallocate T & S structure ==!241 IF( .NOT.ln_tsd_dmp ) THEN !== deallocate T & S structure ==! 241 242 ! (data used only for initialisation) 242 243 IF(lwp) WRITE(numout,*) 'dta_tsd: deallocte T & S arrays as they are only use to initialize the run' -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
r9190 r9490 98 98 !!---------------------------------------------------------------------- 99 99 ! 100 IF(lwp) THEN 101 WRITE(numout,*) 102 WRITE(numout,*) 'dyn_adv_init : choice/control of the momentum advection scheme' 103 WRITE(numout,*) '~~~~~~~~~~~~' 104 ENDIF 105 ! 100 106 REWIND( numnam_ref ) ! Namelist namdyn_adv in reference namelist : Momentum advection scheme 101 107 READ ( numnam_ref, namdyn_adv, IOSTAT = ios, ERR = 901) … … 107 113 108 114 IF(lwp) THEN ! Namelist print 109 WRITE(numout,*)110 WRITE(numout,*) 'dyn_adv_init : choice/control of the momentum advection scheme'111 WRITE(numout,*) '~~~~~~~~~~~~'112 115 WRITE(numout,*) ' Namelist namdyn_adv : chose a advection formulation & scheme for momentum' 113 116 WRITE(numout,*) ' linear dynamics : no momentum advection ln_dynadv_NONE = ', ln_dynadv_NONE -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r9168 r9490 53 53 PUBLIC dyn_hpg_init ! routine called by opa module 54 54 55 ! !!* Namelist namdyn_hpg : hydrostatic pressure gradient 56 LOGICAL , PUBLIC :: ln_hpg_zco !: z-coordinate - full steps 57 LOGICAL , PUBLIC :: ln_hpg_zps !: z-coordinate - partial steps (interpolation) 58 LOGICAL , PUBLIC :: ln_hpg_sco !: s-coordinate (standard jacobian formulation) 59 LOGICAL , PUBLIC :: ln_hpg_djc !: s-coordinate (Density Jacobian with Cubic polynomial) 60 LOGICAL , PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) 61 LOGICAL , PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf 62 63 INTEGER , PUBLIC :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 55 ! !!* Namelist namdyn_hpg : hydrostatic pressure gradient 56 LOGICAL, PUBLIC :: ln_hpg_zco !: z-coordinate - full steps 57 LOGICAL, PUBLIC :: ln_hpg_zps !: z-coordinate - partial steps (interpolation) 58 LOGICAL, PUBLIC :: ln_hpg_sco !: s-coordinate (standard jacobian formulation) 59 LOGICAL, PUBLIC :: ln_hpg_djc !: s-coordinate (Density Jacobian with Cubic polynomial) 60 LOGICAL, PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) 61 LOGICAL, PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf 62 63 ! !! Flag to control the type of hydrostatic pressure gradient 64 INTEGER, PARAMETER :: np_ERROR =-10 ! error in specification of lateral diffusion 65 INTEGER, PARAMETER :: np_zco = 0 ! z-coordinate - full steps 66 INTEGER, PARAMETER :: np_zps = 1 ! z-coordinate - partial steps (interpolation) 67 INTEGER, PARAMETER :: np_sco = 2 ! s-coordinate (standard jacobian formulation) 68 INTEGER, PARAMETER :: np_djc = 3 ! s-coordinate (Density Jacobian with Cubic polynomial) 69 INTEGER, PARAMETER :: np_prj = 4 ! s-coordinate (Pressure Jacobian scheme) 70 INTEGER, PARAMETER :: np_isf = 5 ! s-coordinate similar to sco modify for isf 71 ! 72 INTEGER, PUBLIC :: nhpg !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 64 73 65 74 !! * Substitutions … … 95 104 ! 96 105 SELECT CASE ( nhpg ) ! Hydrostatic pressure gradient computation 97 CASE ( 0) ; CALL hpg_zco ( kt ) ! z-coordinate98 CASE ( 1) ; CALL hpg_zps ( kt ) ! z-coordinate plus partial steps (interpolation)99 CASE ( 2) ; CALL hpg_sco ( kt ) ! s-coordinate (standard jacobian formulation)100 CASE ( 3) ; CALL hpg_djc ( kt ) ! s-coordinate (Density Jacobian with Cubic polynomial)101 CASE ( 4) ; CALL hpg_prj ( kt ) ! s-coordinate (Pressure Jacobian scheme)102 CASE ( 5) ; CALL hpg_isf ( kt ) ! s-coordinate similar to sco modify for ice shelf106 CASE ( np_zco ) ; CALL hpg_zco ( kt ) ! z-coordinate 107 CASE ( np_zps ) ; CALL hpg_zps ( kt ) ! z-coordinate plus partial steps (interpolation) 108 CASE ( np_sco ) ; CALL hpg_sco ( kt ) ! s-coordinate (standard jacobian formulation) 109 CASE ( np_djc ) ; CALL hpg_djc ( kt ) ! s-coordinate (Density Jacobian with Cubic polynomial) 110 CASE ( np_prj ) ; CALL hpg_prj ( kt ) ! s-coordinate (Pressure Jacobian scheme) 111 CASE ( np_isf ) ; CALL hpg_isf ( kt ) ! s-coordinate similar to sco modify for ice shelf 103 112 END SELECT 104 113 ! … … 179 188 ENDIF 180 189 ! 181 ! ! Set nhpg from ln_hpg_... flags 182 IF( ln_hpg_zco ) nhpg = 0 183 IF( ln_hpg_zps ) nhpg = 1 184 IF( ln_hpg_sco ) nhpg = 2 185 IF( ln_hpg_djc ) nhpg = 3 186 IF( ln_hpg_prj ) nhpg = 4 187 IF( ln_hpg_isf ) nhpg = 5 188 ! 189 ! ! Consistency check 190 ! ! Set nhpg from ln_hpg_... flags & consistency check 191 nhpg = np_ERROR 190 192 ioptio = 0 191 IF( ln_hpg_zco ) ioptio = ioptio + 1 192 IF( ln_hpg_zps ) ioptio = ioptio + 1 193 IF( ln_hpg_sco ) ioptio = ioptio + 1 194 IF( ln_hpg_djc ) ioptio = ioptio + 1 195 IF( ln_hpg_prj ) ioptio = ioptio + 1 196 IF( ln_hpg_isf ) ioptio = ioptio + 1 193 IF( ln_hpg_zco ) THEN ; nhpg = np_zco ; ioptio = ioptio +1 ; ENDIF 194 IF( ln_hpg_zps ) THEN ; nhpg = np_zps ; ioptio = ioptio +1 ; ENDIF 195 IF( ln_hpg_sco ) THEN ; nhpg = np_sco ; ioptio = ioptio +1 ; ENDIF 196 IF( ln_hpg_djc ) THEN ; nhpg = np_djc ; ioptio = ioptio +1 ; ENDIF 197 IF( ln_hpg_prj ) THEN ; nhpg = np_prj ; ioptio = ioptio +1 ; ENDIF 198 IF( ln_hpg_isf ) THEN ; nhpg = np_isf ; ioptio = ioptio +1 ; ENDIF 199 ! 197 200 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 198 201 ! 202 IF(lwp) THEN 203 WRITE(numout,*) 204 SELECT CASE( nhpg ) 205 CASE( np_zco ) ; WRITE(numout,*) ' ==>>> z-coord. - full steps ' 206 CASE( np_zps ) ; WRITE(numout,*) ' ==>>> z-coord. - partial steps (interpolation)' 207 CASE( np_sco ) ; WRITE(numout,*) ' ==>>> s-coord. (standard jacobian formulation)' 208 CASE( np_djc ) ; WRITE(numout,*) ' ==>>> s-coord. (Density Jacobian: Cubic polynomial)' 209 CASE( np_prj ) ; WRITE(numout,*) ' ==>>> s-coord. (Pressure Jacobian: Cubic polynomial)' 210 CASE( np_isf ) ; WRITE(numout,*) ' ==>>> s-coord. (standard jacobian formulation) for isf' 211 END SELECT 212 WRITE(numout,*) 213 ENDIF 199 214 ! 200 215 IF ( .NOT. ln_isfcav ) THEN !--- no ice shelf load -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
r9190 r9490 17 17 USE phycst ! physical constants 18 18 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 19 USE ldfslp ! lateral diffusion: slopes of mixing orientation20 19 USE dynldf_lap_blp ! lateral mixing (dyn_ldf_lap & dyn_ldf_blp routines) 21 20 USE dynldf_iso ! lateral mixing (dyn_ldf_iso routine ) … … 34 33 PUBLIC dyn_ldf ! called by step module 35 34 PUBLIC dyn_ldf_init ! called by opa module 36 37 ! ! Parameter to control the type of lateral viscous operator38 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 !: error in setting the operator39 INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 !: without operator (i.e. no lateral viscous trend)40 ! !! laplacian ! bilaplacian !41 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 !: iso-level operator42 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 !: iso-neutral or geopotential operator43 44 INTEGER, PUBLIC :: nldf !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals)45 35 46 36 !! * Substitutions … … 72 62 ENDIF 73 63 74 SELECT CASE ( nldf )! compute lateral mixing trend and add it to the general trend64 SELECT CASE ( nldf_dyn ) ! compute lateral mixing trend and add it to the general trend 75 65 ! 76 CASE ( np_lap ) ; CALL dyn_ldf_lap 77 CASE ( np_lap_i ) ; CALL dyn_ldf_iso 78 CASE ( np_blp ) ; CALL dyn_ldf_blp 66 CASE ( np_lap ) ; CALL dyn_ldf_lap( kt, ub, vb, ua, va, 1 ) ! iso-level laplacian 67 CASE ( np_lap_i ) ; CALL dyn_ldf_iso( kt ) ! rotated laplacian 68 CASE ( np_blp ) ; CALL dyn_ldf_blp( kt, ub, vb, ua, va ) ! iso-level bi-laplacian 79 69 ! 80 70 END SELECT … … 101 91 !! ** Purpose : initializations of the horizontal ocean dynamics physics 102 92 !!---------------------------------------------------------------------- 103 INTEGER :: ioptio, ierr ! temporary integers104 !!----------------------------------------------------------------------105 93 ! 106 ! !== Namelist nam_dynldf ==! already read in ldfdyn module 107 ! 108 IF(lwp) THEN !== Namelist print ==! 94 IF(lwp) THEN !== Namelist print ==! 109 95 WRITE(numout,*) 110 96 WRITE(numout,*) 'dyn_ldf_init : Choice of the lateral diffusive operator on dynamics' 111 97 WRITE(numout,*) '~~~~~~~~~~~~' 112 WRITE(numout,*) ' Namelist nam_dynldf : set lateral mixing parameters (type, direction, coefficients)' 113 WRITE(numout,*) ' Type of operator' 114 WRITE(numout,*) ' no explicit diffusion ln_dynldf_NONE = ', ln_dynldf_NONE 115 WRITE(numout,*) ' laplacian operator ln_dynldf_lap = ', ln_dynldf_lap 116 WRITE(numout,*) ' bilaplacian operator ln_dynldf_blp = ', ln_dynldf_blp 117 WRITE(numout,*) ' Direction of action' 118 WRITE(numout,*) ' iso-level ln_dynldf_lev = ', ln_dynldf_lev 119 WRITE(numout,*) ' horizontal (geopotential) ln_dynldf_hor = ', ln_dynldf_hor 120 WRITE(numout,*) ' iso-neutral ln_dynldf_iso = ', ln_dynldf_iso 121 ENDIF 122 ! !== use of lateral operator or not ==! 123 nldf = np_ERROR 124 ioptio = 0 125 IF( ln_dynldf_NONE ) THEN ; nldf = np_no_ldf ; ioptio = ioptio + 1 ; ENDIF 126 IF( ln_dynldf_lap ) THEN ; ioptio = ioptio + 1 ; ENDIF 127 IF( ln_dynldf_blp ) THEN ; ioptio = ioptio + 1 ; ENDIF 128 IF( ioptio /= 1 ) CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' ) 129 ! 130 IF(.NOT.ln_dynldf_NONE ) THEN !== direction ==>> type of operator ==! 131 ioptio = 0 132 IF( ln_dynldf_lev ) ioptio = ioptio + 1 133 IF( ln_dynldf_hor ) ioptio = ioptio + 1 134 IF( ln_dynldf_iso ) ioptio = ioptio + 1 135 IF( ioptio /= 1 ) CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 direction options (level/hor/iso)' ) 98 WRITE(numout,*) ' Namelist namdyn_ldf: already read in ldfdyn module' 99 WRITE(numout,*) ' see ldf_dyn_init report for lateral mixing parameters' 100 WRITE(numout,*) 136 101 ! 137 ! ! Set nldf, the type of lateral diffusion, from ln_dynldf_... logicals 138 ierr = 0 139 IF( ln_dynldf_lap ) THEN ! laplacian operator 140 IF( ln_zco ) THEN ! z-coordinate 141 IF ( ln_dynldf_lev ) nldf = np_lap ! iso-level = horizontal (no rotation) 142 IF ( ln_dynldf_hor ) nldf = np_lap ! iso-level = horizontal (no rotation) 143 IF ( ln_dynldf_iso ) nldf = np_lap_i ! iso-neutral ( rotation) 144 ENDIF 145 IF( ln_zps ) THEN ! z-coordinate with partial step 146 IF ( ln_dynldf_lev ) nldf = np_lap ! iso-level (no rotation) 147 IF ( ln_dynldf_hor ) nldf = np_lap ! iso-level (no rotation) 148 IF ( ln_dynldf_iso ) nldf = np_lap_i ! iso-neutral ( rotation) 149 ENDIF 150 IF( ln_sco ) THEN ! s-coordinate 151 IF ( ln_dynldf_lev ) nldf = np_lap ! iso-level = horizontal (no rotation) 152 IF ( ln_dynldf_hor ) nldf = np_lap_i ! horizontal ( rotation) 153 IF ( ln_dynldf_iso ) nldf = np_lap_i ! iso-neutral ( rotation) 154 ENDIF 155 ENDIF 156 ! 157 IF( ln_dynldf_blp ) THEN ! bilaplacian operator 158 IF( ln_zco ) THEN ! z-coordinate 159 IF( ln_dynldf_lev ) nldf = np_blp ! iso-level = horizontal (no rotation) 160 IF( ln_dynldf_hor ) nldf = np_blp ! iso-level = horizontal (no rotation) 161 IF( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 162 ENDIF 163 IF( ln_zps ) THEN ! z-coordinate with partial step 164 IF( ln_dynldf_lev ) nldf = np_blp ! iso-level (no rotation) 165 IF( ln_dynldf_hor ) nldf = np_blp ! iso-level (no rotation) 166 IF( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 167 ENDIF 168 IF( ln_sco ) THEN ! s-coordinate 169 IF( ln_dynldf_lev ) nldf = np_blp ! iso-level (no rotation) 170 IF( ln_dynldf_hor ) ierr = 2 ! horizontal ( rotation) 171 IF( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 172 ENDIF 173 ENDIF 174 ! 175 IF( ierr == 2 ) CALL ctl_stop( 'rotated bi-laplacian operator does not exist' ) 176 ! 177 IF( nldf == np_lap_i ) l_ldfslp = .TRUE. ! rotation require the computation of the slopes 178 ! 179 ENDIF 180 181 IF(lwp) THEN 182 WRITE(numout,*) 183 IF( nldf == np_no_ldf ) WRITE(numout,*) ' ==>>> NO lateral viscosity' 184 IF( nldf == np_lap ) WRITE(numout,*) ' ==>>> iso-level laplacian operator' 185 IF( nldf == np_lap_i ) WRITE(numout,*) ' ==>>> rotated laplacian operator with iso-level background' 186 IF( nldf == np_blp ) WRITE(numout,*) ' ==>>> iso-level bi-laplacian operator' 102 SELECT CASE( nldf_dyn ) ! print the choice of operator 103 CASE( np_no_ldf ) ; WRITE(numout,*) ' ==>>> NO lateral viscosity' 104 CASE( np_lap ) ; WRITE(numout,*) ' ==>>> iso-level laplacian operator' 105 CASE( np_lap_i ) ; WRITE(numout,*) ' ==>>> rotated laplacian operator with iso-level background' 106 CASE( np_blp ) ; WRITE(numout,*) ' ==>>> iso-level bi-laplacian operator' 107 END SELECT 187 108 ENDIF 188 109 ! -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
r9124 r9490 108 108 REAL(wp) :: zabe1, zmskt, zmkt, zuav, zuwslpi, zuwslpj ! local scalars 109 109 REAL(wp) :: zabe2, zmskf, zmkf, zvav, zvwslpi, zvwslpj ! - - 110 REAL(wp) :: zcof0, zcof1, zcof2, zcof3, zcof4 110 REAL(wp) :: zcof0, zcof1, zcof2, zcof3, zcof4, zaht_0 ! - - 111 111 REAL(wp), DIMENSION(jpi,jpj) :: ziut, zivf, zdku, zdk1u ! 2D workspace 112 112 REAL(wp), DIMENSION(jpi,jpj) :: zjuf, zjvt, zdkv, zdk1v ! - - … … 139 139 ! 140 140 ENDIF 141 141 142 zaht_0 = 0.5_wp * rn_Ud * rn_Ld ! aht_0 from namtra_ldf = zaht_max 143 142 144 ! ! =============== 143 145 DO jk = 1, jpkm1 ! Horizontal slab … … 174 176 & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ) , 1._wp ) 175 177 176 zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )178 zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 177 179 178 180 ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & … … 189 191 & + umask(ji-1,jj,jk+1) + umask(ji,jj,jk ) , 1._wp ) 190 192 191 zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )193 zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 192 194 193 195 ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & … … 206 208 & + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ) , 1._wp ) 207 209 208 zcof2 = - rn_aht_0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )210 zcof2 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 209 211 210 212 zjuf(ji,jj) = ( zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) & … … 227 229 & + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ) , 1._wp ) 228 230 229 zcof1 = - rn_aht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )231 zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 230 232 231 233 zivf(ji,jj) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) & … … 244 246 & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ) , 1._wp ) 245 247 246 zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )248 zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 247 249 248 250 zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & … … 259 261 & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) 260 262 261 zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )263 zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 262 264 263 265 zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & … … 335 337 DO jk = 2, jpkm1 336 338 DO ji = 2, jpim1 337 zcof0 = 0.5_wp * rn_aht_0 * umask(ji,jj,jk)339 zcof0 = 0.5_wp * zaht_0 * umask(ji,jj,jk) 338 340 ! 339 341 zuwslpi = zcof0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) ) … … 353 355 & + zdj1u(ji,jk ) + zdju (ji ,jk ) ) 354 356 ! vertical mixing coefficient (akzu) 355 ! Note: zcof0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0356 akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0357 ! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0 358 akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / zaht_0 357 359 END DO 358 360 END DO … … 361 363 DO jk = 2, jpkm1 362 364 DO ji = 2, jpim1 363 zcof0 = 0.5_wp * rn_aht_0 * vmask(ji,jj,jk)365 zcof0 = 0.5_wp * zaht_0 * vmask(ji,jj,jk) 364 366 ! 365 367 zvwslpi = zcof0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) ) … … 379 381 & + zdjv (ji,jk ) + zdj1v(ji ,jk ) ) 380 382 ! vertical mixing coefficient (akzv) 381 ! Note: zcof0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0382 akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0383 ! Note: zcof0 include zaht_0, so divided by zaht_0 to obtain slp^2 * zaht_0 384 akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / zaht_0 383 385 END DO 384 386 END DO -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
r9250 r9490 19 19 USE zdfdrg ! vertical physics: top/bottom drag coef. 20 20 USE dynadv ,ONLY: ln_dynadv_vec ! dynamics: advection form 21 USE dynldf ,ONLY: nldf, np_lap_i ! dynamics: type of lateral mixing22 21 USE dynldf_iso,ONLY: akzu, akzv ! dynamics: vertical component of rotated lateral mixing 23 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 22 USE ldfdyn ! lateral diffusion: eddy viscosity coef. and type of operator 24 23 USE trd_oce ! trends: ocean variables 25 24 USE trddyn ! trend manager: dynamics … … 156 155 ! !* Matrix construction 157 156 zdt = r2dt * 0.5 158 IF( nldf == np_lap_i ) THEN ! rotated lateral mixing: add its vertical mixing (akzu) 157 SELECT CASE( nldf_dyn ) 158 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 159 159 DO jk = 1, jpkm1 160 160 DO jj = 2, jpjm1 … … 171 171 END DO 172 172 END DO 173 ELSE ! standard case173 CASE DEFAULT ! iso-level lateral mixing 174 174 DO jk = 1, jpkm1 175 175 DO jj = 2, jpjm1 … … 184 184 END DO 185 185 END DO 186 END IF186 END SELECT 187 187 ! 188 188 DO jj = 2, jpjm1 !* Surface boundary conditions … … 274 274 ! !* Matrix construction 275 275 zdt = r2dt * 0.5 276 IF( nldf == np_lap_i ) THEN ! rotated lateral mixing: add its vertical mixing (akzu) 276 SELECT CASE( nldf_dyn ) 277 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 277 278 DO jk = 1, jpkm1 278 279 DO jj = 2, jpjm1 279 280 DO ji = fs_2, fs_jpim1 ! vector opt. 280 281 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 281 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &282 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 282 283 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 283 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &284 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 284 285 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 285 286 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) … … 289 290 END DO 290 291 END DO 291 ELSE ! standard case292 CASE DEFAULT ! iso-level lateral mixing 292 293 DO jk = 1, jpkm1 293 294 DO jj = 2, jpjm1 294 295 DO ji = fs_2, fs_jpim1 ! vector opt. 295 296 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 296 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk )297 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)297 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 298 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 298 299 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 299 300 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) … … 302 303 END DO 303 304 END DO 304 END IF305 END SELECT 305 306 ! 306 307 DO jj = 2, jpjm1 !* Surface boundary conditions -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90
r9367 r9490 28 28 LOGICAL :: ln_rstart !: start from (F) rest or (T) a restart file 29 29 LOGICAL :: ln_rst_list !: output restarts at list of times (T) or by frequency (F) 30 INTEGER :: nn_no !: job number31 30 INTEGER :: nn_rstctl !: control of the time step (0, 1 or 2) 32 31 INTEGER :: nn_rstssh = 0 !: hand made initilization of ssh or not (1/0) … … 46 45 LOGICAL :: ln_xios_read !: use xios to read single file restart 47 46 INTEGER :: nn_wxios !: write resart using xios 0 - no, 1 - single, 2 - multiple file output 47 INTEGER :: nn_no !: Assimilation cycle 48 48 49 49 #if defined key_netcdf4 … … 74 74 75 75 CHARACTER(lc) :: cexper !: experiment name used for output filename 76 INTEGER :: no !: job number77 76 INTEGER :: nrstdt !: control of the time step (0, 1 or 2) 78 77 INTEGER :: nit000 !: index of the first time step … … 122 121 INTEGER :: numstp = -1 !: logical unit for time step 123 122 INTEGER :: numtime = -1 !: logical unit for timing 124 INTEGER :: numout = 6 !: logical unit for output print; Set to stdout to ensure any early125 !output can be collected; do not change123 INTEGER :: numout = 6 !: logical unit for output print; Set to stdout to ensure any 124 ! ! early output can be collected; do not change 126 125 INTEGER :: numnam_ref = -1 !: logical unit for reference namelist 127 126 INTEGER :: numnam_cfg = -1 !: logical unit for configuration specific namelist … … 159 158 160 159 !!---------------------------------------------------------------------- 161 !! NEMO/OPA 3.3 , NEMO Consortium (2010)160 !! NEMO/OPA 4.0 , NEMO Consortium (2018) 162 161 !! $Id$ 163 162 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldfc1d_c2d.F90
r9094 r9490 26 26 PUBLIC ldf_c2d ! called by ldftra and ldfdyn modules 27 27 28 REAL(wp) :: r1_2 = 0.5_wp ! =1/2 29 REAL(wp) :: r1_4 = 0.25_wp ! =1/4 30 REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 28 31 29 !! * Substitutions30 # include "vectopt_loop_substitute.h90"31 32 !!---------------------------------------------------------------------- 32 33 !! NEMO/OPA 3.7 , NEMO Consortium (2015) … … 36 37 CONTAINS 37 38 38 SUBROUTINE ldf_c1d( cd_type, p rat, pahs1, pahs2, pah1, pah2 )39 SUBROUTINE ldf_c1d( cd_type, pahs1, pahs2, pah1, pah2 ) 39 40 !!---------------------------------------------------------------------- 40 41 !! *** ROUTINE ldf_c1d *** … … 43 44 !! 44 45 !! ** Method : 1D eddy diffusivity coefficients F( depth ) 45 !! Reduction by pratfrom surface to bottom46 !! Reduction by zratio from surface to bottom 46 47 !! hyperbolic tangent profile with inflection point 47 48 !! at zh=500m and a width of zw=200m … … 50 51 !! DYN pah1, pah2 defined at T- and F-points 51 52 !!---------------------------------------------------------------------- 52 CHARACTER(len=2) , INTENT(in ) :: cd_type ! DYNamique or TRAcers 53 REAL(wp) , INTENT(in ) :: prat ! ratio surface/deep values [-] 53 CHARACTER(len=3) , INTENT(in ) :: cd_type ! DYNamique or TRAcers 54 54 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pahs1, pahs2 ! surface value of eddy coefficient [m2/s] 55 55 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pah1 , pah2 ! eddy coefficient [m2/s] … … 58 58 REAL(wp) :: zh, zc, zdep1 ! local scalars 59 59 REAL(wp) :: zw , zdep2 ! - - 60 REAL(wp) :: zratio ! - - 60 61 !!---------------------------------------------------------------------- 61 62 IF(lwp) THEN 63 WRITE(numout,*) 64 WRITE(numout,*) ' ldf_c1d : set a given profile to eddy diffusivity/viscosity coefficients' 65 WRITE(numout,*) ' ~~~~~~~' 66 ENDIF 67 62 ! 63 IF(lwp) WRITE(numout,*) 64 IF(lwp) WRITE(numout,*) ' ldf_c1d : set a given profile to eddy mixing coefficients' 65 ! 68 66 ! initialization of the profile 67 zratio = 0.25_wp ! surface/bottom ratio 69 68 zh = 500._wp ! depth of the inflection point [m] 70 69 zw = 1._wp / 200._wp ! width^-1 - - - [1/m] 71 70 ! ! associated coefficient [-] 72 zc = ( 1._wp - prat) / ( 1._wp + TANH( zh * zw) )71 zc = ( 1._wp - zratio ) / ( 1._wp + TANH( zh * zw) ) 73 72 ! 74 73 ! … … 76 75 ! 77 76 CASE( 'DYN' ) ! T- and F-points 78 DO jk = 1, jpk! pah1 at T-point79 pah1(:,:,jk) = pahs1(:,:) * ( prat + zc * ( 1._wp + TANH( - ( gdept_n(:,:,jk) - zh ) * zw) ) ) * tmask(:,:,jk)77 DO jk = jpkm1, 1, -1 ! pah1 at T-point 78 pah1(:,:,jk) = pahs1(:,:) * ( zratio + zc * ( 1._wp + TANH( - ( gdept_0(:,:,jk) - zh ) * zw) ) ) 80 79 END DO 81 DO jk = 1, jpk! pah2 at F-point (zdep2 is an approximation in zps-coord.)80 DO jk = jpkm1, 1, -1 ! pah2 at F-point (zdep2 is an approximation in zps-coord.) 82 81 DO jj = 1, jpjm1 83 DO ji = 1, fs_jpim184 zdep2 = ( gdept_ n(ji,jj+1,jk) + gdept_n(ji+1,jj+1,jk) &85 & + gdept_ n(ji,jj ,jk) + gdept_n(ji+1,jj ,jk) ) * 0.25_wp86 pah2(ji,jj,jk) = pahs2(ji,jj) * ( prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) * fmask(ji,jj,jk)82 DO ji = 1, jpim1 83 zdep2 = ( gdept_0(ji,jj+1,jk) + gdept_0(ji+1,jj+1,jk) & 84 & + gdept_0(ji,jj ,jk) + gdept_0(ji+1,jj ,jk) ) * r1_4 85 pah2(ji,jj,jk) = pahs2(ji,jj) * ( zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) 87 86 END DO 88 87 END DO … … 91 90 ! 92 91 CASE( 'TRA' ) ! U- and V-points (zdep1 & 2 are an approximation in zps-coord.) 93 DO jk = 1, jpk92 DO jk = jpkm1, 1, -1 94 93 DO jj = 1, jpjm1 95 DO ji = 1, fs_jpim196 zdep1 = ( gdept_ n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) * 0.5_wp97 zdep2 = ( gdept_ n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) * 0.5_wp98 pah1(ji,jj,jk) = pahs1(ji,jj) * ( prat + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) ) ) * umask(ji,jj,jk)99 pah2(ji,jj,jk) = pahs2(ji,jj) * ( prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) * vmask(ji,jj,jk)94 DO ji = 1, jpim1 95 zdep1 = ( gdept_0(ji,jj,jk) + gdept_0(ji+1,jj,jk) ) * 0.5_wp 96 zdep2 = ( gdept_0(ji,jj,jk) + gdept_0(ji,jj+1,jk) ) * 0.5_wp 97 pah1(ji,jj,jk) = pahs1(ji,jj) * ( zratio + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) ) ) 98 pah2(ji,jj,jk) = pahs2(ji,jj) * ( zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) 100 99 END DO 101 100 END DO … … 104 103 CALL lbc_lnk_multi( pah1, 'U', 1. , pah2, 'V', 1. ) 105 104 ! 105 CASE DEFAULT ! error 106 CALL ctl_stop( 'ldf_c1d: ', cd_type, ' Unknown, i.e. /= DYN or TRA' ) 106 107 END SELECT 107 108 ! … … 109 110 110 111 111 SUBROUTINE ldf_c2d( cd_type, cd_op, pah0, pah1, pah2 )112 SUBROUTINE ldf_c2d( cd_type, pUfac, knn, pah1, pah2 ) 112 113 !!---------------------------------------------------------------------- 113 114 !! *** ROUTINE ldf_c2d *** … … 124 125 !! DYN pah1, pah2 defined at T- and F-points 125 126 !!---------------------------------------------------------------------- 126 CHARACTER(len=3) 127 CHARACTER(len=3) , INTENT(in ) :: cd_op ! operator:LAPlacian BiLaPlacian128 REAL(wp) , INTENT(in ) :: pah0 ! eddy coefficient [m2/s or m4/s]129 REAL(wp), DIMENSION( jpi,jpj,jpk), INTENT( out) :: pah1, pah2 ! eddy coefficient[m2/s or m4/s]127 CHARACTER(len=3) , INTENT(in ) :: cd_type ! DYNamique or TRAcers 128 REAL(wp) , INTENT(in ) :: pUfac ! =1/2*Uc LAPlacian BiLaPlacian 129 INTEGER , INTENT(in ) :: knn ! characteristic velocity [m/s] 130 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pah1, pah2 ! eddy coefficients [m2/s or m4/s] 130 131 ! 131 132 INTEGER :: ji, jj, jk ! dummy loop indices 132 REAL(wp) :: za00, zd_max, zemax1, zemax2 ! local scalar133 INTEGER :: inn ! local integer 133 134 !!---------------------------------------------------------------------- 134 135 ! 135 zd_max = ra * rad ! = 1 degree at the equator in meters 136 IF(lwp) WRITE(numout,*) 137 IF(lwp) WRITE(numout,*) ' ldf_c2d : aht = Ufac * max(e1,e2) with Ufac = ', pUfac, ' m/s' 136 138 ! 137 IF(lwp) THEN138 WRITE(numout,*)139 WRITE(numout,*) ' ldf_c2d : aht = rn_aht0 * max(e1,e2)/e_equator ( laplacian) '140 WRITE(numout,*) ' ~~~~~~~ or = rn_bht0 * [max(e1,e2)/e_equator]**3 (bilaplacian)'141 WRITE(numout,*)142 ENDIF143 139 ! 144 SELECT CASE( cd_type ) !== surface values ==! ( depending on DYN/TRA)140 SELECT CASE( cd_type ) !== surface values ==! (chosen grid point function of DYN or TRA) 145 141 ! 146 CASE( 'DYN' ) ! T- and F-points 147 IF( cd_op == 'LAP' ) THEN ! laplacian operator 148 IF(lwp) WRITE(numout,*) ' momentum laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)' 149 za00 = pah0 / zd_max 150 DO jj = 1, jpj 151 DO ji = 1, jpi 152 zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1) 153 zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1) 154 pah1(ji,jj,1) = za00 * zemax1 155 pah2(ji,jj,1) = za00 * zemax2 156 END DO 142 CASE( 'DYN' ) ! T- and F-points 143 DO jj = 1, jpj 144 DO ji = 1, jpi 145 pah1(ji,jj,1) = pUfac * MAX( e1t(ji,jj) , e2t(ji,jj) )**knn 146 pah2(ji,jj,1) = pUfac * MAX( e1f(ji,jj) , e2f(ji,jj) )**knn 157 147 END DO 158 ELSEIF( cd_op == 'BLP' ) THEN ! bilaplacian operator 159 IF(lwp) WRITE(numout,*) ' momentum bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3' 160 za00 = pah0 / ( zd_max * zd_max * zd_max ) 161 DO jj = 1, jpj 162 DO ji = 1, jpi 163 zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1) 164 zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1) 165 pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1 166 pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2 167 END DO 148 END DO 149 CASE( 'TRA' ) ! U- and V-points 150 DO jj = 1, jpj 151 DO ji = 1, jpi 152 pah1(ji,jj,1) = pUfac * MAX( e1u(ji,jj), e2u(ji,jj) )**knn 153 pah2(ji,jj,1) = pUfac * MAX( e1v(ji,jj), e2v(ji,jj) )**knn 168 154 END DO 169 ELSE ! NO diffusion/viscosity170 CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' )171 ENDIF172 ! ! deeper values (LAP and BLP cases)173 DO jk = 2, jpk174 pah1(:,:,jk) = pah1(:,:,1) * tmask(:,:,jk)175 pah2(:,:,jk) = pah2(:,:,1) * fmask(:,:,jk)176 155 END DO 177 ! 178 CASE( 'TRA' ) ! U- and V-points (approximation in zps-coord.) 179 IF( cd_op == 'LAP' ) THEN ! laplacian operator 180 IF(lwp) WRITE(numout,*) ' tracer laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)' 181 za00 = pah0 / zd_max 182 DO jj = 1, jpj 183 DO ji = 1, jpi 184 zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1) 185 zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1) 186 pah1(ji,jj,1) = za00 * zemax1 187 pah2(ji,jj,1) = za00 * zemax2 188 END DO 189 END DO 190 ELSEIF( cd_op == 'BLP' ) THEN ! bilaplacian operator (NB: square root of the coeff) 191 IF(lwp) WRITE(numout,*) ' tracer bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3' 192 za00 = pah0 / ( zd_max * zd_max * zd_max ) 193 DO jj = 1, jpj 194 DO ji = 1, jpi 195 zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1) 196 zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1) 197 pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1 198 pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2 199 END DO 200 END DO 201 ELSE ! NO diffusion/viscosity 202 CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' ) 203 ENDIF 204 ! ! deeper values (LAP and BLP cases) 205 DO jk = 2, jpk 206 pah1(:,:,jk) = pah1(:,:,1) * umask(:,:,jk) 207 pah2(:,:,jk) = pah2(:,:,1) * vmask(:,:,jk) 208 END DO 209 ! 156 CASE DEFAULT ! error 157 CALL ctl_stop( 'ldf_c2d: ', cd_type, ' Unknown, i.e. /= DYN or TRA' ) 210 158 END SELECT 159 ! !== deeper values = surface one ==! (except jpk) 160 DO jk = 2, jpkm1 161 pah1(:,:,jk) = pah1(:,:,1) 162 pah2(:,:,jk) = pah2(:,:,1) 163 END DO 211 164 ! 212 165 END SUBROUTINE ldf_c2d -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90
r9169 r9490 17 17 USE dom_oce ! ocean space and time domain 18 18 USE phycst ! physical constants 19 USE ldfslp ! lateral diffusion: slopes of mixing orientation 19 20 USE ldfc1d_c2d ! lateral diffusion: 1D and 2D cases 20 21 ! … … 31 32 PUBLIC ldf_dyn ! called by step.F90 32 33 33 ! 34 ! !!* Namelist namdyn_ldf : lateral mixing on momentum * 34 35 LOGICAL , PUBLIC :: ln_dynldf_NONE !: No operator (i.e. no explicit diffusion) 35 36 LOGICAL , PUBLIC :: ln_dynldf_lap !: laplacian operator … … 37 38 LOGICAL , PUBLIC :: ln_dynldf_lev !: iso-level direction 38 39 LOGICAL , PUBLIC :: ln_dynldf_hor !: horizontal (geopotential) direction 39 LOGICAL , PUBLIC :: ln_dynldf_iso !: iso-neutral direction 40 ! LOGICAL , PUBLIC :: ln_dynldf_iso !: iso-neutral direction (see ldfslp) 40 41 INTEGER , PUBLIC :: nn_ahm_ijk_t !: choice of time & space variations of the lateral eddy viscosity coef. 41 REAL(wp), PUBLIC :: rn_ahm_0 !: lateral laplacian eddy viscosity [m2/s] 42 REAL(wp), PUBLIC :: rn_ahm_b !: lateral laplacian background eddy viscosity [m2/s] 43 REAL(wp), PUBLIC :: rn_bhm_0 !: lateral bilaplacian eddy viscosity [m4/s] 44 !! If nn_ahm_ijk_t = 32 a time and space varying Smagorinsky viscosity 45 !! will be computed. 46 REAL(wp), PUBLIC :: rn_csmc !: Smagorinsky constant of proportionality 47 REAL(wp), PUBLIC :: rn_minfac !: Multiplicative factor of theorectical minimum Smagorinsky viscosity 48 REAL(wp), PUBLIC :: rn_maxfac !: Multiplicative factor of theorectical maximum Smagorinsky viscosity 49 50 LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef. 51 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahmt, ahmf !: eddy diffusivity coef. at U- and V-points [m2/s or m4/s] 42 ! ! time invariant coefficients: aht = 1/2 Ud*Ld (lap case) 43 ! ! bht = 1/12 Ud*Ld^3 (blp case) 44 REAL(wp), PUBLIC :: rn_Uv !: lateral viscous velocity [m/s] 45 REAL(wp), PUBLIC :: rn_Lv !: lateral viscous length [m] 46 ! ! Smagorinsky viscosity (nn_ahm_ijk_t = 32) 47 REAL(wp), PUBLIC :: rn_csmc !: Smagorinsky constant of proportionality 48 REAL(wp), PUBLIC :: rn_minfac !: Multiplicative factor of theorectical minimum Smagorinsky viscosity 49 REAL(wp), PUBLIC :: rn_maxfac !: Multiplicative factor of theorectical maximum Smagorinsky viscosity 50 ! ! iso-neutral laplacian (ln_dynldf_lap=ln_dynldf_iso=T) 51 REAL(wp), PUBLIC :: rn_ahm_b !: lateral laplacian background eddy viscosity [m2/s] 52 53 ! !!* Parameter to control the type of lateral viscous operator 54 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 !: error in setting the operator 55 INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 !: without operator (i.e. no lateral viscous trend) 56 ! !! laplacian ! bilaplacian ! 57 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 !: iso-level operator 58 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 !: iso-neutral or geopotential operator 59 ! 60 INTEGER , PUBLIC :: nldf_dyn !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 61 LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef. 62 63 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahmt, ahmf !: eddy viscosity coef. at T- and F-points [m2/s or m4/s] 53 64 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dtensq !: horizontal tension squared (Smagorinsky only) 54 65 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dshesq !: horizontal shearing strain squared (Smagorinsky only) 55 66 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: esqt, esqf !: Square of the local gridscale (e1e2/(e1+e2))**2 56 67 57 REAL(wp) :: r1_ 12 = 1._wp / 12._wp ! =1/1268 REAL(wp) :: r1_2 = 0.5_wp ! =1/2 58 69 REAL(wp) :: r1_4 = 0.25_wp ! =1/4 59 70 REAL(wp) :: r1_8 = 0.125_wp ! =1/8 71 REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 60 72 REAL(wp) :: r1_288 = 1._wp / 288._wp ! =1/( 12^2 * 2 ) 61 73 … … 92 104 !! or L^4|D|/8 bilaplacian operator ) 93 105 !!---------------------------------------------------------------------- 94 INTEGER :: ji, jj, jk ! dummy loop indices 95 INTEGER :: ierr, inum, ios ! local integer 96 REAL(wp) :: zah0 ! local scalar 97 ! 106 INTEGER :: ji, jj, jk ! dummy loop indices 107 INTEGER :: ioptio, ierr, inum, ios, inn ! local integer 108 REAL(wp) :: zah0, zah_max, zUfac ! local scalar 109 CHARACTER(len=5) :: cl_Units ! units (m2/s or m4/s) 110 !! 98 111 NAMELIST/namdyn_ldf/ ln_dynldf_NONE, ln_dynldf_lap, ln_dynldf_blp, & ! type of operator 99 & ln_dynldf_lev , ln_dynldf_hor, ln_dynldf_iso, & ! acting direction of the operator100 & nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0, & ! lateral eddy coefficient101 & rn_csmc , rn_minfac, rn_maxfac! Smagorinsky settings112 & ln_dynldf_lev , ln_dynldf_hor, ln_dynldf_iso, & ! acting direction of the operator 113 & nn_ahm_ijk_t , rn_Uv , rn_Lv, rn_ahm_b, & ! lateral eddy coefficient 114 & rn_csmc , rn_minfac , rn_maxfac ! Smagorinsky settings 102 115 !!---------------------------------------------------------------------- 103 116 ! … … 129 142 WRITE(numout,*) ' coefficients :' 130 143 WRITE(numout,*) ' type of time-space variation nn_ahm_ijk_t = ', nn_ahm_ijk_t 131 WRITE(numout,*) ' lateral laplacian eddy viscosity rn_ahm_0 = ', rn_ahm_0, ' m2/s' 132 WRITE(numout,*) ' background viscosity (iso case) rn_ahm_b = ', rn_ahm_b, ' m2/s' 133 WRITE(numout,*) ' lateral bilaplacian eddy viscosity rn_bhm_0 = ', rn_bhm_0, ' m4/s' 144 WRITE(numout,*) ' lateral viscous velocity (if cst) rn_Uv = ', rn_Uv, ' m/s' 145 WRITE(numout,*) ' lateral viscous length (if cst) rn_Lv = ', rn_Lv, ' m' 146 WRITE(numout,*) ' background viscosity (iso-lap case) rn_ahm_b = ', rn_ahm_b, ' m2/s' 147 ! 134 148 WRITE(numout,*) ' Smagorinsky settings (nn_ahm_ijk_t = 32) :' 135 149 WRITE(numout,*) ' Smagorinsky coefficient rn_csmc = ', rn_csmc 136 WRITE(numout,*) ' factor multiplier for theorectical lower limit for ' 137 WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_minfac = ', rn_minfac 138 WRITE(numout,*) ' factor multiplier for theorectical lower upper for ' 139 WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_maxfac = ', rn_maxfac 150 WRITE(numout,*) ' factor multiplier for eddy visc.' 151 WRITE(numout,*) ' lower limit (default 1.0) rn_minfac = ', rn_minfac 152 WRITE(numout,*) ' upper limit (default 1.0) rn_maxfac = ', rn_maxfac 140 153 ENDIF 141 154 142 ! ! Parameter control 155 ! 156 ! !== type of lateral operator used ==! (set nldf_dyn) 157 ! !=====================================! 158 ! 159 nldf_dyn = np_ERROR 160 ioptio = 0 161 IF( ln_dynldf_NONE ) THEN ; nldf_dyn = np_no_ldf ; ioptio = ioptio + 1 ; ENDIF 162 IF( ln_dynldf_lap ) THEN ; ioptio = ioptio + 1 ; ENDIF 163 IF( ln_dynldf_blp ) THEN ; ioptio = ioptio + 1 ; ENDIF 164 IF( ioptio /= 1 ) CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' ) 165 ! 166 IF(.NOT.ln_dynldf_NONE ) THEN !== direction ==>> type of operator ==! 167 ioptio = 0 168 IF( ln_dynldf_lev ) ioptio = ioptio + 1 169 IF( ln_dynldf_hor ) ioptio = ioptio + 1 170 IF( ln_dynldf_iso ) ioptio = ioptio + 1 171 IF( ioptio /= 1 ) CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 direction options (level/hor/iso)' ) 172 ! 173 ! ! Set nldf_dyn, the type of lateral diffusion, from ln_dynldf_... logicals 174 ierr = 0 175 IF( ln_dynldf_lap ) THEN ! laplacian operator 176 IF( ln_zco ) THEN ! z-coordinate 177 IF ( ln_dynldf_lev ) nldf_dyn = np_lap ! iso-level = horizontal (no rotation) 178 IF ( ln_dynldf_hor ) nldf_dyn = np_lap ! iso-level = horizontal (no rotation) 179 IF ( ln_dynldf_iso ) nldf_dyn = np_lap_i ! iso-neutral ( rotation) 180 ENDIF 181 IF( ln_zps ) THEN ! z-coordinate with partial step 182 IF ( ln_dynldf_lev ) nldf_dyn = np_lap ! iso-level (no rotation) 183 IF ( ln_dynldf_hor ) nldf_dyn = np_lap ! iso-level (no rotation) 184 IF ( ln_dynldf_iso ) nldf_dyn = np_lap_i ! iso-neutral ( rotation) 185 ENDIF 186 IF( ln_sco ) THEN ! s-coordinate 187 IF ( ln_dynldf_lev ) nldf_dyn = np_lap ! iso-level = horizontal (no rotation) 188 IF ( ln_dynldf_hor ) nldf_dyn = np_lap_i ! horizontal ( rotation) 189 IF ( ln_dynldf_iso ) nldf_dyn = np_lap_i ! iso-neutral ( rotation) 190 ENDIF 191 ENDIF 192 ! 193 IF( ln_dynldf_blp ) THEN ! bilaplacian operator 194 IF( ln_zco ) THEN ! z-coordinate 195 IF( ln_dynldf_lev ) nldf_dyn = np_blp ! iso-level = horizontal (no rotation) 196 IF( ln_dynldf_hor ) nldf_dyn = np_blp ! iso-level = horizontal (no rotation) 197 IF( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 198 ENDIF 199 IF( ln_zps ) THEN ! z-coordinate with partial step 200 IF( ln_dynldf_lev ) nldf_dyn = np_blp ! iso-level (no rotation) 201 IF( ln_dynldf_hor ) nldf_dyn = np_blp ! iso-level (no rotation) 202 IF( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 203 ENDIF 204 IF( ln_sco ) THEN ! s-coordinate 205 IF( ln_dynldf_lev ) nldf_dyn = np_blp ! iso-level (no rotation) 206 IF( ln_dynldf_hor ) ierr = 2 ! horizontal ( rotation) 207 IF( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 208 ENDIF 209 ENDIF 210 ! 211 IF( ierr == 2 ) CALL ctl_stop( 'rotated bi-laplacian operator does not exist' ) 212 ! 213 IF( nldf_dyn == np_lap_i ) l_ldfslp = .TRUE. ! rotation require the computation of the slopes 214 ! 215 ENDIF 216 ! 217 IF(lwp) THEN 218 WRITE(numout,*) 219 SELECT CASE( nldf_dyn ) 220 CASE( np_no_ldf ) ; WRITE(numout,*) ' ==>>> NO lateral viscosity' 221 CASE( np_lap ) ; WRITE(numout,*) ' ==>>> iso-level laplacian operator' 222 CASE( np_lap_i ) ; WRITE(numout,*) ' ==>>> rotated laplacian operator with iso-level background' 223 CASE( np_blp ) ; WRITE(numout,*) ' ==>>> iso-level bi-laplacian operator' 224 END SELECT 225 WRITE(numout,*) 226 ENDIF 227 228 ! 229 ! !== Space/time variation of eddy coefficients ==! 230 ! !=================================================! 231 ! 232 l_ldfdyn_time = .FALSE. ! no time variation except in case defined below 233 ! 143 234 IF( ln_dynldf_NONE ) THEN 144 235 IF(lwp) WRITE(numout,*) ' ==>>> No viscous operator selected. ahmt and ahmf are not allocated' 145 l_ldfdyn_time = .FALSE.146 236 RETURN 147 ENDIF 148 ! 149 IF( ln_dynldf_blp .AND. ln_dynldf_iso ) THEN ! iso-neutral bilaplacian not implemented 150 CALL ctl_stop( 'dyn_ldf_init: iso-neutral bilaplacian not coded yet' ) 151 ENDIF 152 153 ! ... Space/Time variation of eddy coefficients 154 ! ! allocate the ahm arrays 155 ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr ) 156 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays') 157 ! 158 ahmt(:,:,jpk) = 0._wp ! last level always 0 159 ahmf(:,:,jpk) = 0._wp 160 ! 161 ! ! value of eddy mixing coef. 162 IF ( ln_dynldf_lap ) THEN ; zah0 = rn_ahm_0 ! laplacian operator 163 ELSEIF( ln_dynldf_blp ) THEN ; zah0 = ABS( rn_bhm_0 ) ! bilaplacian operator 164 ELSE ! NO viscous operator 165 CALL ctl_warn( 'ldf_dyn_init: No lateral viscous operator used ' ) 166 ENDIF 167 ! 168 l_ldfdyn_time = .FALSE. ! no time variation except in case defined below 169 ! 170 IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN ! only if a lateral diffusion operator is used 171 ! 172 SELECT CASE( nn_ahm_ijk_t ) ! Specification of space time variations of ahmt, ahmf 237 ! 238 ELSE !== a lateral diffusion operator is used ==! 239 ! 240 ! ! allocate the ahm arrays 241 ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr ) 242 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays') 243 ! 244 ahmt(:,:,jpk) = 0._wp ! last level always 0 245 ahmf(:,:,jpk) = 0._wp 246 ! 247 ! ! value of lap/blp eddy mixing coef. 248 IF( ln_dynldf_lap ) THEN ; zUfac = r1_2 *rn_Uv ; inn = 1 ; cl_Units = ' m2/s' ! laplacian 249 ELSEIF( ln_dynldf_blp ) THEN ; zUfac = r1_12*rn_Uv ; inn = 3 ; cl_Units = ' m4/s' ! bilaplacian 250 ENDIF 251 zah0 = zUfac * rn_Lv**inn ! mixing coefficient 252 zah_max = zUfac * (ra*rad)**inn ! maximum reachable coefficient (value at the Equator) 253 ! 254 SELECT CASE( nn_ahm_ijk_t ) !* Specification of space-time variations of ahmt, ahmf 173 255 ! 174 256 CASE( 0 ) !== constant ==! 175 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef. = constant '176 ahmt(:,:, :) = zah0 * tmask(:,:,:)177 ahmf(:,:, :) = zah0 * fmask(:,:,:)257 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity. = constant = ', zah0, cl_Units 258 ahmt(:,:,1:jpkm1) = zah0 259 ahmf(:,:,1:jpkm1) = zah0 178 260 ! 179 261 CASE( 10 ) !== fixed profile ==! 180 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef. = F( depth )' 181 ahmt(:,:,1) = zah0 * tmask(:,:,1) ! constant surface value 182 ahmf(:,:,1) = zah0 * fmask(:,:,1) 183 CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 262 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F( depth )' 263 IF(lwp) WRITE(numout,*) ' surface viscous coef. = constant = ', zah0, cl_Units 264 ahmt(:,:,1) = zah0 ! constant surface value 265 ahmf(:,:,1) = zah0 266 CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 184 267 ! 185 268 CASE ( -20 ) !== fixed horizontal shape read in file ==! 186 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef.= F(i,j) read in eddy_viscosity.nc file'269 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F(i,j) read in eddy_viscosity.nc file' 187 270 CALL iom_open( 'eddy_viscosity_2D.nc', inum ) 188 271 CALL iom_get ( inum, jpdom_data, 'ahmt_2d', ahmt(:,:,1) ) 189 272 CALL iom_get ( inum, jpdom_data, 'ahmf_2d', ahmf(:,:,1) ) 190 273 CALL iom_close( inum ) 191 !!gm Question : info for LAP or BLP case to take into account the SQRT in the bilaplacian case ???192 !! do we introduce a scaling by the max value of the array, and then multiply by zah0 ????193 !! better: check that the max is <=1 i.e. it is a shape from 0 to 1, not a coef that has physical dimension194 274 DO jk = 2, jpkm1 195 ahmt(:,:,jk) = ahmt(:,:,1) * tmask(:,:,jk)196 ahmf(:,:,jk) = ahmf(:,:,1) * fmask(:,:,jk)275 ahmt(:,:,jk) = ahmt(:,:,1) 276 ahmf(:,:,jk) = ahmf(:,:,1) 197 277 END DO 198 278 ! 199 279 CASE( 20 ) !== fixed horizontal shape ==! 200 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)' 201 IF( ln_dynldf_lap ) CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor 202 IF( ln_dynldf_blp ) CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor^3 280 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)' 281 IF(lwp) WRITE(numout,*) ' using a fixed viscous velocity = ', rn_Uv ,' m/s and Lv = Max(e1,e2)' 282 IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, ' for e1=1°)' 283 CALL ldf_c2d( 'DYN', zUfac , inn , ahmt, ahmf ) ! surface value proportional to scale factor^inn 203 284 ! 204 285 CASE( -30 ) !== fixed 3D shape read in file ==! 205 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file'286 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F(i,j,k) read in eddy_viscosity_3D.nc file' 206 287 CALL iom_open( 'eddy_viscosity_3D.nc', inum ) 207 288 CALL iom_get ( inum, jpdom_data, 'ahmt_3d', ahmt ) 208 289 CALL iom_get ( inum, jpdom_data, 'ahmf_3d', ahmf ) 209 290 CALL iom_close( inum ) 210 !!gm Question : info for LAP or BLP case to take into account the SQRT in the bilaplacian case ????211 !! do we introduce a scaling by the max value of the array, and then multiply by zah0 ????212 DO jk = 1, jpkm1213 ahmt(:,:,jk) = ahmt(:,:,jk) * tmask(:,:,jk)214 ahmf(:,:,jk) = ahmf(:,:,jk) * fmask(:,:,jk)215 END DO216 291 ! 217 292 CASE( 30 ) !== fixed 3D shape ==! 218 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef.= F( latitude, longitude, depth )'219 IF( ln_dynldf_lap ) CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor220 IF( ln_dynldf_blp ) CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor221 ! ! reduction with depth222 CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )293 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F( latitude, longitude, depth )' 294 IF(lwp) WRITE(numout,*) ' using a fixed viscous velocity = ', rn_Uv ,' m/s and Ld = Max(e1,e2)' 295 IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, ' for e1=1°)' 296 CALL ldf_c2d( 'DYN', zUfac , inn , ahmt, ahmf ) ! surface value proportional to scale factor^inn 297 CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) ! reduction with depth 223 298 ! 224 299 CASE( 31 ) !== time varying 3D field ==! 225 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef.= F( latitude, longitude, depth , time )'226 IF(lwp) WRITE(numout,*) ' proportional to the velocity : |u|e/12 or |u|e^3/12'300 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F( latitude, longitude, depth , time )' 301 IF(lwp) WRITE(numout,*) ' proportional to the local velocity : 1/2 |u|e (lap) or 1/12 |u|e^3 (blp)' 227 302 ! 228 303 l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90 229 304 ! 230 305 CASE( 32 ) !== time varying 3D field ==! 231 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef. = F( latitude, longitude, depth , time )' 232 IF(lwp) WRITE(numout,*) ' proportional to the local deformation rate and gridscale (Smagorinsky)' 233 IF(lwp) WRITE(numout,*) ' : L^2|D| or L^4|D|/8' 306 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F( latitude, longitude, depth , time )' 307 IF(lwp) WRITE(numout,*) ' proportional to the local deformation rate and gridscale (Smagorinsky)' 234 308 ! 235 309 l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90 236 310 ! 237 ! allocate arrays used in ldf_dyn.311 ! ! allocate arrays used in ldf_dyn. 238 312 ALLOCATE( dtensq(jpi,jpj) , dshesq(jpi,jpj) , esqt(jpi,jpj) , esqf(jpi,jpj) , STAT=ierr ) 239 313 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays') 240 314 ! 241 ! Set local gridscale values 242 DO jj = 2, jpjm1 315 DO jj = 2, jpjm1 ! Set local gridscale values 243 316 DO ji = fs_2, fs_jpim1 244 317 esqt(ji,jj) = ( e1e2t(ji,jj) /( e1t(ji,jj) + e2t(ji,jj) ) )**2 … … 251 324 END SELECT 252 325 ! 253 IF( ln_dynldf_blp .AND. .NOT. l_ldfdyn_time ) THEN ! bilapcian and no time variation: 254 ahmt(:,:,:) = SQRT( ahmt(:,:,:) ) ! take the square root of the coefficient 255 ahmf(:,:,:) = SQRT( ahmf(:,:,:) ) 326 IF( .NOT.l_ldfdyn_time ) THEN !* No time variation 327 IF( ln_dynldf_lap ) THEN ! laplacian operator (mask only) 328 ahmt(:,:,1:jpkm1) = ahmt(:,:,1:jpkm1) * tmask(:,:,1:jpkm1) 329 ahmf(:,:,1:jpkm1) = ahmf(:,:,1:jpkm1) * fmask(:,:,1:jpkm1) 330 ELSEIF( ln_dynldf_blp ) THEN ! bilaplacian operator (square root + mask) 331 ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1) 332 ahmf(:,:,1:jpkm1) = SQRT( ahmf(:,:,1:jpkm1) ) * fmask(:,:,1:jpkm1) 333 ENDIF 256 334 ENDIF 257 335 ! … … 341 419 & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) & 342 420 & * r1_e2t(ji,jj) * e1t(ji,jj) ) * tmask(ji,jj,jk) 343 dtensq(ji,jj) = zdb *zdb421 dtensq(ji,jj) = zdb * zdb 344 422 END DO 345 423 END DO … … 351 429 & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) & 352 430 & * r1_e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,jk) 353 dshesq(ji,jj) = zdb *zdb431 dshesq(ji,jj) = zdb * zdb 354 432 END DO 355 433 END DO … … 385 463 ! 386 464 IF( ln_dynldf_blp ) THEN ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 |D|/8) 387 ! = sqrt( A_lap_smag L^2/8 ) 388 ! stability limits already applied to laplacian values 389 ! effective default limits are |U|L^3/12 < B_hm < L^4/(32*2dt) 390 ! 465 ! ! = sqrt( A_lap_smag L^2/8 ) 466 ! ! stability limits already applied to laplacian values 467 ! ! effective default limits are 1/12 |U|L^3 < B_hm < 1//(32*2dt) L^4 391 468 DO jk = 1, jpkm1 392 !393 469 DO jj = 2, jpjm1 394 470 DO ji = fs_2, fs_jpim1 395 ! 396 ahmt(ji,jj,jk) = sqrt( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 397 ahmf(ji,jj,jk) = sqrt( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 398 ! 471 ahmt(ji,jj,jk) = SQRT( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 472 ahmf(ji,jj,jk) = SQRT( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 399 473 END DO 400 474 END DO -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r9124 r9490 22 22 USE oce ! ocean dynamics and tracers 23 23 USE dom_oce ! ocean space and time domain 24 USE ldfdyn ! lateral diffusion: eddy viscosity coef.24 ! USE ldfdyn ! lateral diffusion: eddy viscosity coef. 25 25 USE phycst ! physical constants 26 26 USE zdfmxl ! mixed layer depth … … 43 43 LOGICAL , PUBLIC :: l_ldfslp = .FALSE. !: slopes flag 44 44 45 LOGICAL , PUBLIC :: ln_traldf_iso = .TRUE. !: iso-neutral direction 46 LOGICAL , PUBLIC :: ln_traldf_triad = .FALSE. !: griffies triad scheme 47 48 LOGICAL , PUBLIC :: ln_triad_iso = .FALSE. !: pure horizontal mixing in ML 49 LOGICAL , PUBLIC :: ln_botmix_triad = .FALSE. !: mixing on bottom 50 REAL(wp), PUBLIC :: rn_sw_triad = 1._wp !: =1 switching triads ; =0 all four triads used 51 REAL(wp), PUBLIC :: rn_slpmax = 0.01_wp !: slope limit 45 LOGICAL , PUBLIC :: ln_traldf_iso = .TRUE. !: iso-neutral direction (nam_traldf namelist) 46 LOGICAL , PUBLIC :: ln_traldf_triad = .FALSE. !: griffies triad scheme (nam_traldf namelist) 47 LOGICAL , PUBLIC :: ln_dynldf_iso !: iso-neutral direction (nam_dynldf namelist) 48 49 LOGICAL , PUBLIC :: ln_triad_iso = .FALSE. !: pure horizontal mixing in ML (nam_traldf namelist) 50 LOGICAL , PUBLIC :: ln_botmix_triad = .FALSE. !: mixing on bottom (nam_traldf namelist) 51 REAL(wp), PUBLIC :: rn_sw_triad = 1._wp !: =1 switching triads ; =0 all four triads used (nam_traldf namelist) 52 REAL(wp), PUBLIC :: rn_slpmax = 0.01_wp !: slope limit (nam_traldf namelist) 52 53 53 54 LOGICAL , PUBLIC :: l_grad_zps = .FALSE. !: special treatment for Horz Tgradients w partial steps (triad operator) … … 749 750 ! 750 751 IF( ln_traldf_triad ) THEN ! Griffies operator : triad of slopes 751 IF(lwp) WRITE(numout,*) ' Griffies (triad) operator initialisation'752 IF(lwp) WRITE(numout,*) ' ==>>> triad) operator (Griffies)' 752 753 ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , & 753 754 & triadi (jpi,jpj,jpk,0:1,0:1) , triadj (jpi,jpj,jpk,0:1,0:1) , & … … 757 758 ! 758 759 ELSE ! Madec operator : slopes at u-, v-, and w-points 759 IF(lwp) WRITE(numout,*) ' Madec operator initialisation'760 IF(lwp) WRITE(numout,*) ' ==>>> iso operator (Madec)' 760 761 ALLOCATE( omlmask(jpi,jpj,jpk) , & 761 762 & uslp(jpi,jpj,jpk) , uslpml(jpi,jpj) , wslpi(jpi,jpj,jpk) , wslpiml(jpi,jpj) , & -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90
r9169 r9490 50 50 LOGICAL , PUBLIC :: ln_traldf_hor !: horizontal (geopotential) direction 51 51 ! LOGICAL , PUBLIC :: ln_traldf_iso !: iso-neutral direction (see ldfslp) 52 ! != iso-neutral options =! 52 53 ! LOGICAL , PUBLIC :: ln_traldf_triad !: griffies triad scheme (see ldfslp) 53 54 LOGICAL , PUBLIC :: ln_traldf_msc !: Method of Stabilizing Correction … … 58 59 ! != Coefficients =! 59 60 INTEGER , PUBLIC :: nn_aht_ijk_t !: choice of time & space variations of the lateral eddy diffusivity coef. 60 REAL(wp), PUBLIC :: rn_aht_0 !: laplacian lateral eddy diffusivity [m2/s] 61 REAL(wp), PUBLIC :: rn_bht_0 !: bilaplacian lateral eddy diffusivity [m4/s] 62 63 ! !!* Namelist namtra_ldfeiv : eddy induced velocity param. * 61 ! ! time invariant coefficients: aht_0 = 1/2 Ud*Ld (lap case) 62 ! ! bht_0 = 1/12 Ud*Ld^3 (blp case) 63 REAL(wp), PUBLIC :: rn_Ud !: lateral diffusive velocity [m/s] 64 REAL(wp), PUBLIC :: rn_Ld !: lateral diffusive length [m] 65 66 ! !!* Namelist namtra_eiv : eddy induced velocity param. * 64 67 ! != Use/diagnose eiv =! 65 68 LOGICAL , PUBLIC :: ln_ldfeiv !: eddy induced velocity flag … … 67 70 ! != Coefficients =! 68 71 INTEGER , PUBLIC :: nn_aei_ijk_t !: choice of time/space variation of the eiv coeff. 69 REAL(wp), PUBLIC :: rn_aeiv_0 !: eddy induced velocity coefficient [m2/s] 72 REAL(wp), PUBLIC :: rn_Ue !: lateral diffusive velocity [m/s] 73 REAL(wp), PUBLIC :: rn_Le !: lateral diffusive length [m] 70 74 75 ! ! Flag to control the type of lateral diffusive operator 76 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 ! error in specification of lateral diffusion 77 INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 ! without operator (i.e. no lateral diffusive trend) 78 ! !! laplacian ! bilaplacian ! 79 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 ! iso-level operator 80 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 , np_blp_i = 21 ! standard iso-neutral or geopotential operator 81 INTEGER, PARAMETER, PUBLIC :: np_lap_it = 12 , np_blp_it = 22 ! triad iso-neutral or geopotential operator 82 83 INTEGER , PUBLIC :: nldf_tra = 0 !: type of lateral diffusion used defined from ln_traldf_... (namlist logicals) 71 84 LOGICAL , PUBLIC :: l_ldftra_time = .FALSE. !: flag for time variation of the lateral eddy diffusivity coef. 72 LOGICAL , PUBLIC :: l_ldfeiv_time = .FALSE. ! flag for time variation of the eiv coef.85 LOGICAL , PUBLIC :: l_ldfeiv_time = .FALSE. !: flag for time variation of the eiv coef. 73 86 74 87 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahtu, ahtv !: eddy diffusivity coef. at U- and V-points [m2/s] 75 88 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aeiu, aeiv !: eddy induced velocity coeff. [m2/s] 76 89 90 REAL(wp) :: aht0, aei0 ! constant eddy coefficients (deduced from namelist values) [m2/s] 91 REAL(wp) :: r1_2 = 0.5_wp ! =1/2 77 92 REAL(wp) :: r1_4 = 0.25_wp ! =1/4 78 93 REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 … … 108 123 !! =-30 => = F(i,j,k) = shape read in 'eddy_diffusivity.nc' file 109 124 !! = 30 = F(i,j,k) = 2D (case 20) + decrease with depth (case 10) 110 !! = 31 = F(i,j,k,t) = F(local velocity) ( |u|e /12laplacian operator111 !! or |u|e^3/12bilaplacian operator )125 !! = 31 = F(i,j,k,t) = F(local velocity) ( 1/2 |u|e laplacian operator 126 !! or 1/12 |u|e^3 bilaplacian operator ) 112 127 !! * initialisation of the eddy induced velocity coefficient by a call to ldf_eiv_init 113 128 !! 114 !! ** action : ahtu, ahtv initialized once for all or l_ldftra_time set to true 115 !! aeiu, aeiv initialized once for all or l_ldfeiv_time set to true 116 !!---------------------------------------------------------------------- 117 INTEGER :: jk ! dummy loop indices 118 INTEGER :: ierr, inum, ios ! local integer 119 REAL(wp) :: zah0 ! local scalar 129 !! ** action : ahtu, ahtv initialized one for all or l_ldftra_time set to true 130 !! aeiu, aeiv initialized one for all or l_ldfeiv_time set to true 131 !!---------------------------------------------------------------------- 132 INTEGER :: jk ! dummy loop indices 133 INTEGER :: ioptio, ierr, inum, ios, inn ! local integer 134 REAL(wp) :: zah_max, zUfac ! - - 135 CHARACTER(len=5) :: cl_Units ! units (m2/s or m4/s) 120 136 !! 121 137 NAMELIST/namtra_ldf/ ln_traldf_NONE, ln_traldf_lap , ln_traldf_blp , & ! type of operator … … 123 139 & ln_traldf_iso , ln_traldf_msc , rn_slpmax , & ! option for iso-neutral operator 124 140 & ln_triad_iso , ln_botmix_triad, rn_sw_triad , & ! option for triad operator 125 & rn_aht_0 , rn_bht_0 , nn_aht_ijk_t! lateral eddy coefficient141 & nn_aht_ijk_t , rn_Ud , rn_Ld ! lateral eddy coefficient 126 142 !!---------------------------------------------------------------------- 127 143 ! 128 144 IF(lwp) THEN ! control print 129 145 WRITE(numout,*) 130 WRITE(numout,*) 'ldf_tra_init : lateral tracer physics'146 WRITE(numout,*) 'ldf_tra_init : lateral tracer diffusion' 131 147 WRITE(numout,*) '~~~~~~~~~~~~ ' 132 148 ENDIF 149 133 150 ! 134 151 ! Choice of lateral tracer physics … … 154 171 WRITE(numout,*) ' iso-neutral Madec operator ln_traldf_iso = ', ln_traldf_iso 155 172 WRITE(numout,*) ' iso-neutral triad operator ln_traldf_triad = ', ln_traldf_triad 156 WRITE(numout,*) ' iso-neutral (Method of Stab. Corr.)ln_traldf_msc = ', ln_traldf_msc173 WRITE(numout,*) ' use the Method of Stab. Correction ln_traldf_msc = ', ln_traldf_msc 157 174 WRITE(numout,*) ' maximum isoppycnal slope rn_slpmax = ', rn_slpmax 158 175 WRITE(numout,*) ' pure lateral mixing in ML ln_triad_iso = ', ln_triad_iso … … 160 177 WRITE(numout,*) ' lateral mixing on bottom ln_botmix_triad = ', ln_botmix_triad 161 178 WRITE(numout,*) ' coefficients :' 162 WRITE(numout,*) ' lateral eddy diffusivity (lap case) rn_aht_0 = ', rn_aht_0163 WRITE(numout,*) ' lateral eddy diffusivity (bilap case) rn_bht_0 = ', rn_bht_0164 179 WRITE(numout,*) ' type of time-space variation nn_aht_ijk_t = ', nn_aht_ijk_t 165 ENDIF 166 ! 167 ! ! Parameter control 168 ! 169 IF( ln_traldf_NONE ) THEN 170 IF(lwp) WRITE(numout,*) ' ==>>> No diffusive operator selected. ahtu and ahtv are not allocated' 171 l_ldftra_time = .FALSE. 172 RETURN 173 ENDIF 180 WRITE(numout,*) ' lateral diffusive velocity (if cst) rn_Ud = ', rn_Ud, ' m/s' 181 WRITE(numout,*) ' lateral diffusive length (if cst) rn_Ld = ', rn_Ld, ' m' 182 ENDIF 183 ! 184 ! 185 ! Operator and its acting direction (set nldf_tra) 186 ! ================================= 187 ! 188 nldf_tra = np_ERROR 189 ioptio = 0 190 IF( ln_traldf_NONE ) THEN ; nldf_tra = np_no_ldf ; ioptio = ioptio + 1 ; ENDIF 191 IF( ln_traldf_lap ) THEN ; ioptio = ioptio + 1 ; ENDIF 192 IF( ln_traldf_blp ) THEN ; ioptio = ioptio + 1 ; ENDIF 193 IF( ioptio /= 1 ) CALL ctl_stop( 'tra_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' ) 194 ! 195 IF( .NOT.ln_traldf_NONE ) THEN !== direction ==>> type of operator ==! 196 ioptio = 0 197 IF( ln_traldf_lev ) ioptio = ioptio + 1 198 IF( ln_traldf_hor ) ioptio = ioptio + 1 199 IF( ln_traldf_iso ) ioptio = ioptio + 1 200 IF( ioptio /= 1 ) CALL ctl_stop( 'tra_ldf_init: use ONE direction (level/hor/iso)' ) 201 ! 202 ! ! defined the type of lateral diffusion from ln_traldf_... logicals 203 ierr = 0 204 IF ( ln_traldf_lap ) THEN ! laplacian operator 205 IF ( ln_zco ) THEN ! z-coordinate 206 IF ( ln_traldf_lev ) nldf_tra = np_lap ! iso-level = horizontal (no rotation) 207 IF ( ln_traldf_hor ) nldf_tra = np_lap ! iso-level = horizontal (no rotation) 208 IF ( ln_traldf_iso ) nldf_tra = np_lap_i ! iso-neutral: standard ( rotation) 209 IF ( ln_traldf_triad ) nldf_tra = np_lap_it ! iso-neutral: triad ( rotation) 210 ENDIF 211 IF ( ln_zps ) THEN ! z-coordinate with partial step 212 IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed 213 IF ( ln_traldf_hor ) nldf_tra = np_lap ! horizontal (no rotation) 214 IF ( ln_traldf_iso ) nldf_tra = np_lap_i ! iso-neutral: standard (rotation) 215 IF ( ln_traldf_triad ) nldf_tra = np_lap_it ! iso-neutral: triad (rotation) 216 ENDIF 217 IF ( ln_sco ) THEN ! s-coordinate 218 IF ( ln_traldf_lev ) nldf_tra = np_lap ! iso-level (no rotation) 219 IF ( ln_traldf_hor ) nldf_tra = np_lap_i ! horizontal ( rotation) 220 IF ( ln_traldf_iso ) nldf_tra = np_lap_i ! iso-neutral: standard ( rotation) 221 IF ( ln_traldf_triad ) nldf_tra = np_lap_it ! iso-neutral: triad ( rotation) 222 ENDIF 223 ENDIF 224 ! 225 IF( ln_traldf_blp ) THEN ! bilaplacian operator 226 IF ( ln_zco ) THEN ! z-coordinate 227 IF ( ln_traldf_lev ) nldf_tra = np_blp ! iso-level = horizontal (no rotation) 228 IF ( ln_traldf_hor ) nldf_tra = np_blp ! iso-level = horizontal (no rotation) 229 IF ( ln_traldf_iso ) nldf_tra = np_blp_i ! iso-neutral: standard ( rotation) 230 IF ( ln_traldf_triad ) nldf_tra = np_blp_it ! iso-neutral: triad ( rotation) 231 ENDIF 232 IF ( ln_zps ) THEN ! z-coordinate with partial step 233 IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed 234 IF ( ln_traldf_hor ) nldf_tra = np_blp ! horizontal (no rotation) 235 IF ( ln_traldf_iso ) nldf_tra = np_blp_i ! iso-neutral: standard ( rotation) 236 IF ( ln_traldf_triad ) nldf_tra = np_blp_it ! iso-neutral: triad ( rotation) 237 ENDIF 238 IF ( ln_sco ) THEN ! s-coordinate 239 IF ( ln_traldf_lev ) nldf_tra = np_blp ! iso-level (no rotation) 240 IF ( ln_traldf_hor ) nldf_tra = np_blp_it ! horizontal ( rotation) 241 IF ( ln_traldf_iso ) nldf_tra = np_blp_i ! iso-neutral: standard ( rotation) 242 IF ( ln_traldf_triad ) nldf_tra = np_blp_it ! iso-neutral: triad ( rotation) 243 ENDIF 244 ENDIF 245 IF ( ierr == 1 ) CALL ctl_stop( 'iso-level in z-partial step, not allowed' ) 246 ENDIF 247 ! 248 IF( ln_ldfeiv .AND. .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) ) & 249 & CALL ctl_stop( 'ln_ldfeiv=T requires iso-neutral laplacian diffusion' ) 250 IF( ln_isfcav .AND. ln_traldf_triad ) & 251 & CALL ctl_stop( ' ice shelf cavity and traldf_triad not tested' ) 252 ! 253 IF( nldf_tra == np_lap_i .OR. nldf_tra == np_lap_it .OR. & 254 & nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it ) l_ldfslp = .TRUE. ! slope of neutral surfaces required 174 255 ! 175 256 IF( ln_traldf_blp .AND. ( ln_traldf_iso .OR. ln_traldf_triad) ) THEN ! iso-neutral bilaplacian need MSC … … 177 258 ENDIF 178 259 ! 260 IF(lwp) THEN 261 WRITE(numout,*) 262 SELECT CASE( nldf_tra ) 263 CASE( np_no_ldf ) ; WRITE(numout,*) ' ==>>> NO lateral diffusion' 264 CASE( np_lap ) ; WRITE(numout,*) ' ==>>> laplacian iso-level operator' 265 CASE( np_lap_i ) ; WRITE(numout,*) ' ==>>> Rotated laplacian operator (standard)' 266 CASE( np_lap_it ) ; WRITE(numout,*) ' ==>>> Rotated laplacian operator (triad)' 267 CASE( np_blp ) ; WRITE(numout,*) ' ==>>> bilaplacian iso-level operator' 268 CASE( np_blp_i ) ; WRITE(numout,*) ' ==>>> Rotated bilaplacian operator (standard)' 269 CASE( np_blp_it ) ; WRITE(numout,*) ' ==>>> Rotated bilaplacian operator (triad)' 270 END SELECT 271 WRITE(numout,*) 272 ENDIF 273 274 ! 179 275 ! Space/time variation of eddy coefficients 180 276 ! =========================================== 181 ! ! allocate the aht arrays 182 ALLOCATE( ahtu(jpi,jpj,jpk) , ahtv(jpi,jpj,jpk) , STAT=ierr ) 183 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_tra_init: failed to allocate arrays') 184 ! 185 ahtu(:,:,jpk) = 0._wp ! last level always 0 186 ahtv(:,:,jpk) = 0._wp 187 ! 188 ! ! value of eddy mixing coef. 189 IF ( ln_traldf_lap ) THEN ; zah0 = rn_aht_0 ! laplacian operator 190 ELSEIF( ln_traldf_blp ) THEN ; zah0 = ABS( rn_bht_0 ) ! bilaplacian operator 191 ENDIF 192 ! 193 l_ldftra_time = .FALSE. ! no time variation except in case defined below 194 ! 195 IF( ln_traldf_lap .OR. ln_traldf_blp ) THEN ! only if a lateral diffusion operator is used 196 ! 197 SELECT CASE( nn_aht_ijk_t ) ! Specification of space time variations of ehtu, ahtv 277 ! 278 l_ldftra_time = .FALSE. ! no time variation except in case defined below 279 ! 280 IF( ln_traldf_NONE ) THEN !== no explicit diffusive operator ==! 281 ! 282 IF(lwp) WRITE(numout,*) ' ==>>> No diffusive operator selected. ahtu and ahtv are not allocated' 283 RETURN 284 ! 285 ELSE !== a lateral diffusion operator is used ==! 286 ! 287 ! ! allocate the aht arrays 288 ALLOCATE( ahtu(jpi,jpj,jpk) , ahtv(jpi,jpj,jpk) , STAT=ierr ) 289 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_tra_init: failed to allocate arrays') 290 ! 291 ahtu(:,:,jpk) = 0._wp ! last level always 0 292 ahtv(:,:,jpk) = 0._wp 293 !. 294 ! ! value of lap/blp eddy mixing coef. 295 IF( ln_traldf_lap ) THEN ; zUfac = r1_2 *rn_Ud ; inn = 1 ; cl_Units = ' m2/s' ! laplacian 296 ELSEIF( ln_traldf_blp ) THEN ; zUfac = r1_12*rn_Ud ; inn = 3 ; cl_Units = ' m4/s' ! bilaplacian 297 ENDIF 298 aht0 = zUfac * rn_Ld**inn ! mixing coefficient 299 zah_max = zUfac * (ra*rad)**inn ! maximum reachable coefficient (value at the Equator for e1=1 degree) 300 ! 301 ! 302 SELECT CASE( nn_aht_ijk_t ) !* Specification of space-time variations of ahtu, ahtv 198 303 ! 199 304 CASE( 0 ) !== constant ==! 200 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = constant = ', rn_aht_0201 ahtu(:,:, :) = zah0 * umask(:,:,:)202 ahtv(:,:, :) = zah0 * vmask(:,:,:)305 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = constant = ', aht0, cl_Units 306 ahtu(:,:,1:jpkm1) = aht0 307 ahtv(:,:,1:jpkm1) = aht0 203 308 ! 204 309 CASE( 10 ) !== fixed profile ==! 205 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( depth )' 206 ahtu(:,:,1) = zah0 * umask(:,:,1) ! constant surface value 207 ahtv(:,:,1) = zah0 * vmask(:,:,1) 208 CALL ldf_c1d( 'TRA', r1_4, ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv ) 209 ! 210 CASE ( -20 ) !== fixed horizontal shape read in file ==! 211 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F(i,j) read in eddy_diffusivity.nc file' 310 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( depth )' 311 IF(lwp) WRITE(numout,*) ' surface eddy diffusivity = constant = ', aht0, cl_Units 312 ahtu(:,:,1) = aht0 ! constant surface value 313 ahtv(:,:,1) = aht0 314 CALL ldf_c1d( 'TRA', ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv ) 315 ! 316 CASE ( -20 ) !== fixed horizontal shape and magnitude read in file ==! 317 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F(i,j) read in eddy_diffusivity.nc file' 212 318 CALL iom_open( 'eddy_diffusivity_2D.nc', inum ) 213 319 CALL iom_get ( inum, jpdom_data, 'ahtu_2D', ahtu(:,:,1) ) … … 215 321 CALL iom_close( inum ) 216 322 DO jk = 2, jpkm1 217 ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk)218 ahtv(:,:,jk) = ahtv(:,:,1) * vmask(:,:,jk)323 ahtu(:,:,jk) = ahtu(:,:,1) 324 ahtv(:,:,jk) = ahtv(:,:,1) 219 325 END DO 220 326 ! 221 327 CASE( 20 ) !== fixed horizontal shape ==! 222 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or blp case)' 223 IF( ln_traldf_lap ) CALL ldf_c2d( 'TRA', 'LAP', zah0, ahtu, ahtv ) ! surface value proportional to scale factor 224 IF( ln_traldf_blp ) CALL ldf_c2d( 'TRA', 'BLP', zah0, ahtu, ahtv ) ! surface value proportional to scale factor 328 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or blp case)' 329 IF(lwp) WRITE(numout,*) ' using a fixed diffusive velocity = ', rn_Ud,' m/s and Ld = Max(e1,e2)' 330 IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, ' for e1=1°)' 331 CALL ldf_c2d( 'TRA', zUfac , inn , ahtu, ahtv ) ! value proportional to scale factor^inn 225 332 ! 226 333 CASE( 21 ) !== time varying 2D field ==! 227 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( latitude, longitude, time )' 228 IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )' 229 IF(lwp) WRITE(numout,*) ' min value = 0.1 * rn_aht_0' 230 IF(lwp) WRITE(numout,*) ' max value = rn_aht_0 (rn_aeiv_0 if nn_aei_ijk_t=21)' 231 IF(lwp) WRITE(numout,*) ' increased to rn_aht_0 within 20N-20S' 334 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( latitude, longitude, time )' 335 IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )' 336 IF(lwp) WRITE(numout,*) ' min value = 0.2 * aht0 (with aht0= 1/2 rn_Ud*rn_Ld)' 337 IF(lwp) WRITE(numout,*) ' max value = aei0 (with aei0=1/2 rn_Ue*Le increased to aht0 within 20N-20S' 232 338 ! 233 339 l_ldftra_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90 234 340 ! 235 IF( ln_traldf_blp ) THEN 236 CALL ctl_stop( 'ldf_tra_init: aht=F(growth rate of baroc. insta.) incompatible with bilaplacian operator' ) 237 ENDIF 341 IF( ln_traldf_blp ) CALL ctl_stop( 'ldf_tra_init: aht=F( growth rate of baroc. insta .)', & 342 & ' incompatible with bilaplacian operator' ) 238 343 ! 239 344 CASE( -30 ) !== fixed 3D shape read in file ==! 240 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef.= F(i,j,k) read in eddy_diffusivity.nc file'345 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F(i,j,k) read in eddy_diffusivity.nc file' 241 346 CALL iom_open( 'eddy_diffusivity_3D.nc', inum ) 242 347 CALL iom_get ( inum, jpdom_data, 'ahtu_3D', ahtu ) 243 348 CALL iom_get ( inum, jpdom_data, 'ahtv_3D', ahtv ) 244 349 CALL iom_close( inum ) 245 DO jk = 1, jpkm1246 ahtu(:,:,jk) = ahtu(:,:,jk) * umask(:,:,jk)247 ahtv(:,:,jk) = ahtv(:,:,jk) * vmask(:,:,jk)248 END DO249 350 ! 250 351 CASE( 30 ) !== fixed 3D shape ==! 251 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef.= F( latitude, longitude, depth )'252 IF( ln_traldf_lap ) CALL ldf_c2d( 'TRA', 'LAP', zah0, ahtu, ahtv ) ! surface value proportional to scale factor253 IF( ln_traldf_blp ) CALL ldf_c2d( 'TRA', 'BLP', zah0, ahtu, ahtv ) ! surface value proportional to scale factor254 ! ! reduction with depth255 CALL ldf_c1d( 'TRA', r1_4, ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv )352 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( latitude, longitude, depth )' 353 IF(lwp) WRITE(numout,*) ' using a fixed diffusive velocity = ', rn_Ud,' m/s and Ld = Max(e1,e2)' 354 IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, ' for e1=1°)' 355 CALL ldf_c2d( 'TRA', zUfac , inn , ahtu, ahtv ) ! surface value proportional to scale factor^inn 356 CALL ldf_c1d( 'TRA', ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv ) ! reduction with depth 256 357 ! 257 358 CASE( 31 ) !== time varying 3D field ==! 258 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef.= F( latitude, longitude, depth , time )'259 IF(lwp) WRITE(numout,*) ' proportional to the velocity : |u|e/12 or |u|e^3/12'359 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( latitude, longitude, depth , time )' 360 IF(lwp) WRITE(numout,*) ' proportional to the velocity : 1/2 |u|e or 1/12 |u|e^3' 260 361 ! 261 362 l_ldftra_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90 … … 265 366 END SELECT 266 367 ! 267 IF( ln_traldf_blp .AND. .NOT. l_ldftra_time ) THEN 268 ahtu(:,:,:) = SQRT( ahtu(:,:,:) ) 269 ahtv(:,:,:) = SQRT( ahtv(:,:,:) ) 368 IF( .NOT.l_ldftra_time ) THEN !* No time variation 369 IF( ln_traldf_lap ) THEN ! laplacian operator (mask only) 370 ahtu(:,:,1:jpkm1) = ahtu(:,:,1:jpkm1) * umask(:,:,1:jpkm1) 371 ahtv(:,:,1:jpkm1) = ahtv(:,:,1:jpkm1) * vmask(:,:,1:jpkm1) 372 ELSEIF( ln_traldf_blp ) THEN ! bilaplacian operator (square root + mask) 373 ahtu(:,:,1:jpkm1) = SQRT( ahtu(:,:,1:jpkm1) ) * umask(:,:,1:jpkm1) 374 ahtv(:,:,1:jpkm1) = SQRT( ahtv(:,:,1:jpkm1) ) * vmask(:,:,1:jpkm1) 375 ENDIF 270 376 ENDIF 271 377 ! … … 281 387 !! ** Purpose : update at kt the tracer lateral mixing coeff. (aht and aeiv) 282 388 !! 283 !! ** Method : 389 !! ** Method : * time varying eddy diffusivity coefficients: 284 390 !! 285 391 !! nn_aei_ijk_t = 21 aeiu, aeiv = F(i,j, t) = F(growth rate of baroclinic instability) … … 290 396 !! or |u|e^3/12 bilaplacian operator ) 291 397 !! 398 !! * time varying EIV coefficients: call to ldf_eiv routine 399 !! 292 400 !! ** action : ahtu, ahtv update at each time step 293 401 !! aeiu, aeiv - - - - (if ln_ldfeiv=T) … … 296 404 ! 297 405 INTEGER :: ji, jj, jk ! dummy loop indices 298 REAL(wp) :: zaht, zahf, zaht_min, z 1_f20! local scalar406 REAL(wp) :: zaht, zahf, zaht_min, zDaht, z1_f20 ! local scalar 299 407 !!---------------------------------------------------------------------- 300 408 ! 301 409 IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN ! eddy induced velocity coefficients 302 410 ! ! =F(growth rate of baroclinic instability) 303 ! ! max value rn_aeiv_0 ; decreased to 0 within 20N-20S304 CALL ldf_eiv( kt, rn_aeiv_0, aeiu, aeiv )411 ! ! max value aeiv_0 ; decreased to 0 within 20N-20S 412 CALL ldf_eiv( kt, aei0, aeiu, aeiv ) 305 413 ENDIF 306 414 ! … … 308 416 ! 309 417 CASE( 21 ) !== time varying 2D field ==! = F( growth rate of baroclinic instability ) 310 ! ! min value rn_aht_0 / 10311 ! ! max value rn_aht_0 (rn_aeiv_0 if nn_aei_ijk_t=21)312 ! ! increase to rn_aht_0 within 20N-20S418 ! ! min value 0.2*aht0 419 ! ! max value aht0 (aei0 if nn_aei_ijk_t=21) 420 ! ! increase to aht0 within 20N-20S 313 421 IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN ! use the already computed aei. 314 422 ahtu(:,:,1) = aeiu(:,:,1) 315 423 ahtv(:,:,1) = aeiv(:,:,1) 316 424 ELSE ! compute aht. 317 CALL ldf_eiv( kt, rn_aht_0, ahtu, ahtv )425 CALL ldf_eiv( kt, aht0, ahtu, ahtv ) 318 426 ENDIF 319 427 ! 320 z1_f20 = 1._wp / ( 2._wp * omega * SIN( rad * 20._wp ) ) ! 1 / ff(20 degrees) 321 zaht_min = 0.2_wp * rn_aht_0 ! minimum value for aht 428 z1_f20 = 1._wp / ( 2._wp * omega * SIN( rad * 20._wp ) ) ! 1 / ff(20 degrees) 429 zaht_min = 0.2_wp * aht0 ! minimum value for aht 430 zDaht = aht0 - zaht_min 322 431 DO jj = 1, jpj 323 432 DO ji = 1, jpi 324 433 !!gm CAUTION : here we assume lat/lon grid in 20deg N/S band (like all ORCA cfg) 325 434 !! ==>>> The Coriolis value is identical for t- & u_points, and for v- and f-points 326 zaht = ( 1._wp - MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * ( rn_aht_0 - zaht_min )327 zahf = ( 1._wp - MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * ( rn_aht_0 - zaht_min )328 ahtu(ji,jj,1) = ( MAX( zaht_min, ahtu(ji,jj,1) ) + zaht ) * umask(ji,jj,1)! min value zaht_min329 ahtv(ji,jj,1) = ( MAX( zaht_min, ahtv(ji,jj,1) ) + zahf ) * vmask(ji,jj,1)! increase within 20S-20N330 END DO 331 END DO 332 DO jk = 2, jpkm1 ! deeper value = surface value435 zaht = ( 1._wp - MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * zDaht 436 zahf = ( 1._wp - MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * zDaht 437 ahtu(ji,jj,1) = ( MAX( zaht_min, ahtu(ji,jj,1) ) + zaht ) ! min value zaht_min 438 ahtv(ji,jj,1) = ( MAX( zaht_min, ahtv(ji,jj,1) ) + zahf ) ! increase within 20S-20N 439 END DO 440 END DO 441 DO jk = 1, jpkm1 ! deeper value = surface value + mask for all levels 333 442 ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk) 334 443 ahtv(:,:,jk) = ahtv(:,:,1) * vmask(:,:,jk) … … 338 447 IF( ln_traldf_lap ) THEN ! laplacian operator |u| e /12 339 448 DO jk = 1, jpkm1 340 ahtu(:,:,jk) = ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12 449 ahtu(:,:,jk) = ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12 ! n.b. ub,vb are masked 341 450 ahtv(:,:,jk) = ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12 342 451 END DO … … 355 464 CALL iom_put( "ahtv_3d", ahtv(:,:,:) ) ! 3D v-eddy diffusivity coeff. 356 465 ! 357 !!gm : THE IF below is to be checked (comes from Seb)358 466 IF( ln_ldfeiv ) THEN 359 467 CALL iom_put( "aeiu_2d", aeiu(:,:,1) ) ! surface u-EIV coeff. … … 372 480 !! ** Purpose : initialization of the eiv coeff. from namelist choices. 373 481 !! 374 !! ** Method : 375 !! 376 !! ** Action : aeiu , aeiv : EIV coeff. at u- & v-points 377 !! l_ldfeiv_time : =T if EIV coefficients vary with time 378 !!---------------------------------------------------------------------- 379 INTEGER :: jk ! dummy loop indices 380 INTEGER :: ierr, inum, ios ! local integer 381 ! 382 NAMELIST/namtra_ldfeiv/ ln_ldfeiv , ln_ldfeiv_dia, & ! eddy induced velocity (eiv) 383 & nn_aei_ijk_t, rn_aeiv_0 ! eiv coefficient 482 !! ** Method : the eiv diffusivity coef. specification depends on: 483 !! nn_aei_ijk_t = 0 => = constant 484 !! ! 485 !! = 10 => = F(z) : constant with a reduction of 1/4 with depth 486 !! ! 487 !! =-20 => = F(i,j) = shape read in 'eddy_diffusivity.nc' file 488 !! = 20 = F(i,j) = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case) 489 !! = 21 = F(i,j,t) = F(growth rate of baroclinic instability) 490 !! ! 491 !! =-30 => = F(i,j,k) = shape read in 'eddy_diffusivity.nc' file 492 !! = 30 = F(i,j,k) = 2D (case 20) + decrease with depth (case 10) 493 !! 494 !! ** Action : aeiu , aeiv : initialized one for all or l_ldftra_time set to true 495 !! l_ldfeiv_time : =T if EIV coefficients vary with time 496 !!---------------------------------------------------------------------- 497 INTEGER :: jk ! dummy loop indices 498 INTEGER :: ierr, inum, ios, inn ! local integer 499 REAL(wp) :: zah_max, zUfac ! - scalar 500 !! 501 NAMELIST/namtra_eiv/ ln_ldfeiv , ln_ldfeiv_dia, & ! eddy induced velocity (eiv) 502 & nn_aei_ijk_t, rn_Ue, rn_Le ! eiv coefficient 384 503 !!---------------------------------------------------------------------- 385 504 ! … … 390 509 ENDIF 391 510 ! 392 REWIND( numnam_ref ) ! Namelist namtra_ ldfeiv in reference namelist : eddy induced velocity param.393 READ ( numnam_ref, namtra_ ldfeiv, IOSTAT = ios, ERR = 901)394 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_ ldfeiv in reference namelist', lwp )395 ! 396 REWIND( numnam_cfg ) ! Namelist namtra_ ldfeiv in configuration namelist : eddy induced velocity param.397 READ ( numnam_cfg, namtra_ ldfeiv, IOSTAT = ios, ERR = 902 )398 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_ ldfeiv in configuration namelist', lwp )399 IF(lwm) WRITE ( numond, namtra_ ldfeiv )511 REWIND( numnam_ref ) ! Namelist namtra_eiv in reference namelist : eddy induced velocity param. 512 READ ( numnam_ref, namtra_eiv, IOSTAT = ios, ERR = 901) 513 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_eiv in reference namelist', lwp ) 514 ! 515 REWIND( numnam_cfg ) ! Namelist namtra_eiv in configuration namelist : eddy induced velocity param. 516 READ ( numnam_cfg, namtra_eiv, IOSTAT = ios, ERR = 902 ) 517 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_eiv in configuration namelist', lwp ) 518 IF(lwm) WRITE ( numond, namtra_eiv ) 400 519 401 520 IF(lwp) THEN ! control print 402 WRITE(numout,*) ' Namelist namtra_ldfeiv : ' 403 WRITE(numout,*) ' Eddy Induced Velocity (eiv) param. ln_ldfeiv = ', ln_ldfeiv 404 WRITE(numout,*) ' eiv streamfunction & velocity diag. ln_ldfeiv_dia = ', ln_ldfeiv_dia 405 WRITE(numout,*) ' eddy induced velocity coef. rn_aeiv_0 = ', rn_aeiv_0 406 WRITE(numout,*) ' type of time-space variation nn_aei_ijk_t = ', nn_aei_ijk_t 521 WRITE(numout,*) ' Namelist namtra_eiv : ' 522 WRITE(numout,*) ' Eddy Induced Velocity (eiv) param. ln_ldfeiv = ', ln_ldfeiv 523 WRITE(numout,*) ' eiv streamfunction & velocity diag. ln_ldfeiv_dia = ', ln_ldfeiv_dia 524 WRITE(numout,*) ' coefficients :' 525 WRITE(numout,*) ' type of time-space variation nn_aei_ijk_t = ', nn_aht_ijk_t 526 WRITE(numout,*) ' lateral diffusive velocity (if cst) rn_Ue = ', rn_Ue, ' m/s' 527 WRITE(numout,*) ' lateral diffusive length (if cst) rn_Le = ', rn_Le, ' m' 407 528 WRITE(numout,*) 408 529 ENDIF 409 530 ! 410 IF( ln_ldfeiv .AND. ln_traldf_blp ) CALL ctl_stop( 'ldf_eiv_init: eddy induced velocity ONLY with laplacian diffusivity' ) 411 412 ! ! Parameter control 413 l_ldfeiv_time = .FALSE. 414 ! 415 IF( ln_ldfeiv ) THEN ! allocate the aei arrays 531 l_ldfeiv_time = .FALSE. ! no time variation except in case defined below 532 ! 533 ! 534 IF( .NOT.ln_ldfeiv ) THEN !== Parametrization not used ==! 535 ! 536 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity param is NOT used' 537 ln_ldfeiv_dia = .FALSE. 538 ! 539 ELSE !== use the parametrization ==! 540 ! 541 IF(lwp) WRITE(numout,*) ' ==>>> use eddy induced velocity parametrization' 542 IF(lwp) WRITE(numout,*) 543 ! 544 IF( ln_traldf_blp ) CALL ctl_stop( 'ldf_eiv_init: eddy induced velocity ONLY with laplacian diffusivity' ) 545 ! 546 ! != allocate the aei arrays 416 547 ALLOCATE( aeiu(jpi,jpj,jpk), aeiv(jpi,jpj,jpk), STAT=ierr ) 417 548 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'ldf_eiv: failed to allocate arrays') 418 549 ! 419 SELECT CASE( nn_aei_ijk_t ) ! Specification of space time variations of eaiu, aeiv 420 ! 421 CASE( 0 ) !== constant ==! 422 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = constant = ', rn_aeiv_0 423 aeiu(:,:,:) = rn_aeiv_0 424 aeiv(:,:,:) = rn_aeiv_0 425 ! 426 CASE( 10 ) !== fixed profile ==! 550 ! != Specification of space-time variations of eaiu, aeiv 551 ! 552 aeiu(:,:,jpk) = 0._wp ! last level always 0 553 aeiv(:,:,jpk) = 0._wp 554 ! ! value of EIV coef. (laplacian operator) 555 zUfac = r1_2 *rn_Ue ! velocity factor 556 inn = 1 ! L-exponent 557 aei0 = zUfac * rn_Le**inn ! mixing coefficient 558 zah_max = zUfac * (ra*rad)**inn ! maximum reachable coefficient (value at the Equator) 559 560 SELECT CASE( nn_aei_ijk_t ) !* Specification of space-time variations 561 ! 562 CASE( 0 ) !-- constant --! 563 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = constant = ', aei0, ' m2/s' 564 aeiu(:,:,1:jpkm1) = aei0 565 aeiv(:,:,1:jpkm1) = aei0 566 ! 567 CASE( 10 ) !-- fixed profile --! 427 568 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( depth )' 428 aeiu(:,:,1) = rn_aeiv_0 ! constant surface value 429 aeiv(:,:,1) = rn_aeiv_0 430 CALL ldf_c1d( 'TRA', r1_4, aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) 431 ! 432 CASE ( -20 ) !== fixed horizontal shape read in file ==! 433 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F(i,j) read in eddy_diffusivity_2D.nc file' 569 IF(lwp) WRITE(numout,*) ' surface eddy diffusivity = constant = ', aht0, ' m2/s' 570 aeiu(:,:,1) = aei0 ! constant surface value 571 aeiv(:,:,1) = aei0 572 CALL ldf_c1d( 'TRA', aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) 573 ! 574 CASE ( -20 ) !-- fixed horizontal shape read in file --! 575 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F(i,j) read in eddy_diffusivity_2D.nc file' 434 576 CALL iom_open ( 'eddy_induced_velocity_2D.nc', inum ) 435 577 CALL iom_get ( inum, jpdom_data, 'aeiu', aeiu(:,:,1) ) 436 578 CALL iom_get ( inum, jpdom_data, 'aeiv', aeiv(:,:,1) ) 437 579 CALL iom_close( inum ) 438 DO jk = 2, jpk 580 DO jk = 2, jpkm1 439 581 aeiu(:,:,jk) = aeiu(:,:,1) 440 582 aeiv(:,:,jk) = aeiv(:,:,1) 441 583 END DO 442 584 ! 443 CASE( 20 ) !== fixed horizontal shape ==! 444 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or bilap case)' 445 CALL ldf_c2d( 'TRA', 'LAP', rn_aeiv_0, aeiu, aeiv ) ! surface value proportional to scale factor 446 ! 447 CASE( 21 ) !== time varying 2D field ==! 448 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( latitude, longitude, time )' 449 IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )' 585 CASE( 20 ) !-- fixed horizontal shape --! 586 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( e1, e2 )' 587 IF(lwp) WRITE(numout,*) ' using a fixed diffusive velocity = ', rn_Ue, ' m/s and Le = Max(e1,e2)' 588 IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, ' m2/s for e1=1°)' 589 CALL ldf_c2d( 'TRA', zUfac , inn , aeiu, aeiv ) ! value proportional to scale factor^inn 590 ! 591 CASE( 21 ) !-- time varying 2D field --! 592 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( latitude, longitude, time )' 593 IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )' 594 IF(lwp) WRITE(numout,*) ' maximum allowed value: aei0 = ', aei0, ' m2/s' 450 595 ! 451 596 l_ldfeiv_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90 452 597 ! 453 CASE( -30 ) !== fixed 3D shape read in file ==!454 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixingcoef. = F(i,j,k) read in eddy_diffusivity_3D.nc file'598 CASE( -30 ) !-- fixed 3D shape read in file --! 599 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file' 455 600 CALL iom_open ( 'eddy_induced_velocity_3D.nc', inum ) 456 601 CALL iom_get ( inum, jpdom_data, 'aeiu', aeiu ) … … 458 603 CALL iom_close( inum ) 459 604 ! 460 CASE( 30 ) !== fixed 3D shape ==! 461 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( latitude, longitude, depth )' 462 CALL ldf_c2d( 'TRA', 'LAP', rn_aeiv_0, aeiu, aeiv ) ! surface value proportional to scale factor 463 ! ! reduction with depth 464 CALL ldf_c1d( 'TRA', r1_4, aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) 605 CASE( 30 ) !-- fixed 3D shape --! 606 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( latitude, longitude, depth )' 607 CALL ldf_c2d( 'TRA', zUfac , inn , aeiu, aeiv ) ! surface value proportional to scale factor^inn 608 CALL ldf_c1d( 'TRA', aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) ! reduction with depth 465 609 ! 466 610 CASE DEFAULT … … 468 612 END SELECT 469 613 ! 470 ELSE 471 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity param is NOT used neither diagnosed' 472 ln_ldfeiv_dia = .FALSE. 614 IF( .NOT.l_ldfeiv_time ) THEN !* mask if No time variation 615 DO jk = 1, jpkm1 616 aeiu(:,:,jk) = aeiu(:,:,jk) * umask(:,:,jk) 617 ahtv(:,:,jk) = ahtv(:,:,jk) * vmask(:,:,jk) 618 END DO 619 ENDIF 620 ! 473 621 ENDIF 474 622 ! … … 493 641 INTEGER :: ji, jj, jk ! dummy loop indices 494 642 REAL(wp) :: zfw, ze3w, zn2, z1_f20, zaht, zaht_min, zzaei ! local scalars 495 REAL(wp), DIMENSION(jpi,jpj) :: zn, zah, zhw, z ross, zaeiw ! 2D workspace496 !!---------------------------------------------------------------------- 497 ! 498 zn (:,:) = 0._wp! Local initialization499 zhw 500 zah 501 z ross(:,:) = 0._wp643 REAL(wp), DIMENSION(jpi,jpj) :: zn, zah, zhw, zRo, zaeiw ! 2D workspace 644 !!---------------------------------------------------------------------- 645 ! 646 zn (:,:) = 0._wp ! Local initialization 647 zhw(:,:) = 5._wp 648 zah(:,:) = 0._wp 649 zRo(:,:) = 0._wp 502 650 ! ! Compute lateral diffusive coefficient at T-point 503 651 IF( ln_traldf_triad ) THEN … … 538 686 END DO 539 687 END DO 540 END 688 ENDIF 541 689 542 690 DO jj = 2, jpjm1 543 691 DO ji = fs_2, fs_jpim1 ! vector opt. 544 692 zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 ) 545 ! Rossby radius at w-point taken < 40km and > 2km546 z ross(ji,jj) = MAX( MIN( .4 * zn(ji,jj) / zfw, 40.e3 ), 2.e3)693 ! Rossby radius at w-point taken betwenn 2 km and 40km 694 zRo(ji,jj) = MAX( 2.e3 , MIN( .4 * zn(ji,jj) / zfw, 40.e3 ) ) 547 695 ! Compute aeiw by multiplying Ro^2 and T^-1 548 zaeiw(ji,jj) = z ross(ji,jj) * zross(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)696 zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1) 549 697 END DO 550 698 END DO … … 554 702 DO jj = 2, jpjm1 555 703 DO ji = fs_2, fs_jpim1 ! vector opt. 556 zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) 704 zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease 557 705 zaeiw(ji,jj) = MIN( zzaei , paei0 ) ! Max value = paei0 558 706 END DO … … 620 768 DO jj = 1, jpjm1 621 769 DO ji = 1, fs_jpim1 ! vector opt. 622 zpsi_uw(ji,jj,jk) = - 0.25_wp* e2u(ji,jj) * ( wslpi(ji,jj,jk ) + wslpi(ji+1,jj,jk) ) &623 & 624 zpsi_vw(ji,jj,jk) = - 0.25_wp* e1v(ji,jj) * ( wslpj(ji,jj,jk ) + wslpj(ji,jj+1,jk) ) &625 & 770 zpsi_uw(ji,jj,jk) = - r1_4 * e2u(ji,jj) * ( wslpi(ji,jj,jk ) + wslpi(ji+1,jj,jk) ) & 771 & * ( aeiu (ji,jj,jk-1) + aeiu (ji ,jj,jk) ) * umask(ji,jj,jk) 772 zpsi_vw(ji,jj,jk) = - r1_4 * e1v(ji,jj) * ( wslpj(ji,jj,jk ) + wslpj(ji,jj+1,jk) ) & 773 & * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj ,jk) ) * vmask(ji,jj,jk) 626 774 END DO 627 775 END DO -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_oper.F90
r9023 r9490 10 10 !! obs_surf_opt : Compute the model counterpart of surface data 11 11 !!---------------------------------------------------------------------- 12 13 !! * Modules used 14 USE par_kind, ONLY : & ! Precision variables 15 & wp 16 USE in_out_manager ! I/O manager 17 USE obs_inter_sup ! Interpolation support 18 USE obs_inter_h2d, ONLY : & ! Horizontal interpolation to the obs pt 19 & obs_int_h2d, & 20 & obs_int_h2d_init 21 USE obs_averg_h2d, ONLY : & ! Horizontal averaging to the obs footprint 22 & obs_avg_h2d, & 23 & obs_avg_h2d_init, & 24 & obs_max_fpsize 25 USE obs_inter_z1d, ONLY : & ! Vertical interpolation to the obs pt 26 & obs_int_z1d, & 27 & obs_int_z1d_spl 28 USE obs_const, ONLY : & ! Obs fill value 29 & obfillflt 30 USE dom_oce, ONLY : & 31 & glamt, glamf, & 32 & gphit, gphif 33 USE lib_mpp, ONLY : & ! Warning and stopping routines 34 & ctl_warn, ctl_stop 35 USE sbcdcy, ONLY : & ! For calculation of where it is night-time 36 & sbc_dcy, nday_qsr 37 USE obs_grid, ONLY : & 38 & obs_level_search 12 USE obs_inter_sup ! Interpolation support 13 USE obs_inter_h2d, ONLY : obs_int_h2d, obs_int_h2d_init ! Horizontal interpolation to the obs pt 14 USE obs_averg_h2d, ONLY : obs_avg_h2d, obs_avg_h2d_init, obs_max_fpsize ! Horizontal averaging to the obs footprint 15 USE obs_inter_z1d, ONLY : obs_int_z1d, obs_int_z1d_spl ! Vertical interpolation to the obs pt 16 USE obs_const , ONLY : obfillflt ! Obs fill value 17 USE dom_oce, ONLY : glamt, glamf, gphit, gphif ! lat/lon of ocean grid-points 18 USE lib_mpp, ONLY : ctl_warn, ctl_stop ! Warning and stopping routines 19 USE sbcdcy, ONLY : sbc_dcy, nday_qsr ! For calculation of where it is night-time 20 USE obs_grid, ONLY : obs_level_search 21 ! 22 USE par_kind , ONLY : wp ! Precision variables 23 USE in_out_manager ! I/O manager 39 24 40 25 IMPLICIT NONE 41 42 !! * Routine accessibility43 26 PRIVATE 44 27 45 PUBLIC obs_prof_opt, & ! Compute the model counterpart of profile obs 46 & obs_surf_opt ! Compute the model counterpart of surface obs 47 48 INTEGER, PARAMETER, PUBLIC :: & 49 & imaxavtypes = 20 ! Max number of daily avgd obs types 28 PUBLIC obs_prof_opt !: Compute the model counterpart of profile obs 29 PUBLIC obs_surf_opt !: Compute the model counterpart of surface obs 30 31 INTEGER, PARAMETER, PUBLIC :: imaxavtypes = 20 !: Max number of daily avgd obs types 50 32 51 33 !!---------------------------------------------------------------------- 52 !! NEMO/OPA 3.3 , NEMO Consortium (2010)34 !! NEMO/OPA 4.0 , NEMO Consortium (2018) 53 35 !! $Id$ 54 36 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 55 37 !!---------------------------------------------------------------------- 56 57 38 CONTAINS 58 59 39 60 40 SUBROUTINE obs_prof_opt( prodatqc, kt, kpi, kpj, kpk, & … … 64 44 & plam1, plam2, pphi1, pphi2, & 65 45 & k1dint, k2dint, kdailyavtypes ) 66 67 46 !!----------------------------------------------------------------------- 68 !!69 47 !! *** ROUTINE obs_pro_opt *** 70 48 !! … … 114 92 !! ! 17-02 (M. Martin) Include generalised vertical coordinate changes 115 93 !!----------------------------------------------------------------------- 116 117 !! * Modules used118 94 USE obs_profiles_def ! Definition of storage space for profile obs. 119 95 120 96 IMPLICIT NONE 121 97 122 !! * Arguments 123 TYPE(obs_prof), INTENT(INOUT) :: & 124 & prodatqc ! Subset of profile data passing QC 125 INTEGER, INTENT(IN) :: kt ! Time step 126 INTEGER, INTENT(IN) :: kpi ! Model grid parameters 127 INTEGER, INTENT(IN) :: kpj 128 INTEGER, INTENT(IN) :: kpk 129 INTEGER, INTENT(IN) :: kit000 ! Number of the first time step 130 ! (kit000-1 = restart time) 131 INTEGER, INTENT(IN) :: k1dint ! Vertical interpolation type (see header) 132 INTEGER, INTENT(IN) :: k2dint ! Horizontal interpolation type (see header) 133 INTEGER, INTENT(IN) :: kdaystp ! Number of time steps per day 134 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 135 & pvar1, & ! Model field 1 136 & pvar2, & ! Model field 2 137 & pmask1, & ! Land-sea mask 1 138 & pmask2 ! Land-sea mask 2 139 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj) :: & 140 & plam1, & ! Model longitudes for variable 1 141 & plam2, & ! Model longitudes for variable 2 142 & pphi1, & ! Model latitudes for variable 1 143 & pphi2 ! Model latitudes for variable 2 144 REAL(KIND=wp), INTENT(IN), DIMENSION(kpi,kpj,kpk) :: & 145 & pgdept, & ! Model array of depth T levels 146 & pgdepw ! Model array of depth W levels 147 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: & 148 & kdailyavtypes ! Types for daily averages 98 TYPE(obs_prof), INTENT(inout) :: prodatqc ! Subset of profile data passing QC 99 INTEGER , INTENT(in ) :: kt ! Time step 100 INTEGER , INTENT(in ) :: kpi, kpj, kpk ! Model grid parameters 101 INTEGER , INTENT(in ) :: kit000 ! Number of the first time step (kit000-1 = restart time) 102 INTEGER , INTENT(in ) :: k1dint ! Vertical interpolation type (see header) 103 INTEGER , INTENT(in ) :: k2dint ! Horizontal interpolation type (see header) 104 INTEGER , INTENT(in ) :: kdaystp ! Number of time steps per day 105 REAL(KIND=wp) , INTENT(in ), DIMENSION(kpi,kpj,kpk) :: pvar1 , pvar2 ! Model field 1 and 2 106 REAL(KIND=wp) , INTENT(in ), DIMENSION(kpi,kpj,kpk) :: pmask1, pmask2 ! Land-sea mask 1 and 2 107 REAL(KIND=wp) , INTENT(in ), DIMENSION(kpi,kpj) :: plam1 , plam2 ! Model longitude 1 and 2 108 REAL(KIND=wp) , INTENT(in ), DIMENSION(kpi,kpj) :: pphi1 , pphi2 ! Model latitudes 1 and 2 109 REAL(KIND=wp) , INTENT(in ), DIMENSION(kpi,kpj,kpk) :: pgdept, pgdepw ! depth of T and W levels 110 INTEGER, DIMENSION(imaxavtypes), OPTIONAL :: kdailyavtypes ! Types for daily averages 149 111 150 112 !! * Local declarations … … 706 668 !! ! 17-03 (M. Martin) Added horizontal averaging options 707 669 !!----------------------------------------------------------------------- 708 709 !! * Modules used710 670 USE obs_surf_def ! Definition of storage space for surface observations 711 671 712 672 IMPLICIT NONE 713 673 714 !! * Arguments715 674 TYPE(obs_surf), INTENT(INOUT) :: & 716 675 & surfdataqc ! Subset of surface data passing QC … … 866 825 DO ji = 0, imaxifp 867 826 imodi = surfdataqc%mi(jobs) - int(imaxifp/2) + ji - 1 868 827 ! 869 828 !Deal with wrap around in longitude 870 829 IF ( imodi < 1 ) imodi = imodi + jpiglo 871 830 IF ( imodi > jpiglo ) imodi = imodi - jpiglo 872 831 ! 873 832 DO jj = 0, imaxjfp 874 833 imodj = surfdataqc%mj(jobs) - int(imaxjfp/2) + jj - 1 … … 877 836 IF ( imodj < 1 ) imodj = 1 878 837 IF ( imodj > jpjglo ) imodj = jpjglo 879 838 ! 880 839 igrdip1(ji+1,jj+1,iobs) = imodi 881 840 igrdjp1(ji+1,jj+1,iobs) = imodj 882 841 ! 883 842 IF ( ji >= 1 .AND. jj >= 1 ) THEN 884 843 igrdi(ji,jj,iobs) = imodi 885 844 igrdj(ji,jj,iobs) = imodj 886 845 ENDIF 887 846 ! 888 847 END DO 889 848 END DO … … 1010 969 & ) 1011 970 ENDIF 1012 971 ! 1013 972 surfdataqc%nsurfup = surfdataqc%nsurfup + isurf 1014 973 ! 1015 974 END SUBROUTINE obs_surf_opt 1016 975 976 !!====================================================================== 1017 977 END MODULE obs_oper -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/OBS/obs_prep.F90
r9023 r9490 108 108 imin0 = ( nn_time0 - ihou0 * 100 ) 109 109 110 icycle = n o ! Assimilation cycle110 icycle = nn_no ! Assimilation cycle 111 111 112 112 ! Diagnotics counters for various failures. … … 339 339 imin0 = ( nn_time0 - ihou0 * 100 ) 340 340 341 icycle = n o ! Assimilation cycle341 icycle = nn_no ! Assimilation cycle 342 342 343 343 ! Diagnotics counters for various failures. 344 344 345 iotdobs = 0346 igrdobs = 0345 iotdobs = 0 346 igrdobs = 0 347 347 iosdv1obs = 0 348 348 iosdv2obs = 0 … … 884 884 !! ! 2007-01 (K. Mogensen) Original 885 885 !!---------------------------------------------------------------------- 886 !! * Arguments887 886 INTEGER, INTENT(IN) :: kobsno ! Number of observations 888 887 INTEGER, DIMENSION(kobsno), INTENT(IN ) :: & … … 924 923 !! ** Action : 925 924 !! 926 !! History : 927 !! ! 2007-03 (A. Weaver, K. Mogensen) Original 928 !! ! 2007-06 (K. Mogensen et al) Reject obs. near land. 929 !!---------------------------------------------------------------------- 930 !! * Modules used 931 932 !! * Arguments 933 INTEGER, INTENT(IN) :: kobsno ! Total number of observations 934 INTEGER, INTENT(IN) :: kpi ! Number of grid points in (i,j) 935 INTEGER, INTENT(IN) :: kpj 936 INTEGER, DIMENSION(kobsno), INTENT(IN) :: & 937 & kobsi, & ! Observation (i,j) coordinates 938 & kobsj 939 REAL(KIND=wp), DIMENSION(kobsno), INTENT(IN) :: & 940 & pobslam, & ! Observation (lon,lat) coordinates 941 & pobsphi 942 REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: & 943 & plam, pphi ! Model (lon,lat) coordinates 944 REAL(KIND=wp), DIMENSION(kpi,kpj), INTENT(IN) :: & 945 & pmask ! Land mask array 946 INTEGER, DIMENSION(kobsno), INTENT(INOUT) :: & 947 & kobsqc ! Observation quality control 948 INTEGER, INTENT(INOUT) :: kosdobs ! Observations outside space domain 949 INTEGER, INTENT(INOUT) :: klanobs ! Observations within a model land cell 950 INTEGER, INTENT(INOUT) :: knlaobs ! Observations near land 951 INTEGER, INTENT(INOUT) :: kbdyobs ! Observations near boundary 952 LOGICAL, INTENT(IN) :: ld_nea ! Flag observations near land 953 LOGICAL, INTENT(IN) :: ld_bound_reject ! Flag observations near open boundary 954 INTEGER, INTENT(IN) :: kqc_cutoff ! Cutoff QC value 955 956 !! * Local declarations 957 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 958 & zgmsk ! Grid mask 959 960 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 961 & zbmsk ! Boundary mask 962 REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 963 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: & 964 & zglam, & ! Model longitude at grid points 965 & zgphi ! Model latitude at grid points 966 INTEGER, DIMENSION(2,2,kobsno) :: & 967 & igrdi, & ! Grid i,j 968 & igrdj 969 LOGICAL :: lgridobs ! Is observation on a model grid point. 970 INTEGER :: iig, ijg ! i,j of observation on model grid point. 971 INTEGER :: jobs, ji, jj 925 !! History : 2007-03 (A. Weaver, K. Mogensen) Original 926 !! ! 2007-06 (K. Mogensen et al) Reject obs. near land. 927 !!---------------------------------------------------------------------- 928 INTEGER , INTENT(in ) :: kobsno ! Total number of observations 929 INTEGER , INTENT(in ) :: kpi , kpj ! Number of grid points in (i,j) 930 INTEGER , INTENT(in ), DIMENSION(kobsno) :: kobsi , kobsj ! Observation (i,j) coordinates 931 REAL(wp), INTENT(in ), DIMENSION(kobsno) :: pobslam, pobsphi ! Observation (lon,lat) coordinates 932 REAL(wp), INTENT(in ), DIMENSION(kpi,kpj) :: plam , pphi ! Model (lon,lat) coordinates 933 REAL(wp), INTENT(in ), DIMENSION(kpi,kpj) :: pmask ! Land mask array 934 INTEGER , INTENT(inout), DIMENSION(kobsno) :: kobsqc ! Observation quality control 935 INTEGER , INTENT(inout) :: kosdobs ! Observations outside space domain 936 INTEGER , INTENT(inout) :: klanobs ! Observations within a model land cell 937 INTEGER , INTENT(inout) :: knlaobs ! Observations near land 938 INTEGER , INTENT(inout) :: kbdyobs ! Observations near boundary 939 LOGICAL , INTENT(in ) :: ld_nea ! Flag observations near land 940 LOGICAL , INTENT(in ) :: ld_bound_reject ! Flag observations near open boundary 941 INTEGER , INTENT(in ) :: kqc_cutoff ! Cutoff QC value 942 ! 943 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: zgmsk ! Grid mask 944 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: zbmsk ! Boundary mask 945 REAL(KIND=wp), DIMENSION(jpi,jpj) :: zbdymask 946 REAL(KIND=wp), DIMENSION(2,2,kobsno) :: zglam, zgphi ! Model Lon/lat at grid points 947 INTEGER , DIMENSION(2,2,kobsno) :: igrdi, igrdj ! Grid i,j 948 LOGICAL :: lgridobs ! Is observation on a model grid point. 949 INTEGER :: iig, ijg ! i,j of observation on model grid point. 950 INTEGER :: jobs, ji, jj 951 !!---------------------------------------------------------------------- 972 952 973 953 ! Get grid point indices … … 1100 1080 ENDIF 1101 1081 ENDIF 1102 1082 ! 1103 1083 END DO 1104 1084 ! 1105 1085 END SUBROUTINE obs_coo_spc_2d 1086 1106 1087 1107 1088 SUBROUTINE obs_coo_spc_3d( kprofno, kobsno, kpstart, kpend, & … … 1198 1179 INTEGER :: iig, ijg ! i,j of observation on model grid point. 1199 1180 INTEGER :: jobs, jobsp, jk, ji, jj 1181 !!---------------------------------------------------------------------- 1200 1182 1201 1183 ! Get grid point indices … … 1359 1341 ENDIF 1360 1342 ENDIF 1361 1343 ! 1362 1344 END DO 1363 1345 END DO 1364 1346 ! 1365 1347 END SUBROUTINE obs_coo_spc_3d 1348 1366 1349 1367 1350 SUBROUTINE obs_pro_rej( profdata, kqc_cutoff ) … … 1377 1360 !! References : 1378 1361 !! 1379 !! History : 1380 !! ! 2007-10 (K. Mogensen) Original code 1381 !!---------------------------------------------------------------------- 1382 !! * Modules used 1383 !! * Arguments 1384 TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data 1385 INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value 1386 1387 !! * Local declarations 1362 !! History : 2007-10 (K. Mogensen) Original code 1363 !!---------------------------------------------------------------------- 1364 TYPE(obs_prof), INTENT(inout) :: profdata ! Profile data 1365 INTEGER , INTENT(in ) :: kqc_cutoff ! QC cutoff value 1366 ! 1388 1367 INTEGER :: jprof 1389 1368 INTEGER :: jvar 1390 1369 INTEGER :: jobs 1370 !!---------------------------------------------------------------------- 1391 1371 1392 1372 ! Loop over profiles … … 1411 1391 1412 1392 END DO 1413 1393 ! 1414 1394 END SUBROUTINE obs_pro_rej 1395 1415 1396 1416 1397 SUBROUTINE obs_uv_rej( profdata, knumu, knumv, kqc_cutoff ) … … 1426 1407 !! References : 1427 1408 !! 1428 !! History : 1429 !! ! 2009-2 (K. Mogensen) Original code 1430 !!---------------------------------------------------------------------- 1431 !! * Modules used 1432 !! * Arguments 1409 !! History : 2009-2 (K. Mogensen) Original code 1410 !!---------------------------------------------------------------------- 1433 1411 TYPE(obs_prof), INTENT(INOUT) :: profdata ! Profile data 1434 1412 INTEGER, INTENT(INOUT) :: knumu ! Number of u rejected 1435 1413 INTEGER, INTENT(INOUT) :: knumv ! Number of v rejected 1436 1414 INTEGER, INTENT(IN) :: kqc_cutoff ! QC cutoff value 1437 1438 !! * Local declarations 1415 ! 1439 1416 INTEGER :: jprof 1440 1417 INTEGER :: jvar 1441 1418 INTEGER :: jobs 1442 1443 ! Loop over profiles 1444 1445 DO jprof = 1, profdata%nprof 1446 1419 !!---------------------------------------------------------------------- 1420 1421 DO jprof = 1, profdata%nprof !== Loop over profiles ==! 1422 ! 1447 1423 IF ( ( profdata%npvsta(jprof,1) /= profdata%npvsta(jprof,2) ) .OR. & 1448 1424 & ( profdata%npvend(jprof,1) /= profdata%npvend(jprof,2) ) ) THEN 1449 1425 ! 1450 1426 CALL ctl_stop('U,V profiles inconsistent in obs_uv_rej') 1451 1427 RETURN 1452 1453 ENDIF 1454 1428 ! 1429 ENDIF 1430 ! 1455 1431 DO jobs = profdata%npvsta(jprof,1), profdata%npvend(jprof,1) 1456 1432 ! 1457 1433 IF ( ( profdata%var(1)%nvqc(jobs) > kqc_cutoff ) .AND. & 1458 1434 & ( profdata%var(2)%nvqc(jobs) <= kqc_cutoff) ) THEN … … 1465 1441 knumu = knumu + 1 1466 1442 ENDIF 1467 1443 ! 1468 1444 END DO 1469 1445 ! 1470 1446 END DO 1471 1447 ! 1472 1448 END SUBROUTINE obs_uv_rej 1473 1449 1450 !!===================================================================== 1474 1451 END MODULE obs_prep -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/tradmp.F90
r9168 r9490 197 197 ENDIF 198 198 ! 199 IF( ln_tradmp ) THEN199 IF( ln_tradmp ) THEN 200 200 ! ! Allocate arrays 201 201 IF( tra_dmp_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'tra_dmp_init: unable to allocate arrays' ) … … 212 212 ! so can damp to something other than intitial conditions files? 213 213 !!gm: In principle yes. Nevertheless, we can't anticipate demands that have never been formulated. 214 IF( .NOT.ln_tsd_ tradmp ) THEN214 IF( .NOT.ln_tsd_dmp ) THEN 215 215 IF(lwp) WRITE(numout,*) 216 IF(lwp) WRITE(numout, *) ' read T-S data not initialized, we force ln_tsd_ tradmp=T'216 IF(lwp) WRITE(numout, *) ' read T-S data not initialized, we force ln_tsd_dmp=T' 217 217 CALL dta_tsd_init( ld_tradmp=ln_tradmp ) ! forces the initialisation of T-S data 218 218 ENDIF -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90
r9190 r9490 37 37 PUBLIC tra_ldf ! called by step.F90 38 38 PUBLIC tra_ldf_init ! called by nemogcm.F90 39 !40 INTEGER :: nldf = 0 ! type of lateral diffusion used defined from ln_traldf_... (namlist logicals)41 39 42 40 !! * Substitutions 43 41 # include "vectopt_loop_substitute.h90" 44 42 !!---------------------------------------------------------------------- 45 !! NEMO/OPA 3.7 , NEMO Consortium (2015)43 !! NEMO/OPA 4.0 , NEMO Consortium (2018) 46 44 !! $Id$ 47 45 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 68 66 ENDIF 69 67 ! 70 SELECT CASE ( nldf )!* compute lateral mixing trend and add it to the general trend68 SELECT CASE ( nldf_tra ) !* compute lateral mixing trend and add it to the general trend 71 69 CASE ( np_lap ) ! laplacian: iso-level operator 72 70 CALL tra_ldf_lap ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsa, jpts, 1 ) … … 76 74 CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts, 1 ) 77 75 CASE ( np_blp , np_blp_i , np_blp_it ) ! bilaplacian: iso-level & iso-neutral operators 78 CALL tra_ldf_blp ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb , tsa, jpts, nldf )76 CALL tra_ldf_blp ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb , tsa, jpts, nldf_tra ) 79 77 END SELECT 80 78 ! … … 101 99 !! ** Purpose : Choice of the operator for the lateral tracer diffusion 102 100 !! 103 !! ** Method : set nldf from the namtra_ldf logicals101 !! ** Method : set nldf_tra from the namtra_ldf logicals 104 102 !!---------------------------------------------------------------------- 105 103 INTEGER :: ioptio, ierr ! temporary integers … … 112 110 WRITE(numout,*) ' Namelist namtra_ldf: already read in ldftra module' 113 111 WRITE(numout,*) ' see ldf_tra_init report for lateral mixing parameters' 114 ENDIF 115 ! !== use of lateral operator or not ==! 116 nldf = np_ERROR 117 ioptio = 0 118 IF( ln_traldf_NONE ) THEN ; nldf = np_no_ldf ; ioptio = ioptio + 1 ; ENDIF 119 IF( ln_traldf_lap ) THEN ; ioptio = ioptio + 1 ; ENDIF 120 IF( ln_traldf_blp ) THEN ; ioptio = ioptio + 1 ; ENDIF 121 IF( ioptio /= 1 ) CALL ctl_stop( 'tra_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' ) 122 ! 123 IF( .NOT.ln_traldf_NONE ) THEN !== direction ==>> type of operator ==! 124 ioptio = 0 125 IF( ln_traldf_lev ) ioptio = ioptio + 1 126 IF( ln_traldf_hor ) ioptio = ioptio + 1 127 IF( ln_traldf_iso ) ioptio = ioptio + 1 128 IF( ioptio /= 1 ) CALL ctl_stop( 'tra_ldf_init: use ONE direction (level/hor/iso)' ) 112 WRITE(numout,*) 129 113 ! 130 ! ! defined the type of lateral diffusion from ln_traldf_... logicals 131 ierr = 0 132 IF( ln_traldf_lap ) THEN ! laplacian operator 133 IF ( ln_zco ) THEN ! z-coordinate 134 IF ( ln_traldf_lev ) nldf = np_lap ! iso-level = horizontal (no rotation) 135 IF ( ln_traldf_hor ) nldf = np_lap ! iso-level = horizontal (no rotation) 136 IF ( ln_traldf_iso ) nldf = np_lap_i ! iso-neutral: standard ( rotation) 137 IF ( ln_traldf_triad ) nldf = np_lap_it ! iso-neutral: triad ( rotation) 138 ENDIF 139 IF ( ln_zps ) THEN ! z-coordinate with partial step 140 IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed 141 IF ( ln_traldf_hor ) nldf = np_lap ! horizontal (no rotation) 142 IF ( ln_traldf_iso ) nldf = np_lap_i ! iso-neutral: standard (rotation) 143 IF ( ln_traldf_triad ) nldf = np_lap_it ! iso-neutral: triad (rotation) 144 ENDIF 145 IF ( ln_sco ) THEN ! s-coordinate 146 IF ( ln_traldf_lev ) nldf = np_lap ! iso-level (no rotation) 147 IF ( ln_traldf_hor ) nldf = np_lap_i ! horizontal ( rotation) 148 IF ( ln_traldf_iso ) nldf = np_lap_i ! iso-neutral: standard ( rotation) 149 IF ( ln_traldf_triad ) nldf = np_lap_it ! iso-neutral: triad ( rotation) 150 ENDIF 151 ENDIF 152 ! 153 IF( ln_traldf_blp ) THEN ! bilaplacian operator 154 IF ( ln_zco ) THEN ! z-coordinate 155 IF ( ln_traldf_lev ) nldf = np_blp ! iso-level = horizontal (no rotation) 156 IF ( ln_traldf_hor ) nldf = np_blp ! iso-level = horizontal (no rotation) 157 IF ( ln_traldf_iso ) nldf = np_blp_i ! iso-neutral: standard ( rotation) 158 IF ( ln_traldf_triad ) nldf = np_blp_it ! iso-neutral: triad ( rotation) 159 ENDIF 160 IF ( ln_zps ) THEN ! z-coordinate with partial step 161 IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed 162 IF ( ln_traldf_hor ) nldf = np_blp ! horizontal (no rotation) 163 IF ( ln_traldf_iso ) nldf = np_blp_i ! iso-neutral: standard ( rotation) 164 IF ( ln_traldf_triad ) nldf = np_blp_it ! iso-neutral: triad ( rotation) 165 ENDIF 166 IF ( ln_sco ) THEN ! s-coordinate 167 IF ( ln_traldf_lev ) nldf = np_blp ! iso-level (no rotation) 168 IF ( ln_traldf_hor ) nldf = np_blp_it ! horizontal ( rotation) 169 IF ( ln_traldf_iso ) nldf = np_blp_i ! iso-neutral: standard ( rotation) 170 IF ( ln_traldf_triad ) nldf = np_blp_it ! iso-neutral: triad ( rotation) 171 ENDIF 172 ENDIF 173 IF( ierr == 1 ) CALL ctl_stop( 'iso-level in z-partial step, not allowed' ) 174 ENDIF 175 ! 176 IF( ln_ldfeiv .AND. .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) ) & 177 & CALL ctl_stop( 'eddy induced velocity on tracers requires iso-neutral laplacian diffusion' ) 178 IF( ln_isfcav .AND. ln_traldf_triad ) & 179 & CALL ctl_stop( ' ice shelf cavity and traldf_triad not tested' ) 180 ! 181 IF( nldf == np_lap_i .OR. nldf == np_lap_it .OR. & 182 & nldf == np_blp_i .OR. nldf == np_blp_it ) l_ldfslp = .TRUE. ! slope of neutral surfaces required 183 ! 184 IF(lwp) THEN 185 WRITE(numout,*) 186 SELECT CASE( nldf ) 114 SELECT CASE( nldf_tra ) ! print the choice of operator 187 115 CASE( np_no_ldf ) ; WRITE(numout,*) ' ==>>> NO lateral diffusion' 188 116 CASE( np_lap ) ; WRITE(numout,*) ' ==>>> laplacian iso-level operator' -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_blp.F90
r9124 r9490 33 33 PUBLIC tra_ldf_blp ! called by traldf.F90 34 34 35 ! ! Flag to control the type of lateral diffusive operator36 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 ! error in specification of lateral diffusion37 INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 ! without operator (i.e. no lateral diffusive trend)38 ! !! laplacian ! bilaplacian !39 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 ! iso-level operator40 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 , np_blp_i = 21 ! standard iso-neutral or geopotential operator41 INTEGER, PARAMETER, PUBLIC :: np_lap_it = 12 , np_blp_it = 22 ! triad iso-neutral or geopotential operator42 43 35 LOGICAL :: l_ptr ! flag to compute poleward transport 44 36 LOGICAL :: l_hst ! flag to compute heat transport -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r9436 r9490 234 234 INTEGER :: ji ! dummy loop indices 235 235 INTEGER :: ios, ilocal_comm ! local integers 236 CHARACTER(len=120), DIMENSION( 30) :: cltxt, cltxt2, clnam236 CHARACTER(len=120), DIMENSION(60) :: cltxt, cltxt2, clnam 237 237 !! 238 238 NAMELIST/namctl/ ln_ctl , nn_print, nn_ictls, nn_ictle, & -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/step.F90
r9485 r9490 122 122 123 123 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 124 ! Ocean physics update (ua, va, tsa used as workspace)124 ! Ocean physics update 125 125 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 126 126 ! THERMODYNAMICS … … 208 208 209 209 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 210 ! diagnostics and outputs (ua, va, tsa used as workspace)211 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 212 IF( lk_floats ) CALL flo_stp ( kstp )! drifting Floats213 IF( ln_diacfl ) CALL dia_cfl ( kstp )! Courant number diagnostics214 IF( lk_diahth ) CALL dia_hth ( kstp )! Thermocline depth (20 degres isotherm depth)215 IF( lk_diadct ) CALL dia_dct ( kstp )! Transports216 CALL dia_ar5 ( kstp )! ar5 diag210 ! diagnostics and outputs 211 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 212 IF( lk_floats ) CALL flo_stp ( kstp ) ! drifting Floats 213 IF( ln_diacfl ) CALL dia_cfl ( kstp ) ! Courant number diagnostics 214 IF( lk_diahth ) CALL dia_hth ( kstp ) ! Thermocline depth (20 degres isotherm depth) 215 IF( lk_diadct ) CALL dia_dct ( kstp ) ! Transports 216 CALL dia_ar5 ( kstp ) ! ar5 diag 217 217 IF( lk_diaharm ) CALL dia_harm( kstp ) ! Tidal harmonic analysis 218 CALL dia_wri ( kstp )! ocean model: outputs218 CALL dia_wri ( kstp ) ! ocean model: outputs 219 219 ! 220 220 IF( ln_crs ) CALL crs_fld ( kstp ) ! ocean model: online field coarsening & output
Note: See TracChangeset
for help on using the changeset viewer.