Changeset 6598
- Timestamp:
- 2016-05-23T14:22:42+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_surft_in_asminc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r5540 r6598 1 2 1 3 MODULE asminc 2 4 !!====================================================================== … … 14 16 15 17 !!---------------------------------------------------------------------- 16 !! 'key_asminc' : Switch on the assimilation increment interface 17 !!---------------------------------------------------------------------- 18 !! asm_inc_init : Initialize the increment arrays and IAU weights 19 !! calc_date : Compute the calendar date YYYYMMDD on a given step 20 !! tra_asm_inc : Apply the tracer (T and S) increments 21 !! dyn_asm_inc : Apply the dynamic (u and v) increments 22 !! ssh_asm_inc : Apply the SSH increment 23 !! seaice_asm_inc : Apply the seaice increment 18 !! 'key_asminc' : Switch on the assimilation increment interface 19 !!---------------------------------------------------------------------- 20 !! asm_inc_init : Initialize the increment arrays and IAU weights 21 !! calc_date : Compute the calendar date YYYYMMDD on a given step 22 !! tra_asm_inc : Apply the tracer (T and S) increments 23 !! dyn_asm_inc : Apply the dynamic (u and v) increments 24 !! ssh_asm_inc : Apply the SSH increment 25 !! seaice_asm_inc : Apply the seaice increment 26 !! logchl_asm_inc : Apply the logchl increment 27 !! pco2_asm_inc : Apply the pco2 or fco2 increment 24 28 !!---------------------------------------------------------------------- 25 29 USE wrk_nemo ! Memory Allocation … … 32 36 USE zpshde ! Partial step : Horizontal Derivative 33 37 USE iom ! Library to read input files 34 USE asmpar ! Parameters for the ass milation interface38 USE asmpar ! Parameters for the assimilation interface 35 39 USE c1d ! 1D initialization 36 40 USE in_out_manager ! I/O manager 37 41 USE lib_mpp ! MPP library 38 #if defined key_lim2 39 USE ice_2 ! LIM240 #endif 42 USE bioanal ! LogChl balancing 43 USE pco2analysis ! pCO2 balancing 44 USE phycst, ONLY : rt0 ! Freezing point of water in Kelvin 41 45 USE sbc_oce ! Surface boundary condition variables. 46 47 USE eosbn2, only: tfreez 48 49 USE zdfmxl, ONLY : & 50 & hmld_tref, & 51 & hmld, & 52 & hmlp, & 53 & hmlpt 54 USE bdy_oce, ONLY: bdytmask 55 USE histcom 42 56 43 57 IMPLICIT NONE … … 50 64 PUBLIC ssh_asm_inc !: Apply the SSH increment 51 65 PUBLIC seaice_asm_inc !: Apply the seaice increment 52 53 #if defined key_asminc 66 PUBLIC logchl_asm_inc !: Apply the logchl increment (including balancing) 67 PUBLIC pco2_asm_inc !: Apply the pco2 or fco2 increment (including balancing) 68 54 69 LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .TRUE. !: Logical switch for assimilation increment interface 55 #else 56 LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE. !: No assimilation increments 57 #endif 58 LOGICAL, PUBLIC :: ln_bkgwri = .FALSE. !: No output of the background state fields 59 LOGICAL, PUBLIC :: ln_asmiau = .FALSE. !: No applying forcing with an assimilation increment 60 LOGICAL, PUBLIC :: ln_asmdin = .FALSE. !: No direct initialization 61 LOGICAL, PUBLIC :: ln_trainc = .FALSE. !: No tracer (T and S) assimilation increments 62 LOGICAL, PUBLIC :: ln_dyninc = .FALSE. !: No dynamics (u and v) assimilation increments 63 LOGICAL, PUBLIC :: ln_sshinc = .FALSE. !: No sea surface height assimilation increment 64 LOGICAL, PUBLIC :: ln_seaiceinc !: No sea ice concentration increment 65 LOGICAL, PUBLIC :: ln_salfix = .FALSE. !: Apply minimum salinity check 70 LOGICAL, PUBLIC :: ln_bkgwri = .FALSE. !: No output of the background state fields 71 LOGICAL, PUBLIC :: ln_balwri = .FALSE. !: No output of assimilation balancing increments 72 LOGICAL, PUBLIC :: ln_asmiau = .FALSE. !: No applying forcing with an assimilation increment 73 LOGICAL, PUBLIC :: ln_asmdin = .FALSE. !: No direct initialization 74 LOGICAL, PUBLIC :: ln_trainc = .FALSE. !: No tracer (T and S) assimilation increments 75 LOGICAL, PUBLIC :: ln_dyninc = .FALSE. !: No dynamics (u and v) assimilation increments 76 LOGICAL, PUBLIC :: ln_sshinc = .FALSE. !: No sea surface height assimilation increment 77 LOGICAL, PUBLIC :: ln_diuinc = .FALSE. !: No diurnal SST assimilation 78 LOGICAL, PUBLIC :: ln_seaiceinc = .FALSE. !: No sea ice concentration increment 79 LOGICAL, PUBLIC :: ln_logchlinc = .FALSE. !: No logchl increment 80 LOGICAL, PUBLIC :: ln_pco2inc = .FALSE. !: No pco2 increment 81 LOGICAL, PUBLIC :: ln_fco2inc = .FALSE. !: No fco2 increment 82 LOGICAL, PUBLIC :: ln_salfix = .FALSE. !: Apply minimum salinity check 66 83 LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing 67 INTEGER, PUBLIC :: nn_divdmp 84 INTEGER, PUBLIC :: nn_divdmp = 0 !: Apply divergence damping filter nn_divdmp times 68 85 69 86 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: t_bkg , s_bkg !: Background temperature and salinity … … 72 89 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: u_bkginc, v_bkginc !: Increment to the u- & v-components 73 90 REAL(wp), PUBLIC, DIMENSION(:) , ALLOCATABLE :: wgtiau !: IAU weights for each time step 74 #if defined key_asminc75 91 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ssh_iau !: IAU-weighted sea surface height increment 76 #endif77 92 ! !!! time steps relative to the cycle interval [0,nitend-nit000-1] 78 93 INTEGER , PUBLIC :: nitbkg !: Time step of the background state used in the Jb term … … 85 100 REAL(wp), PUBLIC :: salfixmin !: Ensure that the salinity is larger than this value if (ln_salfix) 86 101 87 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ssh_bkg, ssh_bkginc ! Background sea surface height and its increment 88 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: seaice_bkginc ! Increment to the background sea ice conc 89 102 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ssh_bkg, ssh_bkginc !: Background sea surface height and its increment 103 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: seaice_bkginc !: Increment to the background sea ice conc 104 105 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: logchl_bkginc !: Increment to background logchl 106 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: nutrient_bkginc_logchl !: Increment to nutrient due to logchl 107 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: phyto_bkginc_logchl !: Increment to phytoplankton due to logchl 108 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: zoo_bkginc_logchl !: Increment to zooplankton due to logchl 109 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: dn_bkginc_logchl !: Increment to detritus due to logchl 110 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: tco2_bkginc_logchl !: Increment to DIC due to logchl 111 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: alk_bkginc_logchl !: Increment to alkalinity due to logchl 112 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: chl_bkg !: Background chlorophyll (non-log) 113 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: cchl_p_bkg !: Background carbon to chlorophyll ratio 114 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: pgrow_avg_bkg !: Background phytoplankton growth 115 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ploss_avg_bkg !: Background phytoplankton loss 116 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: phyt_avg_bkg !: Background phytoplankton 117 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: mld_max_bkg !: Background maximum mixed layer depth 118 119 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: pco2_bkginc !: Increment to background pco2 120 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: fco2_bkginc !: Increment to background fco2 121 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: tco2_bkginc_pco2 !: Increment to DIC due to pco2/fco2 122 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: alk_bkginc_pco2 !: Increment to alkalinity due to pco2/fco2 123 124 INTEGER :: mld_choice = 4 !: choice of mld criteria to use for physics assimilation 125 !: 1) hmld - Turbocline/mixing depth [W points] 126 !: 2) hmlp - Density criterion (0.01 kg/m^3 change from 10m) [W points] 127 !: 3) hmld_kara - Kara MLD [Interpolated] 128 !: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points] 129 130 INTEGER :: mld_choice_logchl = 5 !: choice of mld criteria to use for logchl assimilation 131 !: 1) hmld - Turbocline/mixing depth [W points] 132 !: 2) hmlp - Density criterion (0.01 kg/m^3 change from 10m) [W points] 133 !: 3) hmld_kara - Kara MLD [Interpolated] 134 !: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points] 135 !: 5) hmlpt - Density criterion (0.01 kg/m^3 change from 10m) [T points] 136 137 INTEGER :: mld_choice_pco2 = 5 !: choice of mld criteria to use for pco2/fco2 assimilation 138 !: 1) hmld - Turbocline/mixing depth [W points] 139 !: 2) hmlp - Density criterion (0.01 kg/m^3 change from 10m) [W points] 140 !: 3) hmld_kara - Kara MLD [Interpolated] 141 !: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points] 142 !: 5) hmlpt - Density criterion (0.01 kg/m^3 change from 10m) [T points] 143 90 144 !! * Substitutions 91 # include "domzgr_substitute.h90" 92 # include "ldfdyn_substitute.h90" 93 # include "vectopt_loop_substitute.h90" 145 !!---------------------------------------------------------------------- 146 !! *** domzgr_substitute.h90 *** 147 !!---------------------------------------------------------------------- 148 !! ** purpose : substitute fsdep. and fse.., the vert. depth and scale 149 !! factors depending on the vertical coord. used, using CPP macro. 150 !!---------------------------------------------------------------------- 151 !! History : 1.0 ! 2005-10 (A. Beckmann, G. Madec) generalisation to all coord. 152 !! 3.1 ! 2009-02 (G. Madec, M. Leclair) pure z* coordinate 153 !!---------------------------------------------------------------------- 154 ! reference for s- or zps-coordinate (3D no time dependency) 155 ! s* or z*-coordinate (3D + time dependency) + use of additional now arrays (..._n) 156 157 158 159 160 161 !!---------------------------------------------------------------------- 162 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 163 !! $Id$ 164 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 165 !!---------------------------------------------------------------------- 166 !!---------------------------------------------------------------------- 167 !! *** ldfdyn_substitute.h90 *** 168 !!---------------------------------------------------------------------- 169 !! ** purpose : substitute fsahm., the lateral eddy viscosity coeff. 170 !! with a constant, or 1D, or 2D or 3D array, using CPP macro. 171 !!---------------------------------------------------------------------- 172 !!---------------------------------------------------------------------- 173 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 174 !! $Id$ 175 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 176 !!---------------------------------------------------------------------- 177 !! 178 !! fsahmt, fsahmf - used for laplaian operator only 179 !! fsahmu, fsahmv - used for bilaplacian operator only 180 !! 181 ! default option : Constant coefficient 182 !!---------------------------------------------------------------------- 183 !! *** vectopt_loop_substitute *** 184 !!---------------------------------------------------------------------- 185 !! ** purpose : substitute the inner loop starting and inding indices 186 !! to allow unrolling of do-loop using CPP macro. 187 !!---------------------------------------------------------------------- 188 !!---------------------------------------------------------------------- 189 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 190 !! $Id$ 191 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 192 !!---------------------------------------------------------------------- 193 94 194 !!---------------------------------------------------------------------- 95 195 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 109 209 !! ** Action : 110 210 !!---------------------------------------------------------------------- 111 INTEGER :: ji, jj, jk, jt ! dummy loop indices 112 INTEGER :: imid, inum ! local integers 113 INTEGER :: ios ! Local integer output status for namelist read 211 !! 212 !! 213 INTEGER :: ji,jj,jk 214 INTEGER :: jt 215 INTEGER :: imid 216 INTEGER :: inum 114 217 INTEGER :: iiauper ! Number of time steps in the IAU period 115 218 INTEGER :: icycper ! Number of time steps in the cycle … … 119 222 INTEGER :: iitiaustr_date ! Date YYYYMMDD of IAU interval start time step 120 223 INTEGER :: iitiaufin_date ! Date YYYYMMDD of IAU interval final time step 121 ! 224 122 225 REAL(wp) :: znorm ! Normalization factor for IAU weights 123 REAL(wp) :: ztotwgt ! Value of time-integrated IAU weights (should be equal to one) 226 REAL(wp) :: ztotwgt ! Value of time-integrated IAU weights 227 ! ! (should be equal to one) 124 228 REAL(wp) :: z_inc_dateb ! Start date of interval on which increment is valid 125 229 REAL(wp) :: z_inc_datef ! End date of interval on which increment is valid 126 230 REAL(wp) :: zdate_bkg ! Date in background state file for DI 127 231 REAL(wp) :: zdate_inc ! Time axis in increments file 128 ! 129 REAL(wp), POINTER, DIMENSION(:,:) :: hdiv ! 2D workspace 130 !! 131 NAMELIST/nam_asminc/ ln_bkgwri, & 232 233 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: & 234 & t_bkginc_2d ! file for reading in 2D 235 ! ! temperature increments 236 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: & 237 & z_mld ! Mixed layer depth 238 239 REAL(wp), POINTER, DIMENSION(:,:) :: hdiv 240 241 !! 242 NAMELIST/nam_asminc/ ln_bkgwri, ln_balwri, & 132 243 & ln_trainc, ln_dyninc, ln_sshinc, & 133 & ln_asmdin, ln_asmiau, & 244 & ln_seaiceinc, ln_logchlinc, ln_pco2inc, & 245 & ln_fco2inc, ln_asmdin, ln_asmiau, & 134 246 & nitbkg, nitdin, nitiaustr, nitiaufin, niaufn, & 135 & ln_salfix, salfixmin, nn_divdmp 247 & nittrjfrq, ln_salfix, salfixmin, & 248 & nn_divdmp, mld_choice, mld_choice_logchl, & 249 & mld_choice_pco2 250 136 251 !!---------------------------------------------------------------------- 137 252 … … 139 254 ! Read Namelist nam_asminc : assimilation increment interface 140 255 !----------------------------------------------------------------------- 256 257 ! Set default values 258 ln_bkgwri = .FALSE. 259 ln_balwri = .FALSE. 260 ln_trainc = .FALSE. 261 ln_dyninc = .FALSE. 262 ln_sshinc = .FALSE. 141 263 ln_seaiceinc = .FALSE. 264 ln_logchlinc = .FALSE. 265 ln_pco2inc = .FALSE. 266 ln_fco2inc = .FALSE. 267 ln_asmdin = .FALSE. 268 ln_asmiau = .TRUE. 269 ln_salfix = .FALSE. 142 270 ln_temnofreeze = .FALSE. 143 144 REWIND( numnam_ref ) ! Namelist nam_asminc in reference namelist : Assimilation increment 145 READ ( numnam_ref, nam_asminc, IOSTAT = ios, ERR = 901) 146 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_asminc in reference namelist', lwp ) 147 148 REWIND( numnam_cfg ) ! Namelist nam_asminc in configuration namelist : Assimilation increment 149 READ ( numnam_cfg, nam_asminc, IOSTAT = ios, ERR = 902 ) 150 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_asminc in configuration namelist', lwp ) 151 IF(lwm) WRITE ( numond, nam_asminc ) 152 271 salfixmin = -9999 272 nitbkg = 0 273 nitdin = 0 274 nitiaustr = 1 275 nitiaufin = 150 ! = 10 days with ORCA2 276 niaufn = 0 277 nittrjfrq = 1 278 279 REWIND ( numnam ) 280 READ ( numnam, nam_asminc ) 281 282 ! Set the data time for diagnostics to the end of the IAU period 283 ! and multiply by the timestep to get seconds from start of run 284 data_time = rdt * nitiaufin 285 286 IF( ln_sco .AND. (ln_sshinc .OR. ln_seaiceinc .OR. ln_asmdin & 287 & .OR. ln_dyninc .OR. ln_logchlinc .OR. ln_pco2inc & 288 & .OR. ln_fco2inc ) )THEN 289 CALL ctl_warn("Only SST assimilation currently supported in "//& 290 & "s-coordinates") 291 ln_sshinc = .FALSE. 292 ln_seaiceinc = .FALSE. 293 ln_asmdin = .FALSE. 294 ln_dyninc = .FALSE. 295 ln_logchlinc = .FALSE. 296 ln_pco2inc = .FALSE. 297 ln_fco2inc = .FALSE. 298 ENDIF 299 300 IF ( ( ln_balwri ).AND.( .NOT. ( ( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) ) ) ) THEN 301 CALL ctl_warn( ' Balancing increments are only calculated for logchl, pco2, fco2', & 302 & ' Not assimilating any of these, so ln_balwri will be set to .false.') 303 ln_balwri = .FALSE. 304 ENDIF 305 153 306 ! Control print 154 307 IF(lwp) THEN … … 158 311 WRITE(numout,*) ' Namelist namasm : set assimilation increment parameters' 159 312 WRITE(numout,*) ' Logical switch for writing out background state ln_bkgwri = ', ln_bkgwri 313 WRITE(numout,*) ' Logical switch for writing out balancing increments ln_balwri = ', ln_balwri 160 314 WRITE(numout,*) ' Logical switch for applying tracer increments ln_trainc = ', ln_trainc 161 315 WRITE(numout,*) ' Logical switch for applying velocity increments ln_dyninc = ', ln_dyninc … … 163 317 WRITE(numout,*) ' Logical switch for Direct Initialization (DI) ln_asmdin = ', ln_asmdin 164 318 WRITE(numout,*) ' Logical switch for applying sea ice increments ln_seaiceinc = ', ln_seaiceinc 319 WRITE(numout,*) ' Logical switch for applying logchl increments ln_logchlinc = ', ln_logchlinc 320 WRITE(numout,*) ' Logical switch for applying pco2 increments ln_pco2inc = ', ln_pco2inc 321 WRITE(numout,*) ' Logical switch for applying fco2 increments ln_fco2inc = ', ln_fco2inc 165 322 WRITE(numout,*) ' Logical switch for Incremental Analysis Updating (IAU) ln_asmiau = ', ln_asmiau 166 323 WRITE(numout,*) ' Timestep of background in [0,nitend-nit000-1] nitbkg = ', nitbkg … … 169 326 WRITE(numout,*) ' Timestep of end of IAU interval in [0,nitend-nit000-1] nitiaufin = ', nitiaufin 170 327 WRITE(numout,*) ' Type of IAU weighting function niaufn = ', niaufn 328 WRITE(numout,*) ' Frequency of trajectory output for 4D-VAR nittrjfrq = ', nittrjfrq 171 329 WRITE(numout,*) ' Logical switch for ensuring that the sa > salfixmin ln_salfix = ', ln_salfix 172 330 WRITE(numout,*) ' Minimum salinity after applying the increments salfixmin = ', salfixmin 173 ENDIF 174 175 nitbkg_r = nitbkg + nit000 - 1 ! Background time referenced to nit000 176 nitdin_r = nitdin + nit000 - 1 ! Background time for DI referenced to nit000 177 nitiaustr_r = nitiaustr + nit000 - 1 ! Start of IAU interval referenced to nit000 178 nitiaufin_r = nitiaufin + nit000 - 1 ! End of IAU interval referenced to nit000 179 180 iiauper = nitiaufin_r - nitiaustr_r + 1 ! IAU interval length 181 icycper = nitend - nit000 + 1 ! Cycle interval length 182 183 CALL calc_date( nit000, nitend , ndate0, iitend_date ) ! Date of final time step 184 CALL calc_date( nit000, nitbkg_r , ndate0, iitbkg_date ) ! Background time for Jb referenced to ndate0 185 CALL calc_date( nit000, nitdin_r , ndate0, iitdin_date ) ! Background time for DI referenced to ndate0 186 CALL calc_date( nit000, nitiaustr_r, ndate0, iitiaustr_date ) ! IAU start time referenced to ndate0 187 CALL calc_date( nit000, nitiaufin_r, ndate0, iitiaufin_date ) ! IAU end time referenced to ndate0 188 ! 331 WRITE(numout,*) ' Choice of MLD for physics assimilation mld_choice = ', mld_choice 332 WRITE(numout,*) ' Choice of MLD for logchl assimilation mld_choice_logchl = ', mld_choice_logchl 333 WRITE(numout,*) ' Choice of MLD for pCO2/fCO2 assimilation mld_choice_pco2 = ', mld_choice_pco2 334 ENDIF 335 336 nitbkg_r = nitbkg + nit000 - 1 ! Background time referenced to nit000 337 nitdin_r = nitdin + nit000 - 1 ! Background time for DI referenced to nit000 338 nitiaustr_r = nitiaustr + nit000 - 1 ! Start of IAU interval referenced to nit000 339 nitiaufin_r = nitiaufin + nit000 - 1 ! End of IAU interval referenced to nit000 340 341 iiauper = nitiaufin_r - nitiaustr_r + 1 ! IAU interval length 342 icycper = nitend - nit000 + 1 ! Cycle interval length 343 344 ! Date of final time step 345 CALL calc_date( nit000, nitend, ndate0, iitend_date ) 346 347 ! Background time for Jb referenced to ndate0 348 CALL calc_date( nit000, nitbkg_r, ndate0, iitbkg_date ) 349 350 ! Background time for DI referenced to ndate0 351 CALL calc_date( nit000, nitdin_r, ndate0, iitdin_date ) 352 353 ! IAU start time referenced to ndate0 354 CALL calc_date( nit000, nitiaustr_r, ndate0, iitiaustr_date ) 355 356 ! IAU end time referenced to ndate0 357 CALL calc_date( nit000, nitiaufin_r, ndate0, iitiaufin_date ) 358 189 359 IF(lwp) THEN 190 360 WRITE(numout,*) … … 197 367 WRITE(numout,*) ' nitiaustr_r = ', nitiaustr_r 198 368 WRITE(numout,*) ' nitiaufin_r = ', nitiaufin_r 369 WRITE(numout,*) ' nittrjfrq = ', nittrjfrq 199 370 WRITE(numout,*) 200 371 WRITE(numout,*) ' Dates referenced to current cycle:' … … 218 389 219 390 IF ( ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) & 220 .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) & 221 & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', & 391 & .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc ) .OR. & 392 & ( ln_logchlinc ) .OR. ( ln_pco2inc ) .OR. ( ln_fco2inc ) )) & 393 & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', & 394 & ' ln_logchlinc, ln_pco2inc and ln_fco2inc is set to .true.', & 222 395 & ' but ln_asmdin and ln_asmiau are both set to .false. :', & 223 396 & ' Inconsistent options') 397 398 IF ( ( ln_bkgwri ).AND.( ( ln_asmdin ).OR.( ln_asmiau ) ) ) & 399 & CALL ctl_stop( ' ln_bkgwri and either ln_asmdin or ln_asmiau are set to .true.:', & 400 & ' The background state must be written before applying the increments') 224 401 225 402 IF ( ( niaufn /= 0 ).AND.( niaufn /= 1 ) ) & … … 227 404 & ' Type IAU weighting function is invalid') 228 405 229 IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) & 230 & ) & 231 & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', & 406 IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ).AND. & 407 & ( .NOT. ln_logchlinc ).AND.( .NOT. ln_pco2inc ).AND.( .NOT. ln_fco2inc ) ) & 408 & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc, ln_logchlinc, ln_pco2inc, ln_fco2inc', & 409 & ' are set to .false. :', & 232 410 & ' The assimilation increments are not applied') 233 411 … … 250 428 & ' the cycle interval') 251 429 430 IF ( ( ln_pco2inc ).AND.( ln_fco2inc ) ) & 431 & CALL ctl_stop( ' ln_pco2inc and ln_fco2inc are both set :', & 432 & ' Can only assimilate either one or the other') 433 434 IF ( ( ln_logchlinc ).AND.( ( ln_pco2inc ).OR.( ln_fco2inc ) ) ) & 435 & CALL ctl_warn( ' ln_logchlinc and either ln_pco2inc or ln_fco2inc are set :', & 436 & ' Only increments from pCO2/fCO2 will be applied to the DIC and alkalinity', & 437 & ' DIC and alkalinity increments due to chlorophyll will be set to zero', & 438 & ' This means that carbon will not be conserved' ) 439 252 440 IF ( nstop > 0 ) RETURN ! if there are any errors then go no further 253 254 441 !-------------------------------------------------------------------- 255 442 ! Initialize the Incremental Analysis Updating weighting function … … 280 467 ! Compute the normalization factor 281 468 znorm = 0.0 282 IF ( MOD( iiauper, 2 ) == 0 ) THEN ! Even number of time steps in IAU interval469 IF ( MOD( iiauper, 2 ) == 0 ) THEN ! Even number of time steps in IAU interval 283 470 imid = iiauper / 2 284 471 DO jt = 1, imid … … 327 514 !-------------------------------------------------------------------- 328 515 329 ALLOCATE( t_bkginc(jpi,jpj,jpk) ) 330 ALLOCATE( s_bkginc(jpi,jpj,jpk) ) 331 ALLOCATE( u_bkginc(jpi,jpj,jpk) ) 332 ALLOCATE( v_bkginc(jpi,jpj,jpk) ) 333 ALLOCATE( ssh_bkginc(jpi,jpj) ) 334 ALLOCATE( seaice_bkginc(jpi,jpj)) 335 #if defined key_asminc 516 IF ( ln_trainc ) THEN 517 ALLOCATE( t_bkginc(jpi,jpj,jpk) ) 518 ALLOCATE( s_bkginc(jpi,jpj,jpk) ) 519 t_bkginc(:,:,:) = 0.0 520 s_bkginc(:,:,:) = 0.0 521 ENDIF 522 IF ( ln_dyninc ) THEN 523 ALLOCATE( u_bkginc(jpi,jpj,jpk) ) 524 ALLOCATE( v_bkginc(jpi,jpj,jpk) ) 525 u_bkginc(:,:,:) = 0.0 526 v_bkginc(:,:,:) = 0.0 527 ENDIF 528 IF ( ln_sshinc ) THEN 529 ALLOCATE( ssh_bkginc(jpi,jpj) ) 530 ssh_bkginc(:,:) = 0.0 531 ENDIF 532 IF ( ln_seaiceinc ) THEN 533 ALLOCATE( seaice_bkginc(jpi,jpj)) 534 seaice_bkginc(:,:) = 0.0 535 ENDIF 536 IF ( ln_logchlinc ) THEN 537 ALLOCATE( logchl_bkginc(jpi,jpj) ) 538 ALLOCATE( nutrient_bkginc_logchl(jpi,jpj,jpk) ) 539 ALLOCATE( phyto_bkginc_logchl(jpi,jpj,jpk) ) 540 ALLOCATE( zoo_bkginc_logchl(jpi,jpj,jpk) ) 541 ALLOCATE( dn_bkginc_logchl(jpi,jpj,jpk) ) 542 ALLOCATE( tco2_bkginc_logchl(jpi,jpj,jpk) ) 543 ALLOCATE( alk_bkginc_logchl(jpi,jpj,jpk) ) 544 logchl_bkginc(:,:) = 0.0 545 nutrient_bkginc_logchl(:,:,:) = 0.0 546 phyto_bkginc_logchl(:,:,:) = 0.0 547 zoo_bkginc_logchl(:,:,:) = 0.0 548 dn_bkginc_logchl(:,:,:) = 0.0 549 tco2_bkginc_logchl(:,:,:) = 0.0 550 alk_bkginc_logchl(:,:,:) = 0.0 551 ENDIF 552 IF ( ln_pco2inc ) THEN 553 ALLOCATE( pco2_bkginc(jpi,jpj) ) 554 pco2_bkginc(:,:) = 0.0 555 ENDIF 556 IF ( ln_fco2inc ) THEN 557 ALLOCATE( fco2_bkginc(jpi,jpj) ) 558 fco2_bkginc(:,:) = 0.0 559 ENDIF 560 IF ( ( ln_pco2inc ).OR.( ln_fco2inc ) ) THEN 561 ALLOCATE( tco2_bkginc_pco2(jpi,jpj) ) 562 ALLOCATE( alk_bkginc_pco2(jpi,jpj) ) 563 tco2_bkginc_pco2(:,:) = 0.0 564 alk_bkginc_pco2(:,:) = 0.0 565 ENDIF 336 566 ALLOCATE( ssh_iau(jpi,jpj) ) 337 #endif338 t_bkginc(:,:,:) = 0.0339 s_bkginc(:,:,:) = 0.0340 u_bkginc(:,:,:) = 0.0341 v_bkginc(:,:,:) = 0.0342 ssh_bkginc(:,:) = 0.0343 seaice_bkginc(:,:) = 0.0344 #if defined key_asminc345 567 ssh_iau(:,:) = 0.0 346 #endif 347 IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN 348 568 IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) & 569 & .OR.( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) ) THEN 349 570 !-------------------------------------------------------------------- 350 571 ! Read the increments from file … … 378 599 379 600 IF ( ln_trainc ) THEN 380 CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 ) 381 CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 ) 382 ! Apply the masks 383 t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:) 384 s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:) 385 ! Set missing increments to 0.0 rather than 1e+20 386 ! to allow for differences in masks 387 WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0 388 WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0 601 602 IF (ln_sco) THEN 603 604 ALLOCATE(z_mld(jpi,jpj)) 605 SELECT CASE(mld_choice) 606 CASE(1) 607 z_mld = hmld 608 CASE(2) 609 z_mld = hmlp 610 CASE(3) 611 CALL ctl_stop("Kara mixed layer not defined in current version of NEMO") ! JW: Safety feature, should be removed 612 ! once the Kara mixed layer is available 613 CASE(4) 614 z_mld = hmld_tref 615 END SELECT 616 617 ALLOCATE( t_bkginc_2d(jpi,jpj) ) 618 CALL iom_get( inum, jpdom_autoglo, 'bckinsurft', t_bkginc_2d, 1) 619 DO jk = 1,jpkm1 620 WHERE( z_mld(:,:) > gdepw_1(:,:,jk) ) 621 t_bkginc(:,:,jk) = t_bkginc_2d(:,:) * 0.5 * & 622 & ( 1 + cos( (gdept_1(:,:,jk)/z_mld(:,:) ) * rpi ) ) 623 624 t_bkginc(:,:,jk) = t_bkginc(:,:,jk) * bdytmask(:,:) 625 ELSEWHERE 626 t_bkginc(:,:,jk) = 0. 627 ENDWHERE 628 ENDDO 629 s_bkginc(:,:,:) = 0. 630 631 ! DEALLOCATE temporary arrays 632 DEALLOCATE(z_mld, t_bkginc_2d) 633 634 ELSE 635 636 CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 ) 637 CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 ) 638 ! Apply the masks 639 t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:) 640 s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:) 641 ! Set missing increments to 0.0 rather than 1e+20 642 ! to allow for differences in masks 643 WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0 644 WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0 645 646 ENDIF 647 389 648 ENDIF 390 649 … … 419 678 ENDIF 420 679 680 IF ( ln_logchlinc ) THEN 681 CALL iom_get( inum, jpdom_autoglo, 'bckinlogchl', logchl_bkginc, 1 ) 682 ! Apply the masks 683 logchl_bkginc(:,:) = logchl_bkginc(:,:) * tmask(:,:,1) 684 ! Set missing increments to 0.0 rather than 1e+20 685 ! to allow for differences in masks 686 WHERE( ABS( logchl_bkginc(:,:) ) > 1.0e+10 ) logchl_bkginc(:,:) = 0.0 687 ENDIF 688 689 IF ( ln_pco2inc ) THEN 690 CALL iom_get( inum, jpdom_autoglo, 'bckinpco2', pco2_bkginc, 1 ) 691 ! Apply the masks 692 pco2_bkginc(:,:) = pco2_bkginc(:,:) * tmask(:,:,1) 693 ! Set missing increments to 0.0 rather than 1e+20 694 ! to allow for differences in masks 695 WHERE( ABS( pco2_bkginc(:,:) ) > 1.0e+10 ) pco2_bkginc(:,:) = 0.0 696 ENDIF 697 698 IF ( ln_fco2inc ) THEN 699 CALL iom_get( inum, jpdom_autoglo, 'bckinfco2', fco2_bkginc, 1 ) 700 ! Apply the masks 701 fco2_bkginc(:,:) = fco2_bkginc(:,:) * tmask(:,:,1) 702 ! Set missing increments to 0.0 rather than 1e+20 703 ! to allow for differences in masks 704 WHERE( ABS( fco2_bkginc(:,:) ) > 1.0e+10 ) fco2_bkginc(:,:) = 0.0 705 ENDIF 706 421 707 CALL iom_close( inum ) 422 708 … … 438 724 439 725 DO jj = 2, jpjm1 440 DO ji = fs_2, fs_jpim1 ! vector opt.726 DO ji = 2, jpim1 ! vector opt. 441 727 hdiv(ji,jj) = & 442 ( e2u(ji ,jj ) * fse3u(ji ,jj ,jk) * u_bkginc(ji ,jj ,jk) &443 - e2u(ji-1,jj ) * fse3u(ji-1,jj ,jk) * u_bkginc(ji-1,jj ,jk) &444 + e1v(ji ,jj ) * fse3v(ji ,jj ,jk) * v_bkginc(ji ,jj ,jk) &445 - e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) * v_bkginc(ji ,jj-1,jk) ) &446 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )728 ( e2u(ji ,jj ) * e3u_1(ji ,jj ,jk) * u_bkginc(ji ,jj ,jk) & 729 - e2u(ji-1,jj ) * e3u_1(ji-1,jj ,jk) * u_bkginc(ji-1,jj ,jk) & 730 + e1v(ji ,jj ) * e3v_1(ji ,jj ,jk) * v_bkginc(ji ,jj ,jk) & 731 - e1v(ji ,jj-1) * e3v_1(ji ,jj-1,jk) * v_bkginc(ji ,jj-1,jk) ) & 732 / ( e1t(ji,jj) * e2t(ji,jj) * e3t_1(ji,jj,jk) ) 447 733 END DO 448 734 END DO … … 451 737 452 738 DO jj = 2, jpjm1 453 DO ji = fs_2, fs_jpim1 ! vector opt.739 DO ji = 2, jpim1 ! vector opt. 454 740 u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji+1,jj)*e2t(ji+1,jj) * hdiv(ji+1,jj) & 455 741 - e1t(ji ,jj)*e2t(ji ,jj) * hdiv(ji ,jj) ) & … … 473 759 !----------------------------------------------------------------------- 474 760 ! Allocate and initialize the background state arrays 761 ! For logchl this is required for both direct initialisation and IAU, 762 ! as this information is needed for the nitrogen balancing 475 763 !----------------------------------------------------------------------- 476 764 … … 529 817 CALL iom_close( inum ) 530 818 819 ENDIF 820 821 IF ( ln_logchlinc ) THEN 822 823 ALLOCATE( chl_bkg(jpi,jpj) ) 824 ALLOCATE( cchl_p_bkg(jpi,jpj) ) 825 ALLOCATE( pgrow_avg_bkg(jpi,jpj) ) 826 ALLOCATE( ploss_avg_bkg(jpi,jpj) ) 827 ALLOCATE( phyt_avg_bkg(jpi,jpj) ) 828 ALLOCATE( mld_max_bkg(jpi,jpj) ) 829 830 ! Just use initial model state rather than reading in from background file 831 ! Keep with the separate arrays as in future we may want to save a daily mean 832 ! from the obs operator step and use that instead, as was done at v3.0/v3.2 531 833 ENDIF 532 834 ! … … 650 952 INTEGER :: ji,jj,jk 651 953 INTEGER :: it 652 REAL(wp) :: zincwgt ! IAU weight for current time step653 REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values954 REAL(wp) :: zincwgt ! IAU weight for current time step 955 REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values 654 956 !!---------------------------------------------------------------------- 655 957 … … 657 959 ! used to prevent the applied increments taking the temperature below the local freezing point 658 960 659 DO jk = 1, jpkm1 660 CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), fsdept(:,:,jk) ) 661 END DO 961 ! Note: For NEMO/CICE this will be identical to nTf (for the surface), but defined at the now point. 962 963 DO jk=1, jpkm1 964 fzptnz (:,:,jk) = tfreez(tsn(:,:,jk,jp_sal )) - 7.53e-4_wp * gdepw_1(:,:,jk) 965 ENDDO 662 966 663 967 IF ( ln_asmiau ) THEN … … 674 978 IF(lwp) THEN 675 979 WRITE(numout,*) 676 WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it) 980 WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', & 981 & kt,' with IAU weight = ', wgtiau(it) 677 982 WRITE(numout,*) '~~~~~~~~~~~~' 678 983 ENDIF … … 708 1013 ENDIF 709 1014 710 711 1015 ELSEIF ( ln_asmdin ) THEN 712 1016 … … 717 1021 IF ( kt == nitdin_r ) THEN 718 1022 719 neuler = 0 ! Force Euler forward step1023 neuler = 0 ! Force Euler forward step 720 1024 721 1025 ! Initialize the now fields with the background + increment 722 1026 IF (ln_temnofreeze) THEN 723 1027 ! Do not apply negative increments if the temperature will fall below freezing 724 WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 1028 WHERE(t_bkginc(:,:,:) > 0.0_wp .OR. & 1029 & tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 725 1030 tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:) 726 1031 END WHERE … … 731 1036 ! Do not apply negative increments if the salinity will fall below a specified 732 1037 ! minimum value salfixmin 733 WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 1038 WHERE(s_bkginc(:,:,:) > 0.0_wp .OR. & 1039 & tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 734 1040 tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:) 735 1041 END WHERE … … 738 1044 ENDIF 739 1045 740 tsb(:,:,:,:) = tsn(:,:,:,:) ! Update before fields 741 742 CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) ) ! Before potential and in situ densities 743 !!gm fabien 744 ! CALL eos( tsb, rhd, rhop ) ! Before potential and in situ densities 745 !!gm 746 747 748 IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav) & 749 & CALL zps_hde ( kt, jpts, tsb, gtsu, gtsv, & ! Partial steps: before horizontal gradient 750 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 751 IF( ln_zps .AND. .NOT. lk_c1d .AND. ln_isfcav) & 752 & CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv, & ! Partial steps for top cell (ISF) 753 & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , & 754 & gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi ) ! of t, s, rd at the last ocean level 755 756 #if defined key_zdfkpp 757 CALL eos( tsn, rhd, fsdept_n(:,:,:) ) ! Compute rhd 758 !!gm fabien CALL eos( tsn, rhd ) ! Compute rhd 759 #endif 1046 tsb(:,:,:,:) = tsn(:,:,:,:) ! Update before fields 1047 1048 CALL eos( tsb, rhd, rhop ) ! Before potential and in situ densities 1049 1050 IF( ln_zps .AND. .NOT. lk_c1d ) & 1051 & CALL zps_hde( nit000, jpts, tsb, & ! Partial steps: before horizontal derivative 1052 & gtsu, gtsv, rhd, & ! of T, S, rd at the bottom ocean level 1053 & gru , grv ) 1054 760 1055 761 1056 DEALLOCATE( t_bkginc ) … … 766 1061 ! 767 1062 ENDIF 1063 1064 1065 768 1066 ! Perhaps the following call should be in step 769 1067 IF ( ln_seaiceinc ) CALL seaice_asm_inc ( kt ) ! apply sea ice concentration increment … … 786 1084 INTEGER :: jk 787 1085 INTEGER :: it 788 REAL(wp) :: zincwgt ! IAU weight for current time step1086 REAL(wp) :: zincwgt ! IAU weight for current time step 789 1087 !!---------------------------------------------------------------------- 790 1088 … … 862 1160 INTEGER :: it 863 1161 INTEGER :: jk 864 REAL(wp) :: zincwgt ! IAU weight for current time step1162 REAL(wp) :: zincwgt ! IAU weight for current time step 865 1163 !!---------------------------------------------------------------------- 866 1164 … … 885 1183 ! Save the tendency associated with the IAU weighted SSH increment 886 1184 ! (applied in dynspg.*) 887 #if defined key_asminc888 1185 ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt 889 #endif890 1186 IF ( kt == nitiaufin_r ) THEN 891 1187 DEALLOCATE( ssh_bkginc ) … … 908 1204 909 1205 ! Update before fields 910 sshb(:,:) = sshn(:,:) 911 912 IF( lk_vvl ) THEN 913 DO jk = 1, jpk 914 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 915 END DO 916 ENDIF 1206 sshb(:,:) = sshn(:,:) 917 1207 918 1208 DEALLOCATE( ssh_bkg ) … … 924 1214 ! 925 1215 END SUBROUTINE ssh_asm_inc 926 927 1216 928 1217 SUBROUTINE seaice_asm_inc( kt, kindic ) … … 936 1225 !! ** Action : 937 1226 !! 938 !!---------------------------------------------------------------------- 1227 !! History : 1228 !! ! 07-2011 (D. Lea) Initial version based on ssh_asm_inc 1229 !!---------------------------------------------------------------------- 1230 939 1231 IMPLICIT NONE 940 ! 941 INTEGER, INTENT(in) :: kt ! Current time step 942 INTEGER, INTENT(in), OPTIONAL :: kindic ! flag for disabling the deallocation 943 ! 944 INTEGER :: it 945 REAL(wp) :: zincwgt ! IAU weight for current time step 946 #if defined key_lim2 947 REAL(wp), DIMENSION(jpi,jpj) :: zofrld, zohicif, zseaicendg, zhicifinc ! LIM 948 REAL(wp) :: zhicifmin = 0.5_wp ! ice minimum depth in metres 949 #endif 950 !!---------------------------------------------------------------------- 1232 1233 !! * Arguments 1234 INTEGER, INTENT(IN) :: kt ! Current time step 1235 INTEGER, OPTIONAL, INTENT(IN) :: kindic ! flag for disabling the deallocation 1236 1237 !! * Local declarations 1238 INTEGER :: it 1239 REAL(wp) :: zincwgt ! IAU weight for current time step 1240 1241 951 1242 952 1243 IF ( ln_asmiau ) THEN … … 969 1260 ENDIF 970 1261 971 ! Sea-ice : LIM-3 case (to add) 972 973 #if defined key_lim2 974 ! Sea-ice : LIM-2 case 975 zofrld (:,:) = frld(:,:) 976 zohicif(:,:) = hicif(:,:) 977 ! 978 frld = MIN( MAX( frld (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 979 pfrld = MIN( MAX( pfrld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp) 980 fr_i(:,:) = 1.0_wp - frld(:,:) ! adjust ice fraction 981 ! 982 zseaicendg(:,:) = zofrld(:,:) - frld(:,:) ! find out actual sea ice nudge applied 983 ! 984 ! Nudge sea ice depth to bring it up to a required minimum depth 985 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 986 zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt 987 ELSEWHERE 988 zhicifinc(:,:) = 0.0_wp 989 END WHERE 990 ! 991 ! nudge ice depth 992 hicif (:,:) = hicif (:,:) + zhicifinc(:,:) 993 phicif(:,:) = phicif(:,:) + zhicifinc(:,:) 994 ! 995 ! seaice salinity balancing (to add) 996 #endif 997 998 #if defined key_cice && defined key_asminc 999 ! Sea-ice : CICE case. Pass ice increment tendency into CICE 1000 ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt 1001 #endif 1262 1002 1263 1003 1264 IF ( kt == nitiaufin_r ) THEN … … 1007 1268 ELSE 1008 1269 1009 #if defined key_cice && defined key_asminc1010 ! Sea-ice : CICE case. Zero ice increment tendency into CICE1011 ndaice_da(:,:) = 0.0_wp1012 #endif1013 1270 1014 1271 ENDIF … … 1024 1281 neuler = 0 ! Force Euler forward step 1025 1282 1026 ! Sea-ice : LIM-3 case (to add)1027 1028 #if defined key_lim21029 ! Sea-ice : LIM-2 case.1030 zofrld(:,:)=frld(:,:)1031 zohicif(:,:)=hicif(:,:)1032 !1033 ! Initialize the now fields the background + increment1034 frld (:,:) = MIN( MAX( frld(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp)1035 pfrld(:,:) = frld(:,:)1036 fr_i (:,:) = 1.0_wp - frld(:,:) ! adjust ice fraction1037 zseaicendg(:,:) = zofrld(:,:) - frld(:,:) ! find out actual sea ice nudge applied1038 !1039 ! Nudge sea ice depth to bring it up to a required minimum depth1040 WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin )1041 zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt1042 ELSEWHERE1043 zhicifinc(:,:) = 0.0_wp1044 END WHERE1045 !1046 ! nudge ice depth1047 hicif (:,:) = hicif (:,:) + zhicifinc(:,:)1048 phicif(:,:) = phicif(:,:)1049 !1050 ! seaice salinity balancing (to add)1051 #endif1052 1283 1053 #if defined key_cice && defined key_asminc1054 ! Sea-ice : CICE case. Pass ice increment tendency into CICE1055 ndaice_da(:,:) = seaice_bkginc(:,:) / rdt1056 #endif1057 1284 IF ( .NOT. PRESENT(kindic) ) THEN 1058 1285 DEALLOCATE( seaice_bkginc ) … … 1061 1288 ELSE 1062 1289 1063 #if defined key_cice && defined key_asminc1064 ! Sea-ice : CICE case. Zero ice increment tendency into CICE1065 ndaice_da(:,:) = 0.0_wp1066 #endif1067 1290 1068 1291 ENDIF 1069 1292 1070 !#if defined defined key_lim2 || defined key_cice 1071 ! 1072 ! IF (ln_seaicebal ) THEN 1073 ! !! balancing salinity increments 1074 ! !! simple case from limflx.F90 (doesn't include a mass flux) 1075 ! !! assumption is that as ice concentration is reduced or increased 1076 ! !! the snow and ice depths remain constant 1077 ! !! note that snow is being created where ice concentration is being increased 1078 ! !! - could be more sophisticated and 1079 ! !! not do this (but would need to alter h_snow) 1080 ! 1081 ! usave(:,:,:)=sb(:,:,:) ! use array as a temporary store 1082 ! 1083 ! DO jj = 1, jpj 1084 ! DO ji = 1, jpi 1085 ! ! calculate change in ice and snow mass per unit area 1086 ! ! positive values imply adding salt to the ocean (results from ice formation) 1087 ! ! fwf : ice formation and melting 1088 ! 1089 ! zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rdt 1090 ! 1091 ! ! change salinity down to mixed layer depth 1092 ! mld=hmld_kara(ji,jj) 1093 ! 1094 ! ! prevent small mld 1095 ! ! less than 10m can cause salinity instability 1096 ! IF (mld < 10) mld=10 1097 ! 1098 ! ! set to bottom of a level 1099 ! DO jk = jpk-1, 2, -1 1100 ! IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN 1101 ! mld=gdepw(ji,jj,jk+1) 1102 ! jkmax=jk 1103 ! ENDIF 1104 ! ENDDO 1105 ! 1106 ! ! avoid applying salinity balancing in shallow water or on land 1107 ! ! 1108 ! 1109 ! ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m) 1110 ! 1111 ! dsal_ocn=0.0_wp 1112 ! sal_thresh=5.0_wp ! minimum salinity threshold for salinity balancing 1113 ! 1114 ! if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) & 1115 ! dsal_ocn = zfons / (rhop(ji,jj,1) * mld) 1116 ! 1117 ! ! put increments in for levels in the mixed layer 1118 ! ! but prevent salinity below a threshold value 1119 ! 1120 ! DO jk = 1, jkmax 1121 ! 1122 ! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 1123 ! sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn 1124 ! sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn 1125 ! ENDIF 1126 ! 1127 ! ENDDO 1128 ! 1129 ! ! ! salt exchanges at the ice/ocean interface 1130 ! ! zpmess = zfons / rdt_ice ! rdt_ice is ice timestep 1131 ! ! 1132 ! !! Adjust fsalt. A +ve fsalt means adding salt to ocean 1133 ! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt 1134 ! !! 1135 ! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d) 1136 ! !! ! E-P (kg m-2 s-2) 1137 ! ! emp(ji,jj) = emp(ji,jj) + zpmess ! E-P (kg m-2 s-2) 1138 ! ENDDO !ji 1139 ! ENDDO !jj! 1140 ! 1141 ! ENDIF !ln_seaicebal 1142 ! 1143 !#endif 1293 !#if defined key_lim3 || defined key_lim2 || defined key_cice 1294 ! 1295 ! IF (ln_seaicebal ) THEN 1296 ! !! balancing salinity increments 1297 ! !! simple case from limflx.F90 (doesn't include a mass flux) 1298 ! !! assumption is that as ice concentration is reduced or increased 1299 ! !! the snow and ice depths remain constant 1300 ! !! note that snow is being created where ice concentration is being increased 1301 ! !! - could be more sophisticated and 1302 ! !! not do this (but would need to alter h_snow) 1303 ! 1304 ! usave(:,:,:)=sb(:,:,:) ! use array as a temporary store 1305 ! 1306 ! DO jj = 1, jpj 1307 ! DO ji = 1, jpi 1308 ! ! calculate change in ice and snow mass per unit area 1309 ! ! positive values imply adding salt to the ocean (results from ice formation) 1310 ! ! fwf : ice formation and melting 1311 ! 1312 ! zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rdt 1313 ! 1314 ! ! change salinity down to mixed layer depth 1315 ! mld=hmld_kara(ji,jj) 1316 ! 1317 ! ! prevent small mld 1318 ! ! less than 10m can cause salinity instability 1319 ! IF (mld < 10) mld=10 1320 ! 1321 ! ! set to bottom of a level 1322 ! DO jk = jpk-1, 2, -1 1323 ! IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN 1324 ! mld=gdepw(ji,jj,jk+1) 1325 ! jkmax=jk 1326 ! ENDIF 1327 ! ENDDO 1328 ! 1329 ! ! avoid applying salinity balancing in shallow water or on land 1330 ! ! 1331 ! 1332 ! ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m) 1333 ! 1334 ! dsal_ocn=0.0_wp 1335 ! sal_thresh=5.0_wp ! minimum salinity threshold for salinity balancing 1336 ! 1337 ! if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) & 1338 ! dsal_ocn = zfons / (rhop(ji,jj,1) * mld) 1339 ! 1340 ! ! put increments in for levels in the mixed layer 1341 ! ! but prevent salinity below a threshold value 1342 ! 1343 ! DO jk = 1, jkmax 1344 ! 1345 ! IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN 1346 ! sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn 1347 ! sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn 1348 ! ENDIF 1349 ! 1350 ! ENDDO 1351 ! 1352 ! ! ! salt exchanges at the ice/ocean interface 1353 ! ! zpmess = zfons / rdt_ice ! rdt_ice is ice timestep 1354 ! ! 1355 ! !! Adjust fsalt. A +ve fsalt means adding salt to ocean 1356 ! !! fsalt(ji,jj) = fsalt(ji,jj) + zpmess ! adjust fsalt 1357 ! !! 1358 ! !! emps(ji,jj) = emps(ji,jj) + zpmess ! or adjust emps (see icestp1d) 1359 ! !! ! E-P (kg m-2 s-2) 1360 ! ! emp(ji,jj) = emp(ji,jj) + zpmess ! E-P (kg m-2 s-2) 1361 ! ENDDO ! ji 1362 ! ENDDO ! jj ! 1363 ! 1364 ! ENDIF ! ln_seaicebal 1365 ! 1366 !#endif 1367 1144 1368 1145 1369 ENDIF 1146 1370 1147 1371 END SUBROUTINE seaice_asm_inc 1148 1372 1373 1374 SUBROUTINE logchl_asm_inc( kt ) 1375 ! Dummy routine 1376 INTEGER, INTENT(IN) :: kt ! Current time step 1377 WRITE(*,*) 'logchl_asm_inc: You should not have seen this print! error?', kt 1378 CALL ctl_stop( ' logchl_asm_inc should not be being called' ) 1379 END SUBROUTINE logchl_asm_inc 1380 1381 1382 SUBROUTINE pco2_asm_inc( kt ) 1383 ! Dummy routine 1384 INTEGER, INTENT(IN) :: kt ! Current time step 1385 WRITE(*,*) 'pco2_asm_inc: You should not have seen this print! error?', kt 1386 CALL ctl_stop( ' pco2_asm_inc should not be being called' ) 1387 END SUBROUTINE pco2_asm_inc 1388 1149 1389 !!====================================================================== 1150 1390 END MODULE asminc
Note: See TracChangeset
for help on using the changeset viewer.
