Changeset 9292 for branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
- Timestamp:
- 2018-01-30T19:41:58+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r8400 r9292 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 CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg ) 753 mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1) 754 ENDIF 755 756 CALL iom_close( inum ) 757 758 DO jt = 1, jptra 759 tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:) 760 END DO 761 762 ELSE IF ( ln_pco2inc .OR. ln_fco2inc ) THEN 763 764 ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) ) 765 ALLOCATE( mld_max_bkg(jpi,jpj) ) 766 tracer_bkg(:,:,:,:) = 0.0 767 mld_max_bkg(:,:) = 0.0 768 769 CALL iom_open( c_asmbkg, inum ) 770 771 #if defined key_hadocc 772 CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) ) 773 CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) ) 774 #elif defined key_medusa 775 CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) ) 776 CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) ) 777 #endif 778 CALL iom_get( inum, jpdom_autoglo, 'mld_max', mld_max_bkg ) 779 780 CALL iom_close( inum ) 781 782 DO jt = 1, jptra 783 tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:) 784 END DO 785 mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1) 786 787 ENDIF 788 #endif 537 789 ! 538 790 END SUBROUTINE asm_inc_init … … 1155 1407 1156 1408 END SUBROUTINE seaice_asm_inc 1409 1410 SUBROUTINE logchl_asm_inc( kt ) 1411 !!---------------------------------------------------------------------- 1412 !! *** ROUTINE logchl_asm_inc *** 1413 !! 1414 !! ** Purpose : Apply the chlorophyll assimilation increments. 1415 !! 1416 !! ** Method : Calculate increments to state variables using nitrogen 1417 !! balancing. 1418 !! Direct initialization or Incremental Analysis Updating. 1419 !! 1420 !! ** Action : 1421 !!---------------------------------------------------------------------- 1422 INTEGER, INTENT(IN) :: kt ! Current time step 1423 ! 1424 INTEGER :: jk ! Loop counter 1425 INTEGER :: it ! Index 1426 REAL(wp) :: zincwgt ! IAU weight for current time step 1427 REAL(wp) :: zincper ! IAU interval in seconds 1428 !!---------------------------------------------------------------------- 1429 1430 IF ( kt <= nit000 ) THEN 1431 1432 zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt 1433 1434 #if defined key_medusa && defined key_foam_medusa 1435 CALL asm_logchl_bal_medusa( logchl_bkginc, zincper, mld_choice_bgc, & 1436 & rn_maxchlinc, ln_logchlbal, ln_asmdin, & 1437 & pgrow_avg_bkg, ploss_avg_bkg, & 1438 & phyt_avg_bkg, mld_max_bkg, & 1439 & tracer_bkg, logchl_balinc ) 1440 #elif defined key_hadocc 1441 CALL asm_logchl_bal_hadocc( logchl_bkginc, zincper, mld_choice_bgc, & 1442 & rn_maxchlinc, ln_logchlbal, ln_asmdin, & 1443 & pgrow_avg_bkg, ploss_avg_bkg, & 1444 & phyt_avg_bkg, mld_max_bkg, & 1445 & chl_bkg, cchl_p_bkg, & 1446 & tracer_bkg, logchl_balinc ) 1447 #else 1448 CALL ctl_stop( 'Attempting to assimilate logchl, ', & 1449 & 'but not defined a biogeochemical model' ) 1450 #endif 1451 1452 ENDIF 1453 1454 IF ( ln_asmiau ) THEN 1455 1456 !-------------------------------------------------------------------- 1457 ! Incremental Analysis Updating 1458 !-------------------------------------------------------------------- 1459 1460 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 1461 1462 it = kt - nit000 + 1 1463 zincwgt = wgtiau(it) ! IAU weight for the current time step 1464 ! note this is not a tendency so should not be divided by rdt 1465 1466 IF(lwp) THEN 1467 WRITE(numout,*) 1468 WRITE(numout,*) 'logchl_asm_inc : logchl IAU at time step = ', & 1469 & kt,' with IAU weight = ', wgtiau(it) 1470 WRITE(numout,*) '~~~~~~~~~~~~' 1471 ENDIF 1472 1473 ! Update the biogeochemical variables 1474 ! Add directly to trn and trb, rather than to tra, because tra gets 1475 ! reset to zero at the start of trc_stp, called after this routine 1476 #if defined key_medusa && defined key_foam_medusa 1477 DO jk = 1, jpkm1 1478 trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + & 1479 & logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 1480 trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + & 1481 & logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 1482 END DO 1483 #elif defined key_hadocc 1484 DO jk = 1, jpkm1 1485 trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + & 1486 & logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 1487 trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + & 1488 & logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 1489 END DO 1490 #endif 1491 1492 ! Do not deallocate arrays - needed by asm_bal_wri 1493 ! which is called at end of model run 1494 ENDIF 1495 1496 ELSEIF ( ln_asmdin ) THEN 1497 1498 !-------------------------------------------------------------------- 1499 ! Direct Initialization 1500 !-------------------------------------------------------------------- 1501 1502 IF ( kt == nitdin_r ) THEN 1503 1504 neuler = 0 ! Force Euler forward step 1505 1506 #if defined key_medusa && defined key_foam_medusa 1507 ! Initialize the now fields with the background + increment 1508 ! Background currently is what the model is initialised with 1509 CALL ctl_warn( ' Doing direct initialisation of MEDUSA with chlorophyll assimilation', & 1510 & ' Background state is taken from model rather than background file' ) 1511 trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & 1512 & logchl_balinc(:,:,:,jp_msa0:jp_msa1) 1513 trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) 1514 #elif defined key_hadocc 1515 ! Initialize the now fields with the background + increment 1516 ! Background currently is what the model is initialised with 1517 CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', & 1518 & ' Background state is taken from model rather than background file' ) 1519 trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & 1520 & logchl_balinc(:,:,:,jp_had0:jp_had1) 1521 trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) 1522 #endif 1523 1524 ! Do not deallocate arrays - needed by asm_bal_wri 1525 ! which is called at end of model run 1526 ENDIF 1527 ! 1528 ENDIF 1529 ! 1530 END SUBROUTINE logchl_asm_inc 1531 1532 1533 SUBROUTINE pco2_asm_inc( kt ) 1534 !!---------------------------------------------------------------------- 1535 !! *** ROUTINE pco2_asm_inc *** 1536 !! 1537 !! ** Purpose : Apply the pco2/fco2 assimilation increments. 1538 !! 1539 !! ** Method : Calculate increments to state variables using carbon 1540 !! balancing. 1541 !! Direct initialization or Incremental Analysis Updating. 1542 !! 1543 !! ** Action : 1544 !!---------------------------------------------------------------------- 1545 INTEGER, INTENT(IN) :: kt ! Current time step 1546 ! 1547 INTEGER :: ji, jj, jk ! Loop counters 1548 INTEGER :: it ! Index 1549 INTEGER :: jkmax ! Index of mixed layer depth 1550 REAL(wp) :: zincwgt ! IAU weight for current time step 1551 REAL(wp), DIMENSION(jpi,jpj) :: zmld ! Mixed layer depth for balancing 1552 REAL(wp), DIMENSION(jpi,jpj) :: pco2_bkginc_temp ! For fCO2 conversion 1553 REAL(wp), DIMENSION(jpi,jpj) :: dic_bkg_temp ! DIC background 1554 REAL(wp), DIMENSION(jpi,jpj) :: alk_bkg_temp ! Alkalinity background 1555 REAL(wp), DIMENSION(jpi,jpj) :: tem_bkg_temp ! Temperature background 1556 REAL(wp), DIMENSION(jpi,jpj) :: sal_bkg_temp ! Salinity background 1557 REAL(wp), DIMENSION(1) :: patm ! Atmospheric pressure 1558 REAL(wp), DIMENSION(1) :: phyd ! Hydrostatic pressure 1559 1560 ! Coefficients for fCO2 to pCO2 conversion 1561 REAL(wp) :: zcoef_fco2_1 = -1636.75 1562 REAL(wp) :: zcoef_fco2_2 = 12.0408 1563 REAL(wp) :: zcoef_fco2_3 = 0.0327957 1564 REAL(wp) :: zcoef_fco2_4 = 0.0000316528 1565 REAL(wp) :: zcoef_fco2_5 = 2.0 1566 REAL(wp) :: zcoef_fco2_6 = 57.7 1567 REAL(wp) :: zcoef_fco2_7 = 0.118 1568 REAL(wp) :: zcoef_fco2_8 = 82.0578 1569 !!---------------------------------------------------------------------- 1570 1571 IF ( kt <= nit000 ) THEN 1572 1573 !-------------------------------------------------------------------- 1574 ! DIC and alkalinity balancing 1575 !-------------------------------------------------------------------- 1576 1577 IF ( ln_fco2inc ) THEN 1578 #if defined key_medusa && defined key_roam 1579 ! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine 1580 patm(1) = 1.0 1581 phyd(1) = 0.0 1582 DO jj = 1, jpj 1583 DO ji = 1, jpi 1584 CALL f2pCO2( pco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) ) 1585 END DO 1586 END DO 1587 #else 1588 ! If assimilating fCO2, then convert to pCO2 using temperature 1589 ! See flux_gas.F90 within HadOCC for details of calculation 1590 pco2_bkginc_temp(:,:) = pco2_bkginc(:,:) / & 1591 & EXP((zcoef_fco2_1 + & 1592 & zcoef_fco2_2 * (tsn(:,:,1,1)+rt0) - & 1593 & zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + & 1594 & zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + & 1595 & zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0))) / & 1596 & (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0))) 1597 #endif 1598 ELSE 1599 pco2_bkginc_temp(:,:) = pco2_bkginc(:,:) 1600 ENDIF 1601 1602 ! Account for physics assimilation if required 1603 IF ( ln_trainc ) THEN 1604 tem_bkg_temp(:,:) = tsn(:,:,1,1) + t_bkginc(:,:,1) 1605 sal_bkg_temp(:,:) = tsn(:,:,1,2) + s_bkginc(:,:,1) 1606 ELSE 1607 tem_bkg_temp(:,:) = tsn(:,:,1,1) 1608 sal_bkg_temp(:,:) = tsn(:,:,1,2) 1609 ENDIF 1610 1611 #if defined key_medusa 1612 ! Account for logchl balancing if required 1613 IF ( ln_logchlinc .AND. ln_logchlbal ) THEN 1614 dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + logchl_balinc(:,:,1,jpdic) 1615 alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + logchl_balinc(:,:,1,jpalk) 1616 ELSE 1617 dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) 1618 alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) 1619 ENDIF 1620 1621 CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), & 1622 & tem_bkg_temp(:,:), sal_bkg_temp(:,:), & 1623 & pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) ) 1624 1625 #elif defined key_hadocc 1626 ! Account for logchl balancing if required 1627 IF ( ln_logchlinc .AND. ln_logchlbal ) THEN 1628 dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + logchl_balinc(:,:,1,jp_had_dic) 1629 alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + logchl_balinc(:,:,1,jp_had_alk) 1630 ELSE 1631 dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) 1632 alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) 1633 ENDIF 1634 1635 CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), & 1636 & tem_bkg_temp(:,:), sal_bkg_temp(:,:), & 1637 & pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) ) 1638 1639 #else 1640 CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', & 1641 & 'but not defined a biogeochemical model' ) 1642 #endif 1643 1644 ! Select mixed layer 1645 IF ( ln_asmdin ) THEN 1646 #if defined key_hadocc || (defined key_medusa && defined key_foam_medusa) 1647 CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', & 1648 & ' Mixed layer depth taken to be background maximum mld_max_bkg' ) 1649 zmld(:,:) = mld_max_bkg(:,:) 1650 #else 1651 CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' ) 1652 #endif 1653 ELSE 1654 SELECT CASE( mld_choice_bgc ) 1655 CASE ( 1 ) ! Turbocline/mixing depth [W points] 1656 zmld(:,:) = hmld(:,:) 1657 CASE ( 2 ) ! Density criterion (0.01 kg/m^3 change from 10m) [W points] 1658 zmld(:,:) = hmlp(:,:) 1659 CASE ( 3 ) ! Kara MLD [Interpolated] 1660 #if defined key_karaml 1661 IF ( ln_kara ) THEN 1662 zmld(:,:) = hmld_kara(:,:) 1663 ELSE 1664 CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', & 1665 & ' but ln_kara=.false.' ) 1666 ENDIF 1667 #else 1668 CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', & 1669 & ' but is not defined' ) 1670 #endif 1671 CASE ( 4 ) ! Temperature criterion (0.2 K change from surface) [T points] 1672 !zmld(:,:) = hmld_tref(:,:) 1673 CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', & 1674 & ' but is not available in this version' ) 1675 CASE ( 5 ) ! Density criterion (0.01 kg/m^3 change from 10m) [T points] 1676 zmld(:,:) = hmlpt(:,:) 1677 END SELECT 1678 ENDIF 1679 1680 ! Propagate through mixed layer 1681 DO jj = 1, jpj 1682 DO ji = 1, jpi 1683 ! 1684 jkmax = jpk-1 1685 DO jk = jpk-1, 1, -1 1686 IF ( ( zmld(ji,jj) > gdepw_n(ji,jj,jk) ) .AND. & 1687 & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN 1688 zmld(ji,jj) = gdepw_n(ji,jj,jk+1) 1689 jkmax = jk 1690 ENDIF 1691 END DO 1692 ! 1693 #if defined key_top 1694 DO jk = 2, jkmax 1695 pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:) 1696 END DO 1697 #endif 1698 ! 1699 END DO 1700 END DO 1701 1702 ENDIF 1703 1704 IF ( ln_asmiau ) THEN 1705 1706 !-------------------------------------------------------------------- 1707 ! Incremental Analysis Updating 1708 !-------------------------------------------------------------------- 1709 1710 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 1711 1712 it = kt - nit000 + 1 1713 zincwgt = wgtiau(it) ! IAU weight for the current time step 1714 ! note this is not a tendency so should not be divided by rdt 1715 1716 IF(lwp) THEN 1717 WRITE(numout,*) 1718 IF ( ln_pco2inc ) THEN 1719 WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', & 1720 & kt,' with IAU weight = ', wgtiau(it) 1721 ELSE IF ( ln_fco2inc ) THEN 1722 WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', & 1723 & kt,' with IAU weight = ', wgtiau(it) 1724 ENDIF 1725 WRITE(numout,*) '~~~~~~~~~~~~' 1726 ENDIF 1727 1728 ! Update the biogeochemical variables 1729 ! Add directly to trn and trb, rather than to tra, because tra gets 1730 ! reset to zero at the start of trc_stp, called after this routine 1731 #if defined key_medusa && defined key_foam_medusa 1732 DO jk = 1, jpkm1 1733 trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + & 1734 & pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 1735 trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + & 1736 & pco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 1737 END DO 1738 #elif defined key_hadocc 1739 DO jk = 1, jpkm1 1740 trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + & 1741 & pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 1742 trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + & 1743 & pco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 1744 END DO 1745 #endif 1746 1747 ! Do not deallocate arrays - needed by asm_bal_wri 1748 ! which is called at end of model run 1749 1750 ENDIF 1751 1752 ELSEIF ( ln_asmdin ) THEN 1753 1754 !-------------------------------------------------------------------- 1755 ! Direct Initialization 1756 !-------------------------------------------------------------------- 1757 1758 IF ( kt == nitdin_r ) THEN 1759 1760 neuler = 0 ! Force Euler forward step 1761 1762 #if defined key_medusa && defined key_foam_medusa 1763 ! Initialize the now fields with the background + increment 1764 ! Background currently is what the model is initialised with 1765 CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pCO2 assimilation', & 1766 & ' Background state is taken from model rather than background file' ) 1767 trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & 1768 & pco2_balinc(:,:,:,jp_msa0:jp_msa1) 1769 trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) 1770 #elif defined key_hadocc 1771 ! Initialize the now fields with the background + increment 1772 ! Background currently is what the model is initialised with 1773 CALL ctl_warn( ' Doing direct initialisation of HadOCC with pCO2 assimilation', & 1774 & ' Background state is taken from model rather than background file' ) 1775 trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & 1776 & pco2_balinc(:,:,:,jp_had0:jp_had1) 1777 trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) 1778 #endif 1779 1780 ! Do not deallocate arrays - needed by asm_bal_wri 1781 ! which is called at end of model run 1782 ENDIF 1783 ! 1784 ENDIF 1785 ! 1786 END SUBROUTINE pco2_asm_inc 1157 1787 1158 1788 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
