Changeset 8456 for branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
 Timestamp:
 20170823T14:20:41+02:00 (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
r8440 r8456 23 23 !! seaice_asm_inc : Apply the seaice increment 24 24 !! logchl_asm_inc : Apply the logchl increment 25 !! pco2_asm_inc : Apply the pCO2/fCO2 increment 25 26 !! 26 27 USE wrk_nemo ! Memory Allocation … … 64 65 USE asmlogchlbal_medusa, ONLY: & 65 66 & asm_logchl_bal_medusa 67 USE asmpco2bal, ONLY: & 68 & asm_pco2_bal 66 69 USE par_medusa 70 USE mocsy_mainmod 67 71 #elif defined key_hadocc 68 72 USE asmlogchlbal_hadocc, ONLY: & 69 73 & asm_logchl_bal_hadocc 74 USE asmpco2bal, ONLY: & 75 & asm_pco2_bal 70 76 USE par_hadocc 71 77 #endif … … 81 87 PUBLIC seaice_asm_inc !: Apply the seaice increment 82 88 PUBLIC logchl_asm_inc !: Apply the logchl increment 89 PUBLIC pco2_asm_inc !: Apply the pCO2/fCO2 increment 83 90 84 91 #if defined key_asminc … … 97 104 LOGICAL, PUBLIC :: ln_logchlinc = .FALSE. !: No log10(chlorophyll) increment 98 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 99 108 LOGICAL, PUBLIC :: ln_salfix = .FALSE. !: Apply minimum salinity check 100 109 LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing … … 123 132 124 133 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: logchl_bkginc !: Increment to background logchl 134 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: pco2_bkginc !: Increment to background pCO2 125 135 #if defined key_top 126 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 127 138 #endif 128 139 #if defined key_hadocc  (defined key_medusa && defined key_foam_medusa) … … 192 203 & ln_trainc, ln_dyninc, ln_sshinc, & 193 204 & ln_logchlinc, ln_logchlbal, & 205 & ln_pco2inc, ln_fco2inc, & 194 206 & ln_asmdin, ln_asmiau, & 195 207 & nitbkg, nitdin, nitiaustr, nitiaufin, niaufn, & … … 229 241 WRITE(numout,*) ' Logical switch for applying logchl increments ln_logchlinc = ', ln_logchlinc 230 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 231 245 WRITE(numout,*) ' Logical switch for Incremental Analysis Updating (IAU) ln_asmiau = ', ln_asmiau 232 246 WRITE(numout,*) ' Timestep of background in [0,nitendnit0001] nitbkg = ', nitbkg … … 287 301 IF ( ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) & 288 302 & .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ).OR. & 289 & ( ln_logchlinc ) )) &303 & ( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) )) & 290 304 & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', & 291 & ' and ln_logchlinc is set to .true.', &305 & ' ln_logchlinc, ln_pco2inc, and ln_fco2inc is set to .true.', & 292 306 & ' but ln_asmdin and ln_asmiau are both set to .false. :', & 293 307 & ' Inconsistent options') … … 298 312 299 313 IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) & 300 & .AND.( .NOT. ln_logchlinc ) ) &314 & .AND.( .NOT. ln_logchlinc ).AND.( .NOT. ln_pco2inc ).AND.( .NOT. ln_fco2inc ) ) & 301 315 & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', & 302 & ' and ln_logchlinc are set to .false. :', &316 & ' ln_logchlinc, ln_pco2inc, and ln_fco2inc are set to .false. :', & 303 317 & ' The assimilation increments are not applied') 304 318 … … 321 335 & ' the cycle interval') 322 336 323 IF ( ( ln_balwri ).AND.( .NOT. ( ln_logchlinc )) ) THEN324 CALL ctl_warn( ' Balancing increments are only calculated for logchl ', &325 & ' Not assimilating logchl, so ln_balwri will be set to .false.')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.') 326 340 ln_balwri = .FALSE. 327 341 ENDIF 328 342 329 IF ( ( ln_logchlbal ).AND.( .NOT. ( ln_logchlinc )) ) THEN343 IF ( ( ln_logchlbal ).AND.( .NOT. ln_logchlinc ) ) THEN 330 344 CALL ctl_warn( ' Cannot calculate logchl balancing increments if not assimilating logchl', & 331 345 & ' Not assimilating logchl, so ln_logchlbal will be set to .false.') 332 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' ) 333 351 ENDIF 334 352 … … 436 454 #endif 437 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 438 464 IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) & 439 & .OR.( ln_logchlinc ) ) THEN465 & .OR.( ln_logchlinc ).OR.( ln_pco2inc ).OR.( ln_fco2inc ) ) THEN 440 466 441 467 ! … … 518 544 ! to allow for differences in masks 519 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 520 559 ENDIF 521 560 … … 703 742 CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) ) 704 743 #endif 705 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 706 774 707 775 CALL iom_close( inum ) … … 1454 1522 ! 1455 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 = jpk1 1668 DO jk = jpk1, 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 1456 1767 1457 1768 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.