Changeset 8495 for branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
- Timestamp:
- 2017-09-04T18:50:10+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc_v2/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r8400 r8495 22 22 !! ssh_asm_inc : Apply the SSH increment 23 23 !! seaice_asm_inc : Apply the seaice increment 24 !! logchl_asm_inc : Apply the logchl increment 25 !! pco2_asm_inc : Apply the pCO2/fCO2 increment 24 26 !!---------------------------------------------------------------------- 25 27 USE wrk_nemo ! Memory Allocation … … 44 46 #endif 45 47 USE sbc_oce ! Surface boundary condition variables. 48 USE zdfmxl, ONLY : & 49 ! & hmld_tref, & 50 #if defined key_karaml 51 & hmld_kara, & 52 & ln_kara, & 53 #endif 54 & hmld, & 55 & hmlp, & 56 & hmlpt 57 #if defined key_top 58 USE trc, ONLY: & 59 & trn, & 60 & trb 61 USE par_trc, ONLY: & 62 & jptra 63 #endif 64 #if defined key_medusa && defined key_foam_medusa 65 USE asmlogchlbal_medusa, ONLY: & 66 & asm_logchl_bal_medusa 67 USE asmpco2bal, ONLY: & 68 & asm_pco2_bal 69 USE par_medusa 70 USE mocsy_mainmod 71 #elif defined key_hadocc 72 USE asmlogchlbal_hadocc, ONLY: & 73 & asm_logchl_bal_hadocc 74 USE asmpco2bal, ONLY: & 75 & asm_pco2_bal 76 USE par_hadocc 77 #endif 46 78 47 79 IMPLICIT NONE … … 54 86 PUBLIC ssh_asm_inc !: Apply the SSH increment 55 87 PUBLIC seaice_asm_inc !: Apply the seaice increment 88 PUBLIC logchl_asm_inc !: Apply the logchl increment 89 PUBLIC pco2_asm_inc !: Apply the pCO2/fCO2 increment 56 90 57 91 #if defined key_asminc … … 61 95 #endif 62 96 LOGICAL, PUBLIC :: ln_bkgwri = .FALSE. !: No output of the background state fields 97 LOGICAL, PUBLIC :: ln_balwri = .FALSE. !: No output of the assimilation balancing increments 63 98 LOGICAL, PUBLIC :: ln_asmiau = .FALSE. !: No applying forcing with an assimilation increment 64 99 LOGICAL, PUBLIC :: ln_asmdin = .FALSE. !: No direct initialization … … 67 102 LOGICAL, PUBLIC :: ln_sshinc = .FALSE. !: No sea surface height assimilation increment 68 103 LOGICAL, PUBLIC :: ln_seaiceinc !: No sea ice concentration increment 104 LOGICAL, PUBLIC :: ln_logchlinc = .FALSE. !: No log10(chlorophyll) increment 105 LOGICAL, PUBLIC :: ln_logchlbal = .FALSE. !: Don't apply multivariate log10(chlorophyll) balancing 106 LOGICAL, PUBLIC :: ln_pco2inc = .FALSE. !: No pCO2 increment 107 LOGICAL, PUBLIC :: ln_fco2inc = .FALSE. !: No fCO2 increment 69 108 LOGICAL, PUBLIC :: ln_salfix = .FALSE. !: Apply minimum salinity check 70 109 LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing … … 91 130 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ssh_bkg, ssh_bkginc ! Background sea surface height and its increment 92 131 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: seaice_bkginc ! Increment to the background sea ice conc 132 133 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: logchl_bkginc !: Increment to background logchl 134 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: pco2_bkginc !: Increment to background pCO2 135 #if defined key_top 136 REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: logchl_balinc !: Increment to BGC variables from logchl assim 137 REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: pco2_balinc !: Increment to BGC variables from pCO2 assim 138 #endif 139 #if defined key_hadocc || (defined key_medusa && defined key_foam_medusa) 140 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pgrow_avg_bkg !: Background phyto growth for logchl balancing 141 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ploss_avg_bkg !: Background phyto loss for logchl balancing 142 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: phyt_avg_bkg !: Background phyto for logchl balancing 143 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: mld_max_bkg !: Background max MLD for logchl balancing 144 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg !: Background tracer state variables 145 #endif 146 #if defined key_hadocc 147 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: chl_bkg !: Background surface chlorophyll 148 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: cchl_p_bkg !: Background surface carbon:chlorophyll 149 #endif 150 REAL(wp) :: rn_maxchlinc = -999.0 !: maximum absolute non-log chlorophyll increment from logchl assimilation 151 !: <= 0 implies no maximum applied (switch turned off) 152 !: > 0 implies maximum absolute chl increment capped at this value 153 154 INTEGER :: mld_choice_bgc = 1 !: choice of mld criteria to use for biogeochemistry assimilation 155 !: 1) hmld - Turbocline/mixing depth [W points] 156 !: 2) hmlp - Density criterion (0.01 kg/m^3 change from 10m) [W points] 157 !: 3) hmld_kara - Kara MLD [Interpolated] 158 !: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points] 159 !: 5) hmlpt - Density criterion (0.01 kg/m^3 change from 10m) [T points] 93 160 94 161 !! * Substitutions … … 133 200 REAL(wp), POINTER, DIMENSION(:,:) :: hdiv ! 2D workspace 134 201 !! 135 NAMELIST/nam_asminc/ ln_bkgwri, 202 NAMELIST/nam_asminc/ ln_bkgwri, ln_balwri, & 136 203 & ln_trainc, ln_dyninc, ln_sshinc, & 204 & ln_logchlinc, ln_logchlbal, & 205 & ln_pco2inc, ln_fco2inc, & 137 206 & ln_asmdin, ln_asmiau, & 138 207 & nitbkg, nitdin, nitiaustr, nitiaufin, niaufn, & 139 208 & ln_salfix, salfixmin, nn_divdmp, & 140 & ln_seaiceinc, ln_temnofreeze 209 & ln_seaiceinc, ln_temnofreeze, & 210 & mld_choice_bgc, rn_maxchlinc 141 211 !!---------------------------------------------------------------------- 142 212 … … 163 233 WRITE(numout,*) ' Namelist namasm : set assimilation increment parameters' 164 234 WRITE(numout,*) ' Logical switch for writing out background state ln_bkgwri = ', ln_bkgwri 235 WRITE(numout,*) ' Logical switch for writing out balancing increments ln_balwri = ', ln_balwri 165 236 WRITE(numout,*) ' Logical switch for applying tracer increments ln_trainc = ', ln_trainc 166 237 WRITE(numout,*) ' Logical switch for applying velocity increments ln_dyninc = ', ln_dyninc … … 168 239 WRITE(numout,*) ' Logical switch for Direct Initialization (DI) ln_asmdin = ', ln_asmdin 169 240 WRITE(numout,*) ' Logical switch for applying sea ice increments ln_seaiceinc = ', ln_seaiceinc 241 WRITE(numout,*) ' Logical switch for applying logchl increments ln_logchlinc = ', ln_logchlinc 242 WRITE(numout,*) ' Logical switch for logchl multivariate balancing ln_logchlbal = ', ln_logchlbal 243 WRITE(numout,*) ' Logical switch for applying pCO2 increments ln_pco2inc = ', ln_pco2inc 244 WRITE(numout,*) ' Logical switch for applying fCO2 increments ln_fco2inc = ', ln_fco2inc 170 245 WRITE(numout,*) ' Logical switch for Incremental Analysis Updating (IAU) ln_asmiau = ', ln_asmiau 171 246 WRITE(numout,*) ' Timestep of background in [0,nitend-nit000-1] nitbkg = ', nitbkg … … 176 251 WRITE(numout,*) ' Logical switch for ensuring that the sa > salfixmin ln_salfix = ', ln_salfix 177 252 WRITE(numout,*) ' Minimum salinity after applying the increments salfixmin = ', salfixmin 253 WRITE(numout,*) ' Choice of MLD for BGC assimilation mld_choice_bgc = ', mld_choice_bgc 254 WRITE(numout,*) ' Maximum absolute chlorophyll increment (<=0 = off) rn_maxchlinc = ', rn_maxchlinc 178 255 ENDIF 179 256 … … 223 300 224 301 IF ( ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) & 225 .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) & 226 & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', & 302 & .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ).OR. & 303 & ( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) )) & 304 & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', & 305 & ' ln_logchlinc, ln_pco2inc, and ln_fco2inc is set to .true.', & 227 306 & ' but ln_asmdin and ln_asmiau are both set to .false. :', & 228 307 & ' Inconsistent options') … … 233 312 234 313 IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) & 235 & ) & 236 & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', & 314 & .AND.( .NOT. ln_logchlinc ).AND.( .NOT. ln_pco2inc ).AND.( .NOT. ln_fco2inc ) ) & 315 & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', & 316 & ' ln_logchlinc, ln_pco2inc, and ln_fco2inc are set to .false. :', & 237 317 & ' The assimilation increments are not applied') 238 318 … … 254 334 & ' Background time step for Direct Initialization is outside', & 255 335 & ' the cycle interval') 336 337 IF ( ( ln_balwri ).AND.( .NOT. ln_logchlinc ).AND.( .NOT. ln_pco2inc ).AND.( .NOT. ln_fco2inc ) ) THEN 338 CALL ctl_warn( ' Balancing increments are only calculated for logchl and pCO2/fCO2', & 339 & ' Not assimilating logchl, pCO2 or fCO2, so ln_balwri will be set to .false.') 340 ln_balwri = .FALSE. 341 ENDIF 342 343 IF ( ( ln_logchlbal ).AND.( .NOT. ln_logchlinc ) ) THEN 344 CALL ctl_warn( ' Cannot calculate logchl balancing increments if not assimilating logchl', & 345 & ' Not assimilating logchl, so ln_logchlbal will be set to .false.') 346 ln_logchlbal = .FALSE. 347 ENDIF 348 349 IF ( ( ln_pco2inc ).AND.( ln_fco2inc ) ) THEN 350 CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' ) 351 ENDIF 256 352 257 353 IF ( nstop > 0 ) RETURN ! if there are any errors then go no further … … 350 446 ssh_iau(:,:) = 0.0 351 447 #endif 352 IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN 448 IF ( ln_logchlinc ) THEN 449 ALLOCATE( logchl_bkginc(jpi,jpj) ) 450 logchl_bkginc(:,:) = 0.0 451 #if defined key_top 452 ALLOCATE( logchl_balinc(jpi,jpj,jpk,jptra) ) 453 logchl_balinc(:,:,:,:) = 0.0 454 #endif 455 ENDIF 456 IF ( ln_pco2inc .OR. ln_fco2inc ) THEN 457 ALLOCATE( pco2_bkginc(jpi,jpj) ) 458 pco2_bkginc(:,:) = 0.0 459 #if defined key_top 460 ALLOCATE( pco2_balinc(jpi,jpj,jpk,jptra) ) 461 pco2_balinc(:,:,:,:) = 0.0 462 #endif 463 ENDIF 464 IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) & 465 & .OR.( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) ) THEN 353 466 354 467 !-------------------------------------------------------------------- … … 424 537 ENDIF 425 538 539 IF ( ln_logchlinc ) THEN 540 CALL iom_get( inum, jpdom_autoglo, 'bckinlogchl', logchl_bkginc(:,:), 1 ) 541 ! Apply the masks 542 logchl_bkginc(:,:) = logchl_bkginc(:,:) * tmask(:,:,1) 543 ! Set missing increments to 0.0 rather than 1e+20 544 ! to allow for differences in masks 545 WHERE( ABS( logchl_bkginc(:,:) ) > 1.0e+10 ) logchl_bkginc(:,:) = 0.0 546 ENDIF 547 548 IF ( ln_pco2inc .OR. ln_fco2inc ) THEN 549 IF ( ln_pco2inc ) THEN 550 CALL iom_get( inum, jpdom_autoglo, 'bckinpco2', pco2_bkginc(:,:), 1 ) 551 ELSE IF ( ln_fco2inc ) THEN 552 CALL iom_get( inum, jpdom_autoglo, 'bckinfco2', pco2_bkginc(:,:), 1 ) 553 ENDIF 554 ! Apply the masks 555 pco2_bkginc(:,:) = pco2_bkginc(:,:) * tmask(:,:,1) 556 ! Set missing increments to 0.0 rather than 1e+20 557 ! to allow for differences in masks 558 WHERE( ABS( pco2_bkginc(:,:) ) > 1.0e+10 ) pco2_bkginc(:,:) = 0.0 559 ENDIF 560 426 561 CALL iom_close( inum ) 427 562 … … 535 670 536 671 ENDIF 672 673 #if defined key_hadocc || (defined key_medusa && defined key_foam_medusa) 674 IF ( ln_logchlinc ) THEN 675 676 ALLOCATE( pgrow_avg_bkg(jpi,jpj) ) 677 ALLOCATE( ploss_avg_bkg(jpi,jpj) ) 678 ALLOCATE( phyt_avg_bkg(jpi,jpj) ) 679 ALLOCATE( mld_max_bkg(jpi,jpj) ) 680 ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) ) 681 pgrow_avg_bkg(:,:) = 0.0 682 ploss_avg_bkg(:,:) = 0.0 683 phyt_avg_bkg(:,:) = 0.0 684 mld_max_bkg(:,:) = 0.0 685 tracer_bkg(:,:,:,:) = 0.0 686 687 #if defined key_hadocc 688 ALLOCATE( chl_bkg(jpi,jpj) ) 689 ALLOCATE( cchl_p_bkg(jpi,jpj) ) 690 chl_bkg(:,:) = 0.0 691 cchl_p_bkg(:,:) = 0.0 692 #endif 693 694 !-------------------------------------------------------------------- 695 ! Read background variables for logchl assimilation 696 ! Some only required if performing balancing 697 !-------------------------------------------------------------------- 698 699 CALL iom_open( c_asmbkg, inum ) 700 701 #if defined key_hadocc 702 CALL iom_get( inum, jpdom_autoglo, 'hadocc_chl', chl_bkg ) 703 CALL iom_get( inum, jpdom_autoglo, 'hadocc_cchl', cchl_p_bkg ) 704 chl_bkg(:,:) = chl_bkg(:,:) * tmask(:,:,1) 705 cchl_p_bkg(:,:) = cchl_p_bkg(:,:) * tmask(:,:,1) 706 #elif defined key_medusa 707 CALL iom_get( inum, jpdom_autoglo, 'medusa_chn', tracer_bkg(:,:,:,jpchn) ) 708 CALL iom_get( inum, jpdom_autoglo, 'medusa_chd', tracer_bkg(:,:,:,jpchd) ) 709 #endif 710 711 IF ( ln_logchlbal ) THEN 712 713 CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg', pgrow_avg_bkg ) 714 CALL iom_get( inum, jpdom_autoglo, 'ploss_avg', ploss_avg_bkg ) 715 CALL iom_get( inum, jpdom_autoglo, 'phyt_avg', phyt_avg_bkg ) 716 CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg ) 717 pgrow_avg_bkg(:,:) = pgrow_avg_bkg(:,:) * tmask(:,:,1) 718 ploss_avg_bkg(:,:) = ploss_avg_bkg(:,:) * tmask(:,:,1) 719 phyt_avg_bkg(:,:) = phyt_avg_bkg(:,:) * tmask(:,:,1) 720 mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1) 721 722 #if defined key_hadocc 723 CALL iom_get( inum, jpdom_autoglo, 'hadocc_nut', tracer_bkg(:,:,:,jp_had_nut) ) 724 CALL iom_get( inum, jpdom_autoglo, 'hadocc_phy', tracer_bkg(:,:,:,jp_had_phy) ) 725 CALL iom_get( inum, jpdom_autoglo, 'hadocc_zoo', tracer_bkg(:,:,:,jp_had_zoo) ) 726 CALL iom_get( inum, jpdom_autoglo, 'hadocc_pdn', tracer_bkg(:,:,:,jp_had_pdn) ) 727 CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) ) 728 CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) ) 729 #elif defined key_medusa 730 CALL iom_get( inum, jpdom_autoglo, 'medusa_phn', tracer_bkg(:,:,:,jpphn) ) 731 CALL iom_get( inum, jpdom_autoglo, 'medusa_phd', tracer_bkg(:,:,:,jpphd) ) 732 CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) ) 733 CALL iom_get( inum, jpdom_autoglo, 'medusa_zmi', tracer_bkg(:,:,:,jpzmi) ) 734 CALL iom_get( inum, jpdom_autoglo, 'medusa_zme', tracer_bkg(:,:,:,jpzme) ) 735 CALL iom_get( inum, jpdom_autoglo, 'medusa_din', tracer_bkg(:,:,:,jpdin) ) 736 CALL iom_get( inum, jpdom_autoglo, 'medusa_sil', tracer_bkg(:,:,:,jpsil) ) 737 CALL iom_get( inum, jpdom_autoglo, 'medusa_fer', tracer_bkg(:,:,:,jpfer) ) 738 CALL iom_get( inum, jpdom_autoglo, 'medusa_det', tracer_bkg(:,:,:,jpdet) ) 739 CALL iom_get( inum, jpdom_autoglo, 'medusa_dtc', tracer_bkg(:,:,:,jpdtc) ) 740 CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) ) 741 CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) ) 742 CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) ) 743 #endif 744 ELSE IF ( ln_pco2inc .OR. ln_fco2inc ) THEN 745 #if defined key_hadocc 746 CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) ) 747 CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) ) 748 #elif defined key_medusa 749 CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) ) 750 CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) ) 751 #endif 752 ENDIF 753 754 CALL iom_close( inum ) 755 756 DO jt = 1, jptra 757 tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:) 758 END DO 759 760 ELSE IF ( ln_pco2inc .OR. ln_fco2inc ) THEN 761 762 ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) ) 763 tracer_bkg(:,:,:,:) = 0.0 764 765 CALL iom_open( c_asmbkg, inum ) 766 767 #if defined key_hadocc 768 CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) ) 769 CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) ) 770 #elif defined key_medusa 771 CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) ) 772 CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) ) 773 #endif 774 775 CALL iom_close( inum ) 776 777 DO jt = 1, jptra 778 tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:) 779 END DO 780 781 ENDIF 782 #endif 537 783 ! 538 784 END SUBROUTINE asm_inc_init … … 1155 1401 1156 1402 END SUBROUTINE seaice_asm_inc 1403 1404 SUBROUTINE logchl_asm_inc( kt ) 1405 !!---------------------------------------------------------------------- 1406 !! *** ROUTINE logchl_asm_inc *** 1407 !! 1408 !! ** Purpose : Apply the chlorophyll assimilation increments. 1409 !! 1410 !! ** Method : Calculate increments to state variables using nitrogen 1411 !! balancing. 1412 !! Direct initialization or Incremental Analysis Updating. 1413 !! 1414 !! ** Action : 1415 !!---------------------------------------------------------------------- 1416 INTEGER, INTENT(IN) :: kt ! Current time step 1417 ! 1418 INTEGER :: jk ! Loop counter 1419 INTEGER :: it ! Index 1420 REAL(wp) :: zincwgt ! IAU weight for current time step 1421 REAL(wp) :: zincper ! IAU interval in seconds 1422 !!---------------------------------------------------------------------- 1423 1424 IF ( kt <= nit000 ) THEN 1425 1426 zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt 1427 1428 #if defined key_medusa && defined key_foam_medusa 1429 CALL asm_logchl_bal_medusa( logchl_bkginc, zincper, mld_choice_bgc, & 1430 & rn_maxchlinc, ln_logchlbal, & 1431 & pgrow_avg_bkg, ploss_avg_bkg, & 1432 & phyt_avg_bkg, mld_max_bkg, & 1433 & tracer_bkg, logchl_balinc ) 1434 #elif defined key_hadocc 1435 CALL asm_logchl_bal_hadocc( logchl_bkginc, zincper, mld_choice_bgc, & 1436 & rn_maxchlinc, ln_logchlbal, & 1437 & pgrow_avg_bkg, ploss_avg_bkg, & 1438 & phyt_avg_bkg, mld_max_bkg, & 1439 & chl_bkg, cchl_p_bkg, & 1440 & tracer_bkg, logchl_balinc ) 1441 #else 1442 CALL ctl_stop( 'Attempting to assimilate logchl, ', & 1443 & 'but not defined a biogeochemical model' ) 1444 #endif 1445 1446 ENDIF 1447 1448 IF ( ln_asmiau ) THEN 1449 1450 !-------------------------------------------------------------------- 1451 ! Incremental Analysis Updating 1452 !-------------------------------------------------------------------- 1453 1454 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 1455 1456 it = kt - nit000 + 1 1457 zincwgt = wgtiau(it) ! IAU weight for the current time step 1458 ! note this is not a tendency so should not be divided by rdt 1459 1460 IF(lwp) THEN 1461 WRITE(numout,*) 1462 WRITE(numout,*) 'logchl_asm_inc : logchl IAU at time step = ', & 1463 & kt,' with IAU weight = ', wgtiau(it) 1464 WRITE(numout,*) '~~~~~~~~~~~~' 1465 ENDIF 1466 1467 ! Update the biogeochemical variables 1468 ! Add directly to trn and trb, rather than to tra, as not a tendency 1469 #if defined key_medusa && defined key_foam_medusa 1470 DO jk = 1, jpkm1 1471 trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + & 1472 & logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 1473 trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + & 1474 & logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 1475 END DO 1476 #elif defined key_hadocc 1477 DO jk = 1, jpkm1 1478 trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + & 1479 & logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 1480 trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + & 1481 & logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 1482 END DO 1483 #endif 1484 1485 ! Do not deallocate arrays - needed by asm_bal_wri 1486 ! which is called at end of model run 1487 ENDIF 1488 1489 ELSEIF ( ln_asmdin ) THEN 1490 1491 !-------------------------------------------------------------------- 1492 ! Direct Initialization 1493 !-------------------------------------------------------------------- 1494 1495 IF ( kt == nitdin_r ) THEN 1496 1497 neuler = 0 ! Force Euler forward step 1498 1499 #if defined key_medusa && defined key_foam_medusa 1500 ! Initialize the now fields with the background + increment 1501 ! Background currently is what the model is initialised with 1502 CALL ctl_warn( ' Doing direct initialisation of MEDUSA with chlorophyll assimilation', & 1503 & ' Background state is taken from model rather than background file' ) 1504 trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & 1505 & logchl_balinc(:,:,:,jp_msa0:jp_msa1) 1506 trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) 1507 #elif defined key_hadocc 1508 ! Initialize the now fields with the background + increment 1509 ! Background currently is what the model is initialised with 1510 CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', & 1511 & ' Background state is taken from model rather than background file' ) 1512 trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & 1513 & logchl_balinc(:,:,:,jp_had0:jp_had1) 1514 trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) 1515 #endif 1516 1517 ! Do not deallocate arrays - needed by asm_bal_wri 1518 ! which is called at end of model run 1519 ENDIF 1520 ! 1521 ENDIF 1522 ! 1523 END SUBROUTINE logchl_asm_inc 1524 1525 1526 SUBROUTINE pco2_asm_inc( kt ) 1527 !!---------------------------------------------------------------------- 1528 !! *** ROUTINE pco2_asm_inc *** 1529 !! 1530 !! ** Purpose : Apply the pco2/fco2 assimilation increments. 1531 !! 1532 !! ** Method : Calculate increments to state variables using carbon 1533 !! balancing. 1534 !! Direct initialization or Incremental Analysis Updating. 1535 !! 1536 !! ** Action : 1537 !!---------------------------------------------------------------------- 1538 INTEGER, INTENT(IN) :: kt ! Current time step 1539 ! 1540 INTEGER :: ji, jj, jk ! Loop counters 1541 INTEGER :: it ! Index 1542 INTEGER :: jkmax ! Index of mixed layer depth 1543 REAL(wp) :: zincwgt ! IAU weight for current time step 1544 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! Mixed layer depth for balancing 1545 REAL(wp), DIMENSION(jpi,jpj) :: pco2_bkginc_temp ! For fCO2 conversion 1546 REAL(wp), DIMENSION(jpi,jpj) :: dic_bkg_temp ! DIC background 1547 REAL(wp), DIMENSION(jpi,jpj) :: alk_bkg_temp ! Alkalinity background 1548 REAL(wp), DIMENSION(jpi,jpj) :: tem_bkg_temp ! Temperature background 1549 REAL(wp), DIMENSION(jpi,jpj) :: sal_bkg_temp ! Salinity background 1550 REAL(wp), DIMENSION(1) :: patm ! Atmospheric pressure 1551 REAL(wp), DIMENSION(1) :: phyd ! Hydrostatic pressure 1552 1553 ! Coefficients for fCO2 to pCO2 conversion 1554 REAL(wp) :: zcoef_fco2_1 = -1636.75 1555 REAL(wp) :: zcoef_fco2_2 = 12.0408 1556 REAL(wp) :: zcoef_fco2_3 = 0.0327957 1557 REAL(wp) :: zcoef_fco2_4 = 0.0000316528 1558 REAL(wp) :: zcoef_fco2_5 = 2.0 1559 REAL(wp) :: zcoef_fco2_6 = 57.7 1560 REAL(wp) :: zcoef_fco2_7 = 0.118 1561 REAL(wp) :: zcoef_fco2_8 = 82.0578 1562 !!---------------------------------------------------------------------- 1563 1564 IF ( kt <= nit000 ) THEN 1565 1566 !-------------------------------------------------------------------- 1567 ! DIC and alkalinity balancing 1568 !-------------------------------------------------------------------- 1569 1570 IF ( ln_fco2inc ) THEN 1571 #if defined key_medusa && defined key_roam 1572 ! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine 1573 patm(1) = 1.0 1574 phyd(1) = 0.0 1575 DO jj = 1, jpj 1576 DO ji = 1, jpi 1577 CALL f2pCO2( pco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) ) 1578 END DO 1579 END DO 1580 #else 1581 ! If assimilating fCO2, then convert to pCO2 using temperature 1582 ! See flux_gas.F90 within HadOCC for details of calculation 1583 pco2_bkginc_temp(:,:) = pco2_bkginc(:,:) / & 1584 & EXP((zcoef_fco2_1 + & 1585 & zcoef_fco2_2 * (tsn(:,:,1,1)+rt0) - & 1586 & zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + & 1587 & zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + & 1588 & zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0))) / & 1589 & (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0))) 1590 #endif 1591 ELSE 1592 pco2_bkginc_temp(:,:) = pco2_bkginc(:,:) 1593 ENDIF 1594 1595 ! Account for physics assimilation if required 1596 IF ( ln_trainc ) THEN 1597 tem_bkg_temp(:,:) = tsn(:,:,1,1) + t_bkginc(:,:,1) 1598 sal_bkg_temp(:,:) = tsn(:,:,1,2) + s_bkginc(:,:,1) 1599 ELSE 1600 tem_bkg_temp(:,:) = tsn(:,:,1,1) 1601 sal_bkg_temp(:,:) = tsn(:,:,1,2) 1602 ENDIF 1603 1604 #if defined key_medusa 1605 ! Account for logchl balancing if required 1606 IF ( ln_logchlinc .AND. ln_logchlbal ) THEN 1607 dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + logchl_balinc(:,:,1,jpdic) 1608 alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + logchl_balinc(:,:,1,jpalk) 1609 ELSE 1610 dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) 1611 alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) 1612 ENDIF 1613 1614 CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), & 1615 & tem_bkg_temp(:,:), sal_bkg_temp(:,:), & 1616 & pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) ) 1617 1618 #elif defined key_hadocc 1619 ! Account for logchl balancing if required 1620 IF ( ln_logchlinc .AND. ln_logchlbal ) THEN 1621 dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + logchl_balinc(:,:,1,jp_had_dic) 1622 alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + logchl_balinc(:,:,1,jp_had_alk) 1623 ELSE 1624 dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) 1625 alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) 1626 ENDIF 1627 1628 CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), & 1629 & tem_bkg_temp(:,:), sal_bkg_temp(:,:), & 1630 & pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) ) 1631 1632 #else 1633 CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', & 1634 & 'but not defined a biogeochemical model' ) 1635 #endif 1636 1637 ! Select mixed layer 1638 SELECT CASE( mld_choice_bgc ) 1639 CASE ( 1 ) ! Turbocline/mixing depth [W points] 1640 zmld(:,:) = hmld(:,:) 1641 CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points] 1642 zmld(:,:) = hmlp(:,:) 1643 CASE ( 3 ) ! Kara MLD [Interpolated] 1644 #if defined key_karaml 1645 IF ( ln_kara ) THEN 1646 zmld(:,:) = hmld_kara(:,:) 1647 ELSE 1648 CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', & 1649 & ' but ln_kara=.false.' ) 1650 ENDIF 1651 #else 1652 CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', & 1653 & ' but is not defined' ) 1654 #endif 1655 CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points] 1656 !zmld(:,:) = hmld_tref(:,:) 1657 CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', & 1658 & ' but is not available in this version' ) 1659 CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points] 1660 zmld(:,:) = hmlpt(:,:) 1661 END SELECT 1662 1663 ! Propagate through mixed layer 1664 DO jj = 1, jpj 1665 DO ji = 1, jpi 1666 ! 1667 jkmax = jpk-1 1668 DO jk = jpk-1, 1, -1 1669 IF ( ( zmld(ji,jj) > gdepw_n(ji,jj,jk) ) .AND. & 1670 & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN 1671 zmld(ji,jj) = gdepw_n(ji,jj,jk+1) 1672 jkmax = jk 1673 ENDIF 1674 END DO 1675 ! 1676 DO jk = 2, jkmax 1677 pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:) 1678 END DO 1679 ! 1680 END DO 1681 END DO 1682 1683 ENDIF 1684 1685 IF ( ln_asmiau ) THEN 1686 1687 !-------------------------------------------------------------------- 1688 ! Incremental Analysis Updating 1689 !-------------------------------------------------------------------- 1690 1691 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 1692 1693 it = kt - nit000 + 1 1694 zincwgt = wgtiau(it) ! IAU weight for the current time step 1695 ! note this is not a tendency so should not be divided by rdt 1696 1697 IF(lwp) THEN 1698 WRITE(numout,*) 1699 IF ( ln_pco2inc ) THEN 1700 WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', & 1701 & kt,' with IAU weight = ', wgtiau(it) 1702 ELSE IF ( ln_fco2inc ) THEN 1703 WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', & 1704 & kt,' with IAU weight = ', wgtiau(it) 1705 ENDIF 1706 WRITE(numout,*) '~~~~~~~~~~~~' 1707 ENDIF 1708 1709 ! Update the biogeochemical variables 1710 ! Add directly to trn and trb, rather than to tra, as not a tendency 1711 #if defined key_medusa && defined key_foam_medusa 1712 DO jk = 1, jpkm1 1713 trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + & 1714 & pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 1715 trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + & 1716 & pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 1717 END DO 1718 #elif defined key_hadocc 1719 DO jk = 1, jpkm1 1720 trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + & 1721 & pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 1722 trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + & 1723 & pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 1724 END DO 1725 #endif 1726 1727 ! Do not deallocate arrays - needed by asm_bal_wri 1728 ! which is called at end of model run 1729 1730 ENDIF 1731 1732 ELSEIF ( ln_asmdin ) THEN 1733 1734 !-------------------------------------------------------------------- 1735 ! Direct Initialization 1736 !-------------------------------------------------------------------- 1737 1738 IF ( kt == nitdin_r ) THEN 1739 1740 neuler = 0 ! Force Euler forward step 1741 1742 #if defined key_medusa && defined key_foam_medusa 1743 ! Initialize the now fields with the background + increment 1744 ! Background currently is what the model is initialised with 1745 CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pCO2 assimilation', & 1746 & ' Background state is taken from model rather than background file' ) 1747 trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & 1748 & pco2_balinc(:,:,:,jp_msa0:jp_msa1) 1749 trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) 1750 #elif defined key_hadocc 1751 ! Initialize the now fields with the background + increment 1752 ! Background currently is what the model is initialised with 1753 CALL ctl_warn( ' Doing direct initialisation of HadOCC with pCO2 assimilation', & 1754 & ' Background state is taken from model rather than background file' ) 1755 trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & 1756 & pco2_balinc(:,:,:,jp_had0:jp_had1) 1757 trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) 1758 #endif 1759 1760 ! Do not deallocate arrays - needed by asm_bal_wri 1761 ! which is called at end of model run 1762 ENDIF 1763 ! 1764 ENDIF 1765 ! 1766 END SUBROUTINE pco2_asm_inc 1157 1767 1158 1768 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
