Changeset 8438
- Timestamp:
- 2017-08-15T12:21:39+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_v3.6_asm_nemovar_community_logchlbal/NEMOGCM/NEMO/OPA_SRC/ASM/asmlogchlbal_ersem.F90
r7803 r8438 22 22 USE par_trc, ONLY: jptra ! Tracer parameters 23 23 24 24 25 IMPLICIT NONE 25 26 PRIVATE … … 27 28 PUBLIC asm_logchl_bal_ersem 28 29 30 REAL(wp), PARAMETER :: tiny = 1e-20 31 29 32 CONTAINS 30 31 SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, & 32 & k_maxchlinc, logchl_bkginc, logchl_balinc ) 33 !!--------------------------------------------------------------------------- 34 !! *** ROUTINE asm_logchl_bal_ersem *** 33 SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, & 34 & k_maxchlinc, logchl_bkginc, logchl_balinc ) 35 36 37 38 !! INTEGER, INTENT(in ) :: npfts 39 !!REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: inc_log_chlTot 40 41 LOGICAL, INTENT(in ) :: ld_logchlpftinc 42 INTEGER, INTENT(in ) :: npfts 43 INTEGER, INTENT(in ) :: mld_choice_bgc 44 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,npfts) :: logchl_bkginc 45 REAL(wp), INTENT(in ) :: k_maxchlinc 46 REAL(wp), INTENT( out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc 47 35 48 !! 36 !! ** Purpose : calculate increments to ERSEM from logchl increments 37 !! 38 !! ** Method : convert logchl increments to chl increments 39 !! split between the ERSEM PFTs 40 !! spread through the mixed layer 41 !! [forthcoming: calculate increments to nutrients and zooplankton] 42 !! 43 !! ** Action : populate logchl_balinc 44 !! 45 !! References : forthcoming... 46 !!--------------------------------------------------------------------------- 47 !! 48 LOGICAL, INTENT(in ) :: ld_logchlpftinc 49 INTEGER, INTENT(in ) :: npfts 50 INTEGER, INTENT(in ) :: mld_choice_bgc 51 REAL(wp), INTENT(in ) :: k_maxchlinc 52 REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,npfts) :: logchl_bkginc 53 REAL(wp), INTENT( out), DIMENSION(jpi,jpj,jpk,jptra) :: logchl_balinc 54 !! 55 INTEGER :: ji, jj, jk 56 INTEGER :: jkmax 57 REAL(wp) :: chl_tot, chl_inc 58 REAL(wp), DIMENSION(jpi,jpj) :: zmld 59 !!--------------------------------------------------------------------------- 60 61 ! Split surface logchl incs into surface Chl1-4 incs 62 ! 63 ! In order to transform logchl incs to chl incs, need to account for the background, 64 ! cannot simply do 10^logchl_bkginc. Need to: 65 ! 1) Add logchl inc to log10(background) to get log10(analysis) 66 ! 2) Take 10^log10(analysis) to get analysis 67 ! 3) Subtract background from analysis to get chl incs 68 ! If k_maxchlinc > 0 then cap total absolute chlorophyll increment at that value 69 ! 70 ! Only apply increments if all of Chl1-4 background values are > 0 71 ! In theory, they always will be, and if any are not that's a sign 72 ! that something's going wrong which the assimilation might make worse 73 ! 49 50 INTEGER :: ji, jj, jk, jp 51 INTEGER :: jkmax 52 REAL(wp) :: chl_tot, chl_inc 53 REAL(wp), DIMENSION(jpi,jpj) :: zmld 54 REAL(wp), DIMENSION(jpi,jpj,npfts,jptra) :: cc_chl_field !! CROSS-COV BETWEEN SPECIFIC SIZE CLASS AND GIVEN BIOGEOCHEM VAR 55 REAL(wp), DIMENSION(jptra) :: threshold 56 57 58 59 ! STEP 1: TRANSFORM ASSIMILATED LOG-FIELD INCREMENTS INTO INCREMENTS AND PRODUCE THEM FOR THE SIZE-CLASSES IF NEEDED, ALSO SAVE INCREMENTS OF LOGS IN SEPARATE VARIABLES (THIS IS IN CASE ONLY TOTAL CHLOROPHYLL WAS ASSIMILATED) 60 61 IF (lwp) WRITE(numout,*) 'Code has gotten to point started' 62 CALL flush(numout) 63 64 74 65 IF ( ld_logchlpftinc ) THEN 75 66 ! … … 151 142 ENDIF 152 143 153 ! Propagate surface Chl1-4 incs through mixed layer 154 ! First, choose mixed layer definition 155 ! 144 145 146 147 148 149 150 IF (lwp) WRITE(numout,*) 'Code has gotten to point log_inc transformed to inc' 151 CALL flush(numout) 152 ! STEP 2: READ CROSS_COVARIANCES 153 154 CALL read_crosscovariances(cc_chl_field, npfts) 155 156 157 IF (lwp) WRITE(numout,*) 'Code has gotten to point cc read' 158 CALL flush(numout) 159 !! STEP 3: CALCULATE INCREMENTS OF OTHER BIOGEOCHEMICAL VARIABLES 160 161 162 threshold(jp_fabm_m1+jp_fabm_n3n) = 100 163 threshold(jp_fabm_m1+jp_fabm_n4n) = 10 164 threshold(jp_fabm_m1+jp_fabm_n1p) = 15 165 threshold(jp_fabm_m1+jp_fabm_n5s) = 80 166 threshold(jp_fabm_m1+jp_fabm_z4c) = 20 167 threshold(jp_fabm_m1+jp_fabm_z5c) = 20 168 threshold(jp_fabm_m1+jp_fabm_p1c) = 25 169 threshold(jp_fabm_m1+jp_fabm_p2c) = 25 170 threshold(jp_fabm_m1+jp_fabm_p3c) = 25 171 threshold(jp_fabm_m1+jp_fabm_p4c) = 25 172 threshold(jp_fabm_m1+jp_fabm_p1n) = 10 173 threshold(jp_fabm_m1+jp_fabm_p2n) = 10 174 threshold(jp_fabm_m1+jp_fabm_p3n) = 10 175 threshold(jp_fabm_m1+jp_fabm_p4n) = 10 176 threshold(jp_fabm_m1+jp_fabm_p1p) = 10 177 threshold(jp_fabm_m1+jp_fabm_p2p) = 10 178 threshold(jp_fabm_m1+jp_fabm_p3p) = 10 179 threshold(jp_fabm_m1+jp_fabm_p4p) = 10 180 181 182 183 184 DO jj = 1, jpj 185 DO ji = 1, jpi 186 DO jp=1,jptra 187 188 189 !CALL t;ge'e;fdgfgd 190 191 IF ( (jp /= jp_fabm_m1+jp_fabm_chl1) .AND. & 192 & (jp /= jp_fabm_m1+jp_fabm_chl2) .AND. & 193 & (jp /= jp_fabm_m1+jp_fabm_chl3) .AND. & 194 & (jp /= jp_fabm_m1+jp_fabm_chl4)) THEN 195 196 CALL asm_inc_biogeo_ersem(npfts, jp, trn(ji,jj,1,jp), logchl_bkginc(ji,jj,:), cc_chl_field(ji,jj,:,:), threshold(jp), & 197 & logchl_balinc(ji,jj,1,jp)) 198 199 200 END IF 201 202 203 END DO 204 205 END DO 206 207 END DO 208 209 210 211 IF (lwp) WRITE(numout,*) 'Code has gotten to point surface logs calculated' 212 213 214 215 216 217 218 219 220 221 CALL flush(numout) 222 !! STEP 4: PROPAGATE EACH INCREMENT THROUGH THE MIXED LAYER 223 224 156 225 SELECT CASE( mld_choice_bgc ) 157 226 CASE ( 1 ) ! Turbocline/mixing depth [W points] … … 179 248 ! Now set MLD to bottom of a level and propagate incs equally through mixed layer 180 249 ! 181 DO jj = 1, jpj 250 251 DO jj = 1, jpj 182 252 DO ji = 1, jpi 183 253 ! … … 192 262 ! 193 263 DO jk = 2, jkmax 194 logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl1) 195 logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl2) 196 logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl3) 197 logchl_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) = logchl_balinc(ji,jj,1,jp_fabm_m1+jp_fabm_chl4) 264 265 DO jp = 1, jptra 266 267 !IF (lwp) WRITE(numout,*) logchl_balinc_ersem_chl1(ji,jj,1) 268 !CALL flush(numout) 269 logchl_balinc(ji,jj,jk,jp) = logchl_balinc(ji,jj,1,jp) 270 271 END DO 198 272 END DO 199 ! 273 274 200 275 END DO 201 276 END DO 202 277 203 ! Multivariate balancing forthcoming... 204 205 END SUBROUTINE asm_logchl_bal_ersem 278 279 IF (lwp) WRITE(numout,*) 'Code has gotten to point log transformed into inc and propagated' 280 CALL flush(numout) 281 282 END SUBROUTINE asm_logchl_bal_ersem 283 284 285 286 !CONTAINS 287 SUBROUTINE read_crosscovariances(cc_chl_field, npfts) 288 289 !! status = nf90_open(path = filepath, cmode = nf90_nowrite, ncid = ccfile) 290 291 INTEGER, INTENT(in ) :: npfts 292 REAL(wp), INTENT( out), DIMENSION(jpi,jpj,npfts,jptra) :: cc_chl_field !! CROSS-COV BETWEEN SPECIFIC SIZE CLASS AND GIVEN BIOGEOCHEM VAR 293 INTEGER :: ncid 294 INTEGER :: pft 295 CHARACTER(LEN=256) :: path_chl(npfts) 296 CHARACTER(LEN=256) :: id_chl(npfts) 297 298 IF (npfts /= 4) THEN 299 300 IF (lwp) WRITE(numout,*) 'Wrong number of Pfts for reading crosscovs files' 301 CALL flush(numout) 302 303 ELSE 304 305 path_chl(1) = "crosscovs_chl1.nc" 306 path_chl(2) = "crosscovs_chl2.nc" 307 path_chl(3) = "crosscovs_chl3.nc" 308 path_chl(4) = "crosscovs_chl4.nc" 309 310 id_chl(1) = "cc_Chl1" 311 id_chl(2) = "cc_Chl2" 312 id_chl(3) = "cc_Chl3" 313 id_chl(4) = "cc_Chl4" 314 315 316 DO pft = 1, npfts 317 318 CALL iom_open( path_chl(pft), ncid ) 319 CALL iom_get( ncid, jpdom_autoglo, 'cc_NO3', cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_n3n)) 320 CALL iom_get( ncid, jpdom_autoglo, 'cc_NH4', cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_n4n)) 321 CALL iom_get( ncid, jpdom_autoglo, "cc_PO4", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_n1p)) 322 CALL iom_get( ncid, jpdom_autoglo, "cc_SiO", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_n5s)) 323 CALL iom_get( ncid, jpdom_autoglo, "cc_Z4c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_z4c)) 324 CALL iom_get( ncid, jpdom_autoglo, "cc_Z5c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_z5c)) 325 CALL iom_get( ncid, jpdom_autoglo, "cc_P1c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p1c)) 326 CALL iom_get( ncid, jpdom_autoglo, "cc_P2c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p2c)) 327 CALL iom_get( ncid, jpdom_autoglo, "cc_P3c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p3c)) 328 CALL iom_get( ncid, jpdom_autoglo, "cc_P4c", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p4c)) 329 CALL iom_get( ncid, jpdom_autoglo, "cc_P1n", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p1n)) 330 CALL iom_get( ncid, jpdom_autoglo, "cc_P2n", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p2n)) 331 CALL iom_get( ncid, jpdom_autoglo, "cc_P3n", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p3n)) 332 CALL iom_get( ncid, jpdom_autoglo, "cc_P4n", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p4n)) 333 CALL iom_get( ncid, jpdom_autoglo, "cc_P1p", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p1p)) 334 CALL iom_get( ncid, jpdom_autoglo, "cc_P2p", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p2p)) 335 CALL iom_get( ncid, jpdom_autoglo, "cc_P3p", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p3p)) 336 CALL iom_get( ncid, jpdom_autoglo, "cc_P4p", cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_p4p)) 337 CALL iom_get( ncid, jpdom_autoglo, id_chl(pft), cc_chl_field(:,:,pft,jp_fabm_m1+jp_fabm_chl1)) 338 CALL iom_close( ncid ) 339 340 END DO 341 342 END IF 343 344 345 END SUBROUTINE read_crosscovariances 346 347 348 !CONTAINS 349 SUBROUTINE asm_inc_biogeo_ersem(npfts, jp, bgf_value, chl, cc, threshold, output) 350 INTEGER, INTENT(in ) :: npfts 351 INTEGER, INTENT(in ) :: jp 352 REAL, INTENT(in ) :: bgf_value 353 REAL, INTENT(in ), DIMENSION(npfts) :: chl 354 REAL, INTENT(in ), DIMENSION(npfts,jptra) :: cc !! VARIANCE OF CHL SIZE-CLASS 355 REAL, INTENT(in ) :: threshold 356 REAL, INTENT( out) :: output 357 358 359 REAL :: total, weight_1, weight_2, weight_3, weight_4 360 361 362 output = 0.0 363 364 IF ( cc(1,jp)*cc(2,jp)*cc(3,jp)*cc(4,jp)*cc(1, jp_fabm_m1+jp_fabm_chl1)*cc(2, jp_fabm_m1+jp_fabm_chl2)*cc(3, jp_fabm_m1+jp_fabm_chl3)*cc(4, jp_fabm_m1+jp_fabm_chl4) /= 0.0 ) THEN 365 366 total = ABS(cc(1,jp)/cc(1, jp_fabm_m1+jp_fabm_chl1)) + ABS(cc(2,jp)/cc(2, jp_fabm_m1+jp_fabm_chl2)) + ABS(cc(3,jp)/cc(3, jp_fabm_m1+jp_fabm_chl3)) + ABS(cc(4,jp)/cc(4, jp_fabm_m1+jp_fabm_chl4)) 367 weight_1 = ABS(cc(1,jp)/cc(1, jp_fabm_m1+jp_fabm_chl1))/total 368 weight_2 = ABS(cc(2,jp)/cc(2, jp_fabm_m1+jp_fabm_chl2))/total 369 weight_3 = ABS(cc(3,jp)/cc(3, jp_fabm_m1+jp_fabm_chl3))/total 370 weight_4 = ABS(cc(4,jp)/cc(4, jp_fabm_m1+jp_fabm_chl4))/total 371 372 ! PRINT *, 'weights', weight_1, weight_2, weight_3, weight_4 373 374 output = weight_1*(cc(1,jp)/cc(1, jp_fabm_m1+jp_fabm_chl1))*chl(1) + weight_2*(cc(2,jp)/cc(2, jp_fabm_m1+jp_fabm_chl2))*chl(2) + weight_3*(cc(3,jp)/cc(3, jp_fabm_m1+jp_fabm_chl3))*chl(3) + weight_4*(cc(4,jp)/cc(4, jp_fabm_m1+jp_fabm_chl4))*chl(4) 375 376 377 ! FIRST BREAK = DON'T INCREASE VALUE MORE THAN TENFOLD 378 379 output = MIN(output, LOG10(11.0)) 380 381 382 ! SECOND BREAK = IF VALUE LARGER THAN THRESHOLD DO NOT INCREASE IT 383 384 IF (bgf_value > threshold) THEN 385 output = MIN(output, 0.0) 386 ENDIF 387 388 ! THIRD BREAK = IF VALUE MORE THAN FIFTH OF THRESHOLD THAN DON'T INCREASE IT BY MORE THAN AN INCREMENT EQUAL TO TENTH OF THRESHOLD (INCREASE IT SLOWLY) 389 390 IF (bgf_value > 0.2*threshold) THEN 391 output = MIN(output, LOG10(0.1*threshold/bgf_value + 1)) 392 ENDIF 393 394 output = (10**( output ) - 1)*bgf_value 395 396 PRINT *, 'output', output 397 398 399 ENDIF 400 401 402 !! calculate increments for one or all variables 403 404 END SUBROUTINE asm_inc_biogeo_ersem 405 406 407 408 409 206 410 207 411 #else … … 210 414 !!---------------------------------------------------------------------- 211 415 CONTAINS 416 CONTAINS 212 417 SUBROUTINE asm_logchl_bal_ersem( ld_logchlpftinc, npfts, mld_choice_bgc, & 213 418 & k_maxchlinc, logchl_bkginc, logchl_balinc ) … … 220 425 WRITE(*,*) 'asm_logchl_bal_ersem: You should not have seen this print! error?' 221 426 END SUBROUTINE asm_logchl_bal_ersem 427 428 429 430 431 432 433 222 434 #endif 223 224 435 !!====================================================================== 225 436 END MODULE asmlogchlbal_ersem
Note: See TracChangeset
for help on using the changeset viewer.
