Changeset 9381 for branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90
- Timestamp:
- 2018-03-07T16:19: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/asmbgc.F90
r9333 r9381 62 62 & asm_pco2_bal 63 63 USE par_medusa ! MEDUSA parameters 64 USE mocsy_mainmod, ONLY: & ! MEDUSA carbonate system 64 USE mocsy_wrapper, ONLY: & ! MEDUSA carbonate system 65 & mocsy_interface 66 USE mocsy_mainmod, ONLY: & ! fCO2 to pCO2 conversion 65 67 & f2pCO2 66 68 USE sms_medusa, ONLY: & ! MEDUSA diagnostic variables … … 177 179 #endif 178 180 181 # include "domzgr_substitute.h90" 182 179 183 CONTAINS 180 184 … … 215 219 ENDIF 216 220 217 IF ( ( ln_balwri ).AND.( .NOT. ln_phytobal ).AND. & 218 & ( .NOT. ln_spco2inc ).AND.( .NOT. ln_sfco2inc ).AND.( .NOT. ln_pphinc ) ) THEN 219 CALL ctl_warn( ' Balancing increments are only calculated for ocean colour', & 220 & ' with ln_phytobal, or for pCO2, fCO2, or pH.', & 221 & ' Not the case, so ln_balwri will be set to .false.') 221 IF ( ( ln_balwri ).AND. & 222 & ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. & 223 & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_schltotinc ).AND. & 224 & ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. & 225 & ( .NOT. ln_slphynoninc ).AND.( .NOT. ln_plchltotinc ).AND. & 226 & ( .NOT. ln_pchltotinc ).AND.( .NOT. ln_pphinc ).AND. & 227 & ( .NOT. ln_spco2inc ).AND.( .NOT. ln_sfco2inc ) ) THEN 228 CALL ctl_warn( ' Balancing increments not being calculated', & 229 & ' so ln_balwri will be set to .false.') 222 230 ln_balwri = .FALSE. 223 231 ENDIF … … 664 672 CALL iom_rstput( kt, kt, inum, 'rdastp' , zdate ) 665 673 666 IF ( ln_slchltotinc ) THEN 674 IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. & 675 & ln_schltotinc .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. & 676 & ln_slphynoninc ) THEN 667 677 #if defined key_medusa 668 CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chn', phyto2d_balinc(:,:,:,jpchn) ) 669 CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chd', phyto2d_balinc(:,:,:,jpchd) ) 678 CALL iom_rstput( kt, kt, inum, 'phy2d_chn', phyto2d_balinc(:,:,:,jpchn) ) 679 CALL iom_rstput( kt, kt, inum, 'phy2d_chd', phyto2d_balinc(:,:,:,jpchd) ) 680 CALL iom_rstput( kt, kt, inum, 'phy2d_phn', phyto2d_balinc(:,:,:,jpphn) ) 681 CALL iom_rstput( kt, kt, inum, 'phy2d_phd', phyto2d_balinc(:,:,:,jpphd) ) 682 CALL iom_rstput( kt, kt, inum, 'phy2d_pds', phyto2d_balinc(:,:,:,jppds) ) 670 683 IF ( ln_phytobal ) THEN 671 CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phn', phyto2d_balinc(:,:,:,jpphn) ) 672 CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phd', phyto2d_balinc(:,:,:,jpphd) ) 673 CALL iom_rstput( kt, kt, inum, 'logchl_balinc_pds', phyto2d_balinc(:,:,:,jppds) ) 674 CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zmi', phyto2d_balinc(:,:,:,jpzmi) ) 675 CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zme', phyto2d_balinc(:,:,:,jpzme) ) 676 CALL iom_rstput( kt, kt, inum, 'logchl_balinc_din', phyto2d_balinc(:,:,:,jpdin) ) 677 CALL iom_rstput( kt, kt, inum, 'logchl_balinc_sil', phyto2d_balinc(:,:,:,jpsil) ) 678 CALL iom_rstput( kt, kt, inum, 'logchl_balinc_fer', phyto2d_balinc(:,:,:,jpfer) ) 679 CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', phyto2d_balinc(:,:,:,jpdet) ) 680 CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dtc', phyto2d_balinc(:,:,:,jpdtc) ) 681 CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', phyto2d_balinc(:,:,:,jpdic) ) 682 CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', phyto2d_balinc(:,:,:,jpalk) ) 683 CALL iom_rstput( kt, kt, inum, 'logchl_balinc_oxy', phyto2d_balinc(:,:,:,jpoxy) ) 684 CALL iom_rstput( kt, kt, inum, 'phy2d_zmi', phyto2d_balinc(:,:,:,jpzmi) ) 685 CALL iom_rstput( kt, kt, inum, 'phy2d_zme', phyto2d_balinc(:,:,:,jpzme) ) 686 CALL iom_rstput( kt, kt, inum, 'phy2d_din', phyto2d_balinc(:,:,:,jpdin) ) 687 CALL iom_rstput( kt, kt, inum, 'phy2d_sil', phyto2d_balinc(:,:,:,jpsil) ) 688 CALL iom_rstput( kt, kt, inum, 'phy2d_fer', phyto2d_balinc(:,:,:,jpfer) ) 689 CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jpdet) ) 690 CALL iom_rstput( kt, kt, inum, 'phy2d_dtc', phyto2d_balinc(:,:,:,jpdtc) ) 691 CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jpdic) ) 692 CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jpalk) ) 693 CALL iom_rstput( kt, kt, inum, 'phy2d_oxy', phyto2d_balinc(:,:,:,jpoxy) ) 684 694 ENDIF 685 695 #elif defined key_hadocc 686 CALL iom_rstput( kt, kt, inum, ' logchl_balinc_phy', phyto2d_balinc(:,:,:,jp_had_phy) )696 CALL iom_rstput( kt, kt, inum, 'phy2d_phy', phyto2d_balinc(:,:,:,jp_had_phy) ) 687 697 IF ( ln_phytobal ) THEN 688 CALL iom_rstput( kt, kt, inum, ' logchl_balinc_nut', phyto2d_balinc(:,:,:,jp_had_nut) )689 CALL iom_rstput( kt, kt, inum, ' logchl_balinc_zoo', phyto2d_balinc(:,:,:,jp_had_zoo) )690 CALL iom_rstput( kt, kt, inum, ' logchl_balinc_det', phyto2d_balinc(:,:,:,jp_had_pdn) )691 CALL iom_rstput( kt, kt, inum, ' logchl_balinc_dic', phyto2d_balinc(:,:,:,jp_had_dic) )692 CALL iom_rstput( kt, kt, inum, ' logchl_balinc_alk', phyto2d_balinc(:,:,:,jp_had_alk) )698 CALL iom_rstput( kt, kt, inum, 'phy2d_nut', phyto2d_balinc(:,:,:,jp_had_nut) ) 699 CALL iom_rstput( kt, kt, inum, 'phy2d_zoo', phyto2d_balinc(:,:,:,jp_had_zoo) ) 700 CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jp_had_pdn) ) 701 CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jp_had_dic) ) 702 CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jp_had_alk) ) 693 703 ENDIF 704 #endif 705 ENDIF 706 707 IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN 708 #if defined key_medusa 709 CALL iom_rstput( kt, kt, inum, 'phy3d_chn', phyto3d_balinc(:,:,:,jpchn) ) 710 CALL iom_rstput( kt, kt, inum, 'phy3d_chd', phyto3d_balinc(:,:,:,jpchd) ) 711 CALL iom_rstput( kt, kt, inum, 'phy3d_phn', phyto3d_balinc(:,:,:,jpphn) ) 712 CALL iom_rstput( kt, kt, inum, 'phy3d_phd', phyto3d_balinc(:,:,:,jpphd) ) 713 CALL iom_rstput( kt, kt, inum, 'phy3d_pds', phyto3d_balinc(:,:,:,jppds) ) 714 #elif defined key_hadocc 715 CALL iom_rstput( kt, kt, inum, 'phy3d_phy', phyto3d_balinc(:,:,:,jp_had_phy) ) 694 716 #endif 695 717 ENDIF … … 697 719 IF ( ln_spco2inc ) THEN 698 720 #if defined key_medusa 699 CALL iom_rstput( kt, kt, inum, 'pco2_ balinc_dic', pco2_balinc(:,:,:,jpdic) )700 CALL iom_rstput( kt, kt, inum, 'pco2_ balinc_alk', pco2_balinc(:,:,:,jpalk) )721 CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jpdic) ) 722 CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jpalk) ) 701 723 #elif defined key_hadocc 702 CALL iom_rstput( kt, kt, inum, 'pco2_ balinc_dic', pco2_balinc(:,:,:,jp_had_dic) )703 CALL iom_rstput( kt, kt, inum, 'pco2_ balinc_alk', pco2_balinc(:,:,:,jp_had_alk) )724 CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jp_had_dic) ) 725 CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jp_had_alk) ) 704 726 #endif 705 727 ELSE IF ( ln_sfco2inc ) THEN 706 728 #if defined key_medusa 707 CALL iom_rstput( kt, kt, inum, 'fco2_ balinc_dic', pco2_balinc(:,:,:,jpdic) )708 CALL iom_rstput( kt, kt, inum, 'fco2_ balinc_alk', pco2_balinc(:,:,:,jpalk) )729 CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jpdic) ) 730 CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jpalk) ) 709 731 #elif defined key_hadocc 710 CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', pco2_balinc(:,:,:,jp_had_dic) ) 711 CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', pco2_balinc(:,:,:,jp_had_alk) ) 732 CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jp_had_dic) ) 733 CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jp_had_alk) ) 734 #endif 735 ENDIF 736 737 IF ( ln_pphinc ) THEN 738 #if defined key_medusa 739 CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jpdic) ) 740 CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jpalk) ) 741 #elif defined key_hadocc 742 CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jp_had_dic) ) 743 CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jp_had_alk) ) 712 744 #endif 713 745 ENDIF … … 944 976 REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights 945 977 ! 946 INTEGER :: ji, jj, jk ! Loop counters 947 INTEGER :: it ! Index 948 REAL(wp) :: zincwgt ! IAU weight for time step 949 REAL(wp) :: zfrac_chn ! Fraction of jpchn 950 REAL(wp) :: zfrac_chd ! Fraction of jpchd 951 REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc ! Chlorophyll increments 952 REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl ! Chlorophyll background 978 INTEGER :: ji, jj, jk ! Loop counters 979 INTEGER :: it ! Index 980 REAL(wp) :: zincwgt ! IAU weight for timestep 981 REAL(wp) :: zfrac_chn ! Fraction of jpchn 982 REAL(wp) :: zfrac_chd ! Fraction of jpchd 983 REAL(wp) :: zrat_phn_chn ! jpphn:jpchn ratio 984 REAL(wp) :: zrat_phd_chd ! jpphd:jpchd ratio 985 REAL(wp) :: zrat_pds_chd ! jppds:jpchd ratio 986 REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc ! Chlorophyll increments 987 REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl ! Chlorophyll background 953 988 !!------------------------------------------------------------------------ 954 989 … … 1006 1041 phyto3d_balinc(ji,jj,jk,jpchn) = chl_inc(ji,jj,jk) * zfrac_chn 1007 1042 phyto3d_balinc(ji,jj,jk,jpchd) = chl_inc(ji,jj,jk) * zfrac_chd 1043 zrat_phn_chn = tracer_bkg(ji,jj,jk,jpphn) / tracer_bkg(ji,jj,jk,jpchn) 1044 zrat_phd_chd = tracer_bkg(ji,jj,jk,jpphd) / tracer_bkg(ji,jj,jk,jpchd) 1045 zrat_pds_chd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpchd) 1046 phyto3d_balinc(ji,jj,jk,jpphn) = phyto3d_balinc(ji,jj,jk,jpchn) * zrat_phn_chn 1047 phyto3d_balinc(ji,jj,jk,jpphd) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_phd_chd 1048 phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_pds_chd 1008 1049 ENDIF 1009 1050 END DO … … 1192 1233 1193 1234 #if defined key_medusa 1194 ! Account for logchlbalancing if required1195 IF ( ln_ slchltotinc .AND. ln_phytobal ) THEN1235 ! Account for phytoplankton balancing if required 1236 IF ( ln_phytobal ) THEN 1196 1237 dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic) 1197 1238 alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_balinc(:,:,1,jpalk) … … 1206 1247 1207 1248 #elif defined key_hadocc 1208 ! Account for slchltotbalancing if required1209 IF ( ln_ slchltotinc .AND. ln_phytobal ) THEN1249 ! Account for phytoplankton balancing if required 1250 IF ( ln_phytobal ) THEN 1210 1251 dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto2d_balinc(:,:,1,jp_had_dic) 1211 1252 alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto2d_balinc(:,:,1,jp_had_alk) … … 1379 1420 !! ** Purpose : Apply the pH assimilation increments. 1380 1421 !! 1381 !! ** Method : Calculate increments to state variables using carbon1382 !! balancing.1422 !! ** Method : Calculate increments to DIC and alkalinity from pH 1423 !! Use a similar approach to the pCO2 scheme 1383 1424 !! Direct initialization or Incremental Analysis Updating. 1384 1425 !! 1385 1426 !! ** Action : 1386 1427 !!------------------------------------------------------------------------ 1387 INTEGER, INTENT(IN) :: kt! Current time step1388 LOGICAL, INTENT(IN) :: ll_asmdin! Flag for direct initialisation1389 LOGICAL, INTENT(IN) :: ll_asmiau! Flag for incremental analysis update1390 INTEGER, INTENT(IN) :: kcycper! Dimension of pwgtiau1391 REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights1392 LOGICAL, INTENT(IN) :: ll_trainc! Flag for T/S increments1428 INTEGER, INTENT(IN) :: kt ! Current time step 1429 LOGICAL, INTENT(IN) :: ll_asmdin ! Flag for direct initialisation 1430 LOGICAL, INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update 1431 INTEGER, INTENT(IN) :: kcycper ! Dimension of pwgtiau 1432 REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau ! IAU weights 1433 LOGICAL, INTENT(IN) :: ll_trainc ! Flag for T/S increments 1393 1434 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments 1394 1435 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments 1395 !!------------------------------------------------------------------------ 1396 1397 CALL ctl_stop( ' pH balancing not yet implemented' ) 1398 1399 ! See solve_at_general line 281 of mocsy_phsolvers.F90 1400 ! 1401 ! Or, call to mocsy_interface line 158 of carb_chem.F90 1402 ! Input variables (rest are output) are: 1403 ! ztmp(ji,jj),zsal(ji,jj),zalk(ji,jj), & 1404 ! zdic(ji,jj),zsil(ji,jj),zpho(ji,jj), & 1405 ! f_pp0(ji,jj),fsdept(ji,jj,jk), & 1406 ! gphit(ji,jj),f_kw660(ji,jj), & 1407 ! f_xco2a(ji,jj),1 1408 ! 1409 ! ztmp = tsn(:,:,:,jp_tem) 1410 ! zsal = tsn(:,:,:,jp_sal) 1411 ! zalk = trn(:,:,:,jpalk) 1412 ! zdic = trn(:,:,:,jpdic) 1413 ! zsil = trn(:,:,:,jpsil) 1414 ! zpho = trn(:,:,:,jpdin) / 16.0 1415 ! f_pp0 = 1.0 1416 ! fsdept = fsdept(:,:,:) 1417 ! gphit = gphit(:,:) 1418 ! f_kw660 = 1.0 1419 ! f_xco2a = f_xco2a(:,:) 1420 1436 1437 REAL(wp) :: zsearch = 10.0 ! Increment to DIC/alk in pH calculation 1438 REAL(wp) :: DIC_IN, ALK_IN ! DIC/alk in pH calculation 1439 REAL(wp) :: PH_OUT ! pH from pH calculation 1440 REAL(wp) :: ph_bkg ! pH from pH calculation 1441 REAL(wp) :: ph_dic ! pH from pH calculation 1442 REAL(wp) :: ph_alk ! pH from pH calculation 1443 REAL(wp) :: dph_ddic ! Change in pH wrt DIC 1444 REAL(wp) :: dph_dalk ! Change in pH wrt alk 1445 REAL(wp) :: weight ! Increment weighting 1446 REAL(wp) :: zincwgt ! IAU weight for current time step 1447 REAL(wp), DIMENSION(jpi,jpj,jpk) :: dic_bkg_temp ! DIC background 1448 REAL(wp), DIMENSION(jpi,jpj,jpk) :: alk_bkg_temp ! Alkalinity background 1449 REAL(wp), DIMENSION(jpi,jpj,jpk) :: din_bkg_temp ! DIN background 1450 REAL(wp), DIMENSION(jpi,jpj,jpk) :: sil_bkg_temp ! Silicate background 1451 REAL(wp), DIMENSION(jpi,jpj,jpk) :: tem_bkg_temp ! Temperature background 1452 REAL(wp), DIMENSION(jpi,jpj,jpk) :: sal_bkg_temp ! Salinity background 1453 REAL(wp), DIMENSION(19) :: dummy_out ! Dummy variables 1454 INTEGER :: ji, jj, jk, jx ! Loop counters 1455 INTEGER :: it ! Index 1456 !!------------------------------------------------------------------------ 1457 1458 #if ! (defined key_medusa && defined key_foam_medusa) 1459 CALL ctl_stop( ' pH balancing only implemented for MEDUSA' ) 1460 #else 1461 1462 IF ( kt <= nit000 ) THEN 1463 1464 !---------------------------------------------------------------------- 1465 ! DIC and alkalinity balancing 1466 !---------------------------------------------------------------------- 1467 1468 ! Account for physics assimilation if required 1469 IF ( ll_trainc ) THEN 1470 tem_bkg_temp(:,:,:) = tsn(:,:,:,1) + pt_bkginc(:,:,:) 1471 sal_bkg_temp(:,:,:) = tsn(:,:,:,2) + ps_bkginc(:,:,:) 1472 ELSE 1473 tem_bkg_temp(:,:,:) = tsn(:,:,:,1) 1474 sal_bkg_temp(:,:,:) = tsn(:,:,:,2) 1475 ENDIF 1476 1477 ! Account for phytoplankton balancing if required 1478 IF ( ln_phytobal ) THEN 1479 dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto2d_balinc(:,:,:,jpdic) 1480 alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto2d_balinc(:,:,:,jpalk) 1481 din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto2d_balinc(:,:,:,jpdin) 1482 sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto2d_balinc(:,:,:,jpsil) 1483 ELSE 1484 dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) 1485 alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) 1486 din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) 1487 sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) 1488 ENDIF 1489 1490 ! Account for pCO2 balancing if required 1491 IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN 1492 dic_bkg_temp(:,:,:) = dic_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpdic) 1493 alk_bkg_temp(:,:,:) = alk_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpalk) 1494 ENDIF 1495 1496 ! Loop over grid points and calculate dpH/dDIC and dpH/dAlk 1497 ! This requires three calls to the MOCSY carbonate package 1498 ! One to find pH at background DIC and alkalinity 1499 ! One to find pH when DIC is increased by 10 1500 ! One to find pH when alkalinity is increased by 10 1501 ! Then calculate the gradients 1502 1503 DO jk = 1, jpk 1504 DO jj = 1, jpj 1505 DO ji = 1, jpi 1506 1507 IF ( pph_bkginc(ji,jj,jk) /= 0.0 ) THEN 1508 1509 DO jx = 1, 3 1510 IF ( jx == 1 ) THEN 1511 DIC_IN = dic_bkg_temp(ji,jj,jk) 1512 ALK_IN = alk_bkg_temp(ji,jj,jk) 1513 ELSE IF ( jx == 2 ) THEN 1514 DIC_IN = dic_bkg_temp(ji,jj,jk) + zsearch 1515 ALK_IN = alk_bkg_temp(ji,jj,jk) 1516 ELSE IF ( jx == 3 ) THEN 1517 DIC_IN = dic_bkg_temp(ji,jj,jk) 1518 ALK_IN = alk_bkg_temp(ji,jj,jk) + zsearch 1519 ENDIF 1520 1521 CALL mocsy_interface(tsn(ji,jj,jk,jp_tem), & ! IN: temperature 1522 & tsn(ji,jj,jk,jp_sal), & ! salinity 1523 & ALK_IN, & ! alkalinity 1524 & DIC_IN, & ! DIC 1525 & sil_bkg_temp(ji,jj,jk), & ! silicate 1526 & din_bkg_temp(ji,jj,jk)/16.0, & ! phosphate (Redfield from DIN) 1527 & 1.0, & ! atmospheric pressure (dummy) 1528 & fsdept(ji,jj,jk), & ! depth 1529 & gphit(ji,jj), & ! latitude 1530 & 1.0, & ! gas transfer (dummy) 1531 & 1.0, & ! atmospheric xCO2 (dummy) 1532 & 1, & ! number of points 1533 & PH_OUT, & ! OUT: pH 1534 & dummy_out(1), dummy_out(2), & ! stuff we don't care about 1535 & dummy_out(3), dummy_out(4), & 1536 & dummy_out(5), dummy_out(6), & 1537 & dummy_out(7), dummy_out(8), & 1538 & dummy_out(9), dummy_out(10), & 1539 & dummy_out(11), dummy_out(12), & 1540 & dummy_out(13), dummy_out(14), & 1541 & dummy_out(15), dummy_out(16), & 1542 & dummy_out(17), dummy_out(18), & 1543 & dummy_out(19) ) 1544 1545 IF ( jx == 1 ) THEN 1546 ph_bkg = PH_OUT 1547 ELSE IF ( jx == 2 ) THEN 1548 ph_dic = PH_OUT 1549 ELSE IF ( jx == 3 ) THEN 1550 ph_alk = PH_OUT 1551 ENDIF 1552 END DO 1553 1554 dph_ddic = (ph_dic - ph_bkg) / zsearch 1555 dph_dalk = (ph_alk - ph_bkg) / zsearch 1556 weight = pph_bkginc(ji,jj,jk) / ( (dph_ddic*dph_ddic) + (dph_dalk*dph_dalk) ) 1557 1558 ph_balinc(ji,jj,jk,jpdic) = weight * dph_ddic 1559 ph_balinc(ji,jj,jk,jpalk) = weight * dph_dalk 1560 1561 ENDIF 1562 1563 END DO 1564 END DO 1565 END DO 1566 1567 ENDIF 1568 1569 IF ( ll_asmiau ) THEN 1570 1571 !-------------------------------------------------------------------- 1572 ! Incremental Analysis Updating 1573 !-------------------------------------------------------------------- 1574 1575 IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 1576 1577 it = kt - nit000 + 1 1578 zincwgt = pwgtiau(it) ! IAU weight for the current time step 1579 ! note this is not a tendency so should not be divided by rdt 1580 1581 IF(lwp) THEN 1582 WRITE(numout,*) 1583 WRITE(numout,*) 'ph_asm_inc : pH IAU at time step = ', & 1584 & kt,' with IAU weight = ', pwgtiau(it) 1585 WRITE(numout,*) '~~~~~~~~~~~~' 1586 ENDIF 1587 1588 ! Update the biogeochemical variables 1589 ! Add directly to trn and trb, rather than to tra, because tra gets 1590 ! reset to zero at the start of trc_stp, called after this routine 1591 DO jk = 1, jpkm1 1592 trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + & 1593 & ph_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 1594 trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + & 1595 & ph_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 1596 END DO 1597 1598 ! Do not deallocate arrays - needed by asm_bgc_bal_wri 1599 ! which is called at end of model run 1600 1601 ENDIF 1602 1603 ELSEIF ( ll_asmdin ) THEN 1604 1605 !-------------------------------------------------------------------- 1606 ! Direct Initialization 1607 !-------------------------------------------------------------------- 1608 1609 IF ( kt == nitdin_r ) THEN 1610 1611 neuler = 0 ! Force Euler forward step 1612 1613 ! Initialize the now fields with the background + increment 1614 ! Background currently is what the model is initialised with 1615 CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pH assimilation', & 1616 & ' Background state is taken from model rather than background file' ) 1617 trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & 1618 & ph_balinc(:,:,:,jp_msa0:jp_msa1) 1619 trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) 1620 1621 ! Do not deallocate arrays - needed by asm_bgc_bal_wri 1622 ! which is called at end of model run 1623 ENDIF 1624 ! 1625 ENDIF 1626 #endif 1421 1627 ! 1422 1628 END SUBROUTINE ph_asm_inc
Note: See TracChangeset
for help on using the changeset viewer.
