Changeset 367
- Timestamp:
- 2011-08-01T11:25:14+02:00 (14 years ago)
- Location:
- branches/ORCHIDEE_EXT/ORCHIDEE_OL
- Files:
-
- 3 edited
-
dim2_driver.f90 (modified) (22 diffs)
-
forcesoil.f90 (modified) (18 diffs)
-
teststomate.f90 (modified) (21 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/ORCHIDEE_EXT/ORCHIDEE_OL/dim2_driver.f90
r258 r367 121 121 122 122 CALL init_para(.FALSE.) 123 CALL init_timer 123 124 124 125 ! driver only for process root … … 686 687 ! This means loading the prognostic variables from the restart file. 687 688 !- 688 IF (is_root_prc) & 689 ALLOCATE(fluxsens_g(iim_g,jjm_g)) 689 Flag=.FALSE. 690 690 IF (is_root_prc) THEN 691 ALLOCATE(fluxsens_g(iim_g,jjm_g)) 691 692 var_name= 'fluxsens' 692 693 CALL restget & 693 694 & (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., fluxsens_g) 694 695 IF (ALL(fluxsens_g(:,:) == val_exp)) THEN 695 fluxsens_g(:,:) = zero 696 Flag=.TRUE. 697 ELSE 698 Flag=.FALSE. 696 699 ENDIF 697 ENDIF 698 CALL scatter2D(fluxsens_g,fluxsens) 699 IF (is_root_prc) & 700 DEALLOCATE(fluxsens_g) 701 !- 702 IF (is_root_prc) & 703 ALLOCATE(vevapp_g(iim_g,jjm_g)) 700 ELSE 701 ALLOCATE(fluxsens_g(0,1)) 702 ENDIF 703 CALL bcast(Flag) 704 IF (.NOT. Flag) THEN 705 CALL scatter2D(fluxsens_g,fluxsens) 706 ELSE 707 fluxsens(:,:) = zero 708 ENDIF 709 DEALLOCATE(fluxsens_g) 710 !- 704 711 IF (is_root_prc) THEN 712 ALLOCATE(vevapp_g(iim_g,jjm_g)) 705 713 var_name= 'vevapp' 706 714 CALL restget & 707 715 & (rest_id, var_name, iim_g, jjm_g, 1, istp_old, .TRUE., vevapp_g) 708 716 IF (ALL(vevapp_g(:,:) == val_exp)) THEN 709 vevapp(:,:) = 0. 717 Flag=.TRUE. 718 ELSE 719 Flag=.FALSE. 710 720 ENDIF 711 ENDIF 712 CALL scatter2D(vevapp_g,vevapp) 713 IF (is_root_prc) & 714 DEALLOCATE(vevapp_g) 715 !- 716 IF (is_root_prc) & 717 ALLOCATE(old_zlev_g(iim_g,jjm_g)) 721 ELSE 722 ALLOCATE(vevapp_g(0,1)) 723 ENDIF 724 CALL bcast(Flag) 725 IF (.NOT. Flag) THEN 726 CALL scatter2D(vevapp_g,vevapp) 727 ELSE 728 vevapp(:,:) = zero 729 ENDIF 730 DEALLOCATE(vevapp_g) 731 !- 718 732 IF (is_root_prc) THEN 733 ALLOCATE(old_zlev_g(iim_g,jjm_g)) 719 734 var_name= 'zlev_old' 720 735 CALL restget & … … 725 740 Flag=.FALSE. 726 741 ENDIF 727 ENDIF 728 CALL scatter2D(old_zlev_g,old_zlev) 729 IF (is_root_prc) & 730 DEALLOCATE(old_zlev_g) 742 ELSE 743 ALLOCATE(old_zlev_g(0,1)) 744 ENDIF 731 745 CALL bcast(Flag) 732 IF ( Flag ) old_zlev(:,:)=zlev_vec(:,:) 733 !- 734 IF (is_root_prc) & 735 ALLOCATE(old_qair_g(iim_g,jjm_g)) 746 IF ( .NOT. Flag ) THEN 747 CALL scatter2D(old_zlev_g,old_zlev) 748 ELSE 749 old_zlev(:,:)=zlev_vec(:,:) 750 ENDIF 751 DEALLOCATE(old_zlev_g) 752 !- 736 753 IF (is_root_prc) THEN 754 ALLOCATE(old_qair_g(iim_g,jjm_g)) 737 755 var_name= 'qair_old' 738 756 CALL restget & … … 743 761 Flag=.FALSE. 744 762 ENDIF 745 ENDIF 746 CALL scatter2D(old_qair_g,old_qair) 747 IF (is_root_prc) & 748 DEALLOCATE(old_qair_g) 763 ELSE 764 ALLOCATE(old_qair_g(0,1)) 765 ENDIF 749 766 CALL bcast(Flag) 750 IF (Flag) old_qair(:,:) = qair_obs(:,:) 751 !- 752 IF (is_root_prc) & 753 ALLOCATE(old_eair_g(iim_g,jjm_g)) 767 IF ( .NOT. Flag ) THEN 768 CALL scatter2D(old_qair_g,old_qair) 769 ELSE 770 old_qair(:,:) = qair_obs(:,:) 771 ENDIF 772 DEALLOCATE(old_qair_g) 773 !- 754 774 IF (is_root_prc) THEN 775 ALLOCATE(old_eair_g(iim_g,jjm_g)) 755 776 var_name= 'eair_old' 756 777 CALL restget & … … 761 782 Flag=.FALSE. 762 783 ENDIF 763 ENDIF 764 CALL scatter2D(old_eair_g,old_eair) 765 IF (is_root_prc) & 766 DEALLOCATE(old_eair_g) 784 ELSE 785 ALLOCATE(old_eair_g(0,1)) 786 ENDIF 767 787 CALL bcast(Flag) 768 IF (Flag) THEN 788 IF ( .NOT. Flag ) THEN 789 CALL scatter2D(old_eair_g,old_eair) 790 ELSE 769 791 DO ik=1,nbindex 770 792 i=ilandindex(ik) … … 773 795 ENDDO 774 796 ENDIF 797 DEALLOCATE(old_eair_g) 775 798 !- 776 799 ! old density is also needed because we do not yet have the right pb 777 800 !- 778 801 !=> obsolète ??!! (tjrs calculé après forcing_read) 779 IF (is_root_prc) &780 ALLOCATE(for_rau_g(iim_g,jjm_g))781 802 IF (is_root_prc) THEN 803 ALLOCATE(for_rau_g(iim_g,jjm_g)) 782 804 var_name= 'rau_old' 783 805 CALL restget & … … 788 810 Flag=.FALSE. 789 811 ENDIF 790 ENDIF 791 CALL scatter2D(for_rau_g,for_rau) 792 IF (is_root_prc) & 793 DEALLOCATE(for_rau_g) 812 ELSE 813 ALLOCATE(for_rau_g(0,1)) 814 ENDIF 794 815 CALL bcast(Flag) 795 IF (Flag) THEN 816 IF ( .NOT. Flag ) THEN 817 CALL scatter2D(for_rau_g,for_rau) 818 ELSE 796 819 DO ik=1,nbindex 797 820 i=ilandindex(ik) … … 800 823 ENDDO 801 824 ENDIF 825 DEALLOCATE(for_rau_g) 802 826 !- 803 827 ! For this variable the restart is extracted by SECHIBA … … 809 833 ! This does not yield a correct restart in the case of relaxation 810 834 !- 811 IF (is_root_prc) &812 ALLOCATE(petAcoef_g(iim_g,jjm_g))813 835 IF (is_root_prc) THEN 836 ALLOCATE(petAcoef_g(iim_g,jjm_g)) 814 837 var_name= 'petAcoef' 815 838 CALL restget & … … 820 843 Flag=.FALSE. 821 844 ENDIF 845 ELSE 846 ALLOCATE(petAcoef_g(0,1)) 822 847 ENDIF 823 CALL scatter2D(petAcoef_g,petAcoef)824 IF (is_root_prc) &825 DEALLOCATE(petAcoef_g)826 848 CALL bcast(Flag) 827 IF (Flag) petAcoef(:,:) = zero 849 IF ( .NOT. Flag ) THEN 850 CALL scatter2D(petAcoef_g,petAcoef) 851 ELSE 852 petAcoef(:,:) = zero 853 ENDIF 854 DEALLOCATE(petAcoef_g) 828 855 !-- 829 IF (is_root_prc) &830 ALLOCATE(petBcoef_g(iim_g,jjm_g))831 856 IF (is_root_prc) THEN 857 ALLOCATE(petBcoef_g(iim_g,jjm_g)) 832 858 var_name= 'petBcoef' 833 859 CALL restget & … … 838 864 Flag=.FALSE. 839 865 ENDIF 866 ELSE 867 ALLOCATE(petBcoef_g(0,1)) 840 868 ENDIF 841 CALL scatter2D(petBcoef_g,petBcoef)842 IF (is_root_prc) &843 DEALLOCATE(petBcoef_g)844 869 CALL bcast(Flag) 845 IF (Flag) petBcoef(:,:) = old_eair(:,:) 870 IF ( .NOT. Flag ) THEN 871 CALL scatter2D(petBcoef_g,petBcoef) 872 ELSE 873 petBcoef(:,:) = old_eair(:,:) 874 ENDIF 875 DEALLOCATE(petBcoef_g) 846 876 !-- 847 IF (is_root_prc) &848 ALLOCATE(peqAcoef_g(iim_g,jjm_g))849 877 IF (is_root_prc) THEN 878 ALLOCATE(peqAcoef_g(iim_g,jjm_g)) 850 879 var_name= 'peqAcoef' 851 880 CALL restget & … … 856 885 Flag=.FALSE. 857 886 ENDIF 887 ELSE 888 ALLOCATE(peqAcoef_g(0,1)) 858 889 ENDIF 859 CALL scatter2D(peqAcoef_g,peqAcoef)860 IF (is_root_prc) &861 DEALLOCATE(peqAcoef_g)862 890 CALL bcast(Flag) 863 IF (Flag) peqAcoef(:,:) = zero 891 IF ( .NOT. Flag ) THEN 892 CALL scatter2D(peqAcoef_g,peqAcoef) 893 ELSE 894 peqAcoef(:,:) = zero 895 ENDIF 896 DEALLOCATE(peqAcoef_g) 864 897 !-- 865 IF (is_root_prc) &866 ALLOCATE(peqBcoef_g(iim_g,jjm_g))867 898 IF (is_root_prc) THEN 899 ALLOCATE(peqBcoef_g(iim_g,jjm_g)) 868 900 var_name= 'peqBcoef' 869 901 CALL restget & … … 874 906 Flag=.FALSE. 875 907 ENDIF 908 ELSE 909 ALLOCATE(peqBcoef_g(0,1)) 876 910 ENDIF 877 CALL scatter2D(peqBcoef_g,peqBcoef)878 IF (is_root_prc) &879 DEALLOCATE(peqBcoef_g)880 911 CALL bcast(Flag) 881 IF (Flag) peqBcoef(:,:) = old_qair(:,:) 912 IF ( .NOT. Flag ) THEN 913 CALL scatter2D(peqBcoef_g,peqBcoef) 914 ELSE 915 peqBcoef(:,:) = old_qair(:,:) 916 ENDIF 917 DEALLOCATE(peqBcoef_g) 882 918 ENDIF 883 919 !- … … 951 987 IF (longprint) THEN 952 988 WRITE(numout,*) "dim2_driver 0 ",it_force 953 WRITE(numout,*) ">> Index of land points =",kindex 989 WRITE(numout,*) ">> Index of land points =",kindex(1:nbindex) 954 990 WRITE(numout,*) "Lowest level wind speed North = ", & 955 991 & (/ ( u(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) … … 1160 1196 IF (longprint) THEN 1161 1197 WRITE(numout,*) "dim2_driver first_CALL ",it_force 1162 WRITE(numout,*) ">> Index of land points =",kindex 1198 WRITE(numout,*) ">> Index of land points =",kindex(1:nbindex) 1163 1199 WRITE(numout,*) "Lowest level wind speed North = ", & 1164 1200 & (/ ( for_u(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) … … 1233 1269 !------- 1234 1270 ! albedo 1235 IF (is_root_prc) & 1236 ALLOCATE(albedo_g(iim_g,jjm_g)) 1271 IF (is_root_prc) THEN 1272 ALLOCATE(albedo_g(iim_g,jjm_g)) 1273 ELSE 1274 ALLOCATE(albedo_g(0,1)) 1275 ENDIF 1237 1276 ! 1238 1277 IF (is_root_prc) THEN … … 1246 1285 ENDIF 1247 1286 ENDIF 1248 CALL scatter2D(albedo_g,albedo_vis)1249 1287 CALL bcast(Flag) 1250 IF (.NOT. Flag) albedo(:,:,1)=albedo_vis(:,:) 1288 IF ( .NOT. Flag ) THEN 1289 CALL scatter2D(albedo_g,albedo_vis) 1290 albedo(:,:,1)=albedo_vis(:,:) 1291 ELSE 1292 albedo_vis(:,:)=albedo(:,:,1) 1293 ENDIF 1251 1294 ! 1252 1295 IF (is_root_prc) THEN … … 1260 1303 ENDIF 1261 1304 ENDIF 1262 CALL scatter2D(albedo_g,albedo_nir)1263 1305 CALL bcast(Flag) 1264 IF (.NOT. Flag) albedo(:,:,2)=albedo_nir(:,:) 1306 IF ( .NOT. Flag ) THEN 1307 CALL scatter2D(albedo_g,albedo_nir) 1308 albedo(:,:,2)=albedo_nir(:,:) 1309 ELSE 1310 albedo_nir(:,:)=albedo(:,:,2) 1311 ENDIF 1265 1312 ! 1266 IF (is_root_prc) & 1267 DEALLOCATE(albedo_g) 1313 DEALLOCATE(albedo_g) 1268 1314 !-- 1269 1315 ! z0 1270 IF (is_root_prc) &1271 ALLOCATE(z0_g(iim_g,jjm_g))1272 1316 IF (is_root_prc) THEN 1317 ALLOCATE(z0_g(iim_g,jjm_g)) 1273 1318 var_name= 'z0' 1274 1319 CALL restget & … … 1279 1324 Flag=.FALSE. 1280 1325 ENDIF 1326 ELSE 1327 ALLOCATE(z0_g(0,1)) 1281 1328 ENDIF 1282 1329 CALL bcast(Flag) 1283 IF (.NOT. Flag) CALL scatter2D(z0_g,z0)1284 IF (is_root_prc) &1285 1330 IF (.NOT. Flag) & 1331 CALL scatter2D(z0_g,z0) 1332 DEALLOCATE(z0_g) 1286 1333 !------- 1287 1334 DO ik=1,nbindex … … 1380 1427 IF (longprint) THEN 1381 1428 WRITE(numout,*) "dim2_driver ",it_force 1382 WRITE(numout,*) ">> Index of land points =",kindex 1429 WRITE(numout,*) ">> Index of land points =",kindex(1:nbindex) 1383 1430 WRITE(numout,*) "Lowest level wind speed North = ", & 1384 1431 & (/ ( for_u(ilandindex(ik), jlandindex(ik)),ik=1,nbindex ) /) … … 1546 1593 !- 1547 1594 var_name = 'fluxsens' 1548 IF (is_root_prc) & 1549 ALLOCATE(fluxsens_g(iim_g,jjm_g)) 1595 IF (is_root_prc) THEN 1596 ALLOCATE(fluxsens_g(iim_g,jjm_g)) 1597 ELSE 1598 ALLOCATE(fluxsens_g(0,1)) 1599 ENDIF 1550 1600 CALL gather2D(fluxsens , fluxsens_g) 1551 1601 IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, fluxsens_g) 1552 IF (is_root_prc) & 1553 DEALLOCATE(fluxsens_g) 1602 DEALLOCATE(fluxsens_g) 1554 1603 1555 1604 var_name = 'vevapp' 1556 IF (is_root_prc) & 1557 ALLOCATE(vevapp_g(iim_g,jjm_g)) 1605 IF (is_root_prc) THEN 1606 ALLOCATE(vevapp_g(iim_g,jjm_g)) 1607 ELSE 1608 ALLOCATE(vevapp_g(0,1)) 1609 ENDIF 1558 1610 CALL gather2D( vevapp, vevapp_g) 1559 1611 IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, vevapp_g) 1560 IF (is_root_prc) & 1561 DEALLOCATE(vevapp_g) 1612 DEALLOCATE(vevapp_g) 1562 1613 1563 1614 var_name = 'zlev_old' 1564 IF (is_root_prc) & 1565 ALLOCATE(old_zlev_g(iim_g,jjm_g)) 1615 IF (is_root_prc) THEN 1616 ALLOCATE(old_zlev_g(iim_g,jjm_g)) 1617 ELSE 1618 ALLOCATE(old_zlev_g(0,1)) 1619 ENDIF 1566 1620 CALL gather2D( old_zlev, old_zlev_g) 1567 1621 IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, old_zlev_g) 1568 IF (is_root_prc) & 1569 DEALLOCATE(old_zlev_g) 1622 DEALLOCATE(old_zlev_g) 1570 1623 1571 1624 var_name = 'qair_old' 1572 IF (is_root_prc) & 1573 ALLOCATE(old_qair_g(iim_g,jjm_g)) 1625 IF (is_root_prc) THEN 1626 ALLOCATE(old_qair_g(iim_g,jjm_g)) 1627 ELSE 1628 ALLOCATE(old_qair_g(0,1)) 1629 ENDIF 1574 1630 CALL gather2D( old_qair, old_qair_g) 1575 1631 IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, old_qair_g) 1576 IF (is_root_prc) & 1577 DEALLOCATE(old_qair_g) 1632 DEALLOCATE(old_qair_g) 1578 1633 1579 1634 var_name = 'eair_old' 1580 IF (is_root_prc) & 1581 ALLOCATE(old_eair_g(iim_g,jjm_g)) 1635 IF (is_root_prc) THEN 1636 ALLOCATE(old_eair_g(iim_g,jjm_g)) 1637 ELSE 1638 ALLOCATE(old_eair_g(0,1)) 1639 ENDIF 1582 1640 CALL gather2D( old_eair, old_eair_g) 1583 1641 IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, old_eair_g) 1584 IF (is_root_prc) & 1585 DEALLOCATE(old_eair_g) 1642 DEALLOCATE(old_eair_g) 1586 1643 1587 1644 var_name = 'rau_old' 1588 IF (is_root_prc) & 1589 ALLOCATE(for_rau_g(iim_g,jjm_g)) 1645 IF (is_root_prc) THEN 1646 ALLOCATE(for_rau_g(iim_g,jjm_g)) 1647 ELSE 1648 ALLOCATE(for_rau_g(0,1)) 1649 ENDIF 1590 1650 CALL gather2D( for_rau, for_rau_g) 1591 1651 IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, for_rau_g) 1592 IF (is_root_prc) & 1593 DEALLOCATE(for_rau_g) 1652 DEALLOCATE(for_rau_g) 1594 1653 1595 IF (is_root_prc) & 1596 ALLOCATE(albedo_g(iim_g,jjm_g)) 1654 IF (is_root_prc) THEN 1655 ALLOCATE(albedo_g(iim_g,jjm_g)) 1656 ELSE 1657 ALLOCATE(albedo_g(0,1)) 1658 ENDIF 1597 1659 var_name= 'albedo_vis' 1598 1660 albedo_vis(:,:)=albedo(:,:,1) … … 1604 1666 CALL gather2D(albedo_nir,albedo_g) 1605 1667 IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, albedo_g) 1606 IF (is_root_prc) & 1607 DEALLOCATE(albedo_g) 1668 DEALLOCATE(albedo_g) 1608 1669 1609 IF (is_root_prc) & 1610 ALLOCATE(z0_g(iim_g,jjm_g)) 1670 IF (is_root_prc) THEN 1671 ALLOCATE(z0_g(iim_g,jjm_g)) 1672 ELSE 1673 ALLOCATE(z0_g(0,1)) 1674 ENDIF 1611 1675 var_name= 'z0' 1612 1676 CALL gather2D(z0,z0_g) 1613 1677 IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, z0_g) 1614 IF (is_root_prc) & 1615 DEALLOCATE(z0_g) 1678 DEALLOCATE(z0_g) 1616 1679 1617 1680 if (.NOT. is_watchout) THEN 1618 var_name = 'petAcoef' 1619 IF (is_root_prc) & 1620 ALLOCATE(petAcoef_g(iim_g,jjm_g)) 1681 var_name = 'petAcoef' 1682 IF (is_root_prc) THEN 1683 ALLOCATE(petAcoef_g(iim_g,jjm_g)) 1684 ELSE 1685 ALLOCATE(petAcoef_g(0,1)) 1686 ENDIF 1621 1687 CALL gather2D( petAcoef, petAcoef_g) 1622 1688 IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, petAcoef_g) 1623 IF (is_root_prc) & 1624 DEALLOCATE(petAcoef_g) 1689 DEALLOCATE(petAcoef_g) 1625 1690 1626 var_name = 'petBcoef' 1627 IF (is_root_prc) & 1628 ALLOCATE(petBcoef_g(iim_g,jjm_g)) 1691 var_name = 'petBcoef' 1692 IF (is_root_prc) THEN 1693 ALLOCATE(petBcoef_g(iim_g,jjm_g)) 1694 ELSE 1695 ALLOCATE(petBcoef_g(0,1)) 1696 ENDIF 1629 1697 CALL gather2D( petBcoef, petBcoef_g) 1630 1698 IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, petBcoef_g) 1631 IF (is_root_prc) & 1632 DEALLOCATE(petBcoef_g) 1699 DEALLOCATE(petBcoef_g) 1633 1700 1634 var_name = 'peqAcoef' 1635 IF (is_root_prc) & 1636 ALLOCATE(peqAcoef_g(iim_g,jjm_g)) 1701 var_name = 'peqAcoef' 1702 IF (is_root_prc) THEN 1703 ALLOCATE(peqAcoef_g(iim_g,jjm_g)) 1704 ELSE 1705 ALLOCATE(peqAcoef_g(0,1)) 1706 ENDIF 1637 1707 CALL gather2D( peqAcoef, peqAcoef_g) 1638 1708 IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, peqAcoef_g) 1639 IF (is_root_prc) & 1640 DEALLOCATE(peqAcoef_g) 1709 DEALLOCATE(peqAcoef_g) 1641 1710 1642 var_name = 'peqBcoef' 1643 IF (is_root_prc) & 1644 ALLOCATE(peqBcoef_g(iim_g,jjm_g)) 1711 var_name = 'peqBcoef' 1712 IF (is_root_prc) THEN 1713 ALLOCATE(peqBcoef_g(iim_g,jjm_g)) 1714 ELSE 1715 ALLOCATE(peqBcoef_g(0,1)) 1716 ENDIF 1645 1717 CALL gather2D( peqBcoef, peqBcoef_g) 1646 1718 IF(is_root_prc) CALL restput (rest_id, var_name, iim_g, jjm_g, 1, istp_old, peqBcoef_g) 1647 IF (is_root_prc) & 1648 DEALLOCATE(peqBcoef_g) 1719 DEALLOCATE(peqBcoef_g) 1649 1720 ENDIF 1650 1721 !- -
branches/ORCHIDEE_EXT/ORCHIDEE_OL/forcesoil.f90
r313 r367 20 20 IMPLICIT NONE 21 21 !- 22 CHARACTER(LEN=80) :: sto_restname_in,sto_restname_out ,var_name22 CHARACTER(LEN=80) :: sto_restname_in,sto_restname_out 23 23 INTEGER(i_std) :: iim,jjm 24 24 25 INTEGER(i_std),PARAMETER :: llm = 1 25 26 INTEGER(i_std) :: kjpindex 27 26 28 INTEGER(i_std) :: itau_dep,itau_len 27 29 CHARACTER(LEN=30) :: time_str 28 INTEGER(i_std) :: ier,iret29 30 REAL(r_std) :: dt_files 30 31 REAL(r_std) :: date0 31 32 INTEGER(i_std) :: rest_id_sto 32 INTEGER(i_std) :: ncfid 33 REAL(r_std) :: dt_force,dt_forcesoil 33 34 !- 35 CHARACTER(LEN=100) :: Cforcing_name 36 INTEGER :: Cforcing_id 37 INTEGER :: v_id 38 REAL(r_std) :: dt_forcesoil 34 39 INTEGER :: nparan 35 INTEGER,PARAMETER :: nparanmax=36 36 REAL(r_std) :: xbid1,xbid2 37 INTEGER(i_std) :: ibid 40 38 41 INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices 42 INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices_g 43 REAL(r_std),DIMENSION(:),ALLOCATABLE :: x_indices_g 44 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lon, lat 39 45 REAL(r_std),DIMENSION(llm) :: lev 40 REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: soilcarbon_input 41 REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: & 42 & carbon,control_moist,control_temp 43 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: & 44 & lon,lat,resp_hetero_soil,var_3d 45 REAL(r_std),DIMENSION(:),ALLOCATABLE :: & 46 & x_indices 47 REAL(r_std) :: time 48 INTEGER :: i,j,m,iatt,iv 46 47 48 INTEGER :: i,m,iatt,iv 49 50 CHARACTER(LEN=80) :: var_name 49 51 CHARACTER(LEN=400) :: taboo_vars 50 52 REAL(r_std),DIMENSION(1) :: xtmp … … 56 58 INTEGER,DIMENSION(varnbdim_max) :: vardims 57 59 LOGICAL :: l1d 60 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: var_3d 58 61 REAL(r_std) :: x_tmp 59 ! clay fraction60 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: clay61 !-62 62 ! string suffix indicating an index 63 63 CHARACTER(LEN=10) :: part_str 64 64 ! 65 CHARACTER(LEN=100) :: Cforcing_name 66 INTEGER :: Cforcing_id 67 INTEGER :: v_id 68 69 REAL(r_std),ALLOCATABLE :: clay_loc(:) 70 REAL(r_std),ALLOCATABLE :: soilcarbon_input_loc(:,:,:,:) 71 REAL(r_std),ALLOCATABLE :: control_temp_loc(:,:,:) 72 REAL(r_std),ALLOCATABLE :: control_moist_loc(:,:,:) 73 REAL(r_std),ALLOCATABLE :: carbon_loc(:,:,:) 74 INTEGER :: ierr 65 ! clay fraction 66 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:) :: clay_g 67 REAL(r_std),DIMENSION(:,:,:,:),ALLOCATABLE :: soilcarbon_input_g 68 REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: control_temp_g 69 REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: control_moist_g 70 REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: carbon_g 71 72 REAL(r_std),ALLOCATABLE :: clay(:) 73 REAL(r_std),ALLOCATABLE :: soilcarbon_input(:,:,:,:) 74 REAL(r_std),ALLOCATABLE :: control_temp(:,:,:) 75 REAL(r_std),ALLOCATABLE :: control_moist(:,:,:) 76 REAL(r_std),ALLOCATABLE :: carbon(:,:,:) 77 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: resp_hetero_soil 78 79 INTEGER(i_std) :: ier,iret 80 81 LOGICAL :: debug 82 75 83 !>> DS add for externalization 76 84 LOGICAL :: l_error … … 78 86 79 87 CALL Init_para(.FALSE.) 80 88 CALL init_timer 81 89 ! 82 90 ! DS : For externalization cause we decoupled forcesoil from ORCHIDEE … … 84 92 85 93 ! 1. Read the number of PFTs 86 CALL getin('NVM',nvm) 94 ! 95 !Config Key = NVM 96 !Config Desc = number of PFTs 97 !Config if = ANYTIME 98 !Config Def = 13 99 !Config Help = The number of vegetation types define by the user 100 !Config Units = NONE 101 CALL getin_p('NVM',nvm) 102 87 103 ! 2. Allocation 88 104 l_error = .FALSE. … … 97 113 98 114 ! 4.Reading of the conrrespondance table in the .def file 99 CALL getin('PFT_TO_MTC',pft_to_mtc) 115 ! 116 !Config Key = PFT_TO_MTC 117 !Config Desc = correspondance array linking a PFT to MTC 118 !Config if = ANYTIME 119 !Config Def = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 120 !Config Help = 121 !Config Units = NONE 122 CALL getin_p('PFT_TO_MTC',pft_to_mtc) 100 123 101 124 ! 4.1 if nothing is found, we use the standard configuration … … 147 170 148 171 ! 7.2 Initialisation 149 DO j= 1, nvm150 natural( j) = natural_mtc(pft_to_mtc(j))151 is_c4( j) = is_c4_mtc(pft_to_mtc(j))172 DO i= 1, nvm 173 natural(i) = natural_mtc(pft_to_mtc(i)) 174 is_c4(i) = is_c4_mtc(pft_to_mtc(i)) 152 175 ENDDO 153 176 177 !--------------------------------------------------------------------- 178 !- 179 ! set debug to have more information 180 !- 181 !Config Key = DEBUG_INFO 182 !Config Desc = Flag for debug information 183 !Config Def = n 184 !Config Help = This option allows to switch on the output of debug 185 !Config information without recompiling the code. 186 !- 187 debug = .FALSE. 188 CALL getin_p('DEBUG_INFO',debug) 189 ! 190 !Config Key = LONGPRINT 191 !Config Desc = ORCHIDEE will print more messages 192 !Config Def = n 193 !Config Help = This flag permits to print more debug messages in the run. 194 ! 195 long_print = .FALSE. 196 CALL getin_p('LONGPRINT',long_print) 154 197 !- 155 198 ! Stomate's restart files … … 158 201 sto_restname_in = 'stomate_start.nc' 159 202 CALL getin ('STOMATE_RESTART_FILEIN',sto_restname_in) 160 WRITE( *,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in)161 sto_restname_out = 'stomate_rest art.nc'203 WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in) 204 sto_restname_out = 'stomate_rest_out.nc' 162 205 CALL getin ('STOMATE_RESTART_FILEOUT',sto_restname_out) 163 WRITE( *,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out)164 !- 165 ! We need to know iim , jjm.206 WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out) 207 !- 208 ! We need to know iim_g, jjm. 166 209 ! Get them from the restart files themselves. 167 210 !- 168 iret = NF90_OPEN (sto_restname_in, NF90_NOWRITE, ncfid)169 iret = NF90_INQUIRE_DIMENSION ( ncfid,1,len=iim)170 iret = NF90_INQUIRE_DIMENSION ( ncfid,2,len=jjm)171 iret = NF90_CLOSE ( ncfid)211 iret = NF90_OPEN (sto_restname_in, NF90_NOWRITE, rest_id_sto) 212 iret = NF90_INQUIRE_DIMENSION (rest_id_sto,1,len=iim_g) 213 iret = NF90_INQUIRE_DIMENSION (rest_id_sto,2,len=jjm_g) 214 iret = NF90_CLOSE (rest_id_sto) 172 215 !- 173 216 ! Allocate longitudes and latitudes 174 217 !- 175 ALLOCATE (lon(iim ,jjm))176 ALLOCATE (lat(iim ,jjm))218 ALLOCATE (lon(iim_g,jjm_g)) 219 ALLOCATE (lat(iim_g,jjm_g)) 177 220 lon(:,:) = 0.0 178 221 lat(:,:) = 0.0 … … 180 223 !- 181 224 CALL restini & 182 & (sto_restname_in, iim , jjm, lon, lat, llm, lev, &225 & (sto_restname_in, iim_g, jjm_g, lon, lat, llm, lev, & 183 226 & sto_restname_out, itau_dep, date0, dt_files, rest_id_sto) 184 227 ENDIF … … 188 231 CALL bcast(date0) 189 232 !!! MM : à revoir : choix du calendrier dans forcesoil ?? Il est dans le restart de stomate ! 190 !CALL ioconf_calendar ('noleap')233 CALL ioconf_calendar ('noleap') 191 234 CALL ioget_calendar (one_year,one_day) 192 235 … … 197 240 ! open FORCESOIL's forcing file to read some basic info 198 241 !- 199 Cforcing_name = ' stomate_Cforcing.nc'242 Cforcing_name = 'NONE' 200 243 CALL getin ('STOMATE_CFORCING_NAME',Cforcing_name) 201 244 !- 202 ier = NF90_OPEN (TRIM(Cforcing_name),NF90_NOWRITE,Cforcing_id) 245 iret = NF90_OPEN (TRIM(Cforcing_name),NF90_NOWRITE,Cforcing_id) 246 IF (iret /= NF90_NOERR) THEN 247 CALL ipslerr (3,'forcesoil', & 248 & 'Could not open file : ', & 249 & Cforcing_name,'(Do you have forget it ?)') 250 ENDIF 203 251 !- 204 252 ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'kjpindex',x_tmp) 205 kjpindex= NINT(x_tmp)253 nbp_glo = NINT(x_tmp) 206 254 ier = NF90_GET_ATT (Cforcing_id,NF90_GLOBAL,'nparan',x_tmp) 207 255 nparan = NINT(x_tmp) 208 256 !- 209 ALLOCATE (indices (kjpindex))210 ALLOCATE (clay (kjpindex))211 !- 212 ALLOCATE (x_indices (kjpindex),stat=ier)257 ALLOCATE (indices_g(nbp_glo)) 258 ALLOCATE (clay_g(nbp_glo)) 259 !- 260 ALLOCATE (x_indices_g(nbp_glo),stat=ier) 213 261 ier = NF90_INQ_VARID (Cforcing_id,'index',v_id) 214 ier = NF90_GET_VAR (Cforcing_id,v_id,x_indices) 215 indices(:) = NINT(x_indices(:)) 216 DEALLOCATE (x_indices) 262 ier = NF90_GET_VAR (Cforcing_id,v_id,x_indices_g) 263 indices_g(:) = NINT(x_indices_g(:)) 264 WRITE(numout,*) mpi_rank,"indices globaux : ",indices_g 265 DEALLOCATE (x_indices_g) 217 266 !- 218 267 ier = NF90_INQ_VARID (Cforcing_id,'clay',v_id) 219 ier = NF90_GET_VAR (Cforcing_id,v_id,clay )268 ier = NF90_GET_VAR (Cforcing_id,v_id,clay_g) 220 269 !- 221 270 ! time step of forcesoil 222 271 !- 223 272 dt_forcesoil = one_year / FLOAT(nparan) 224 WRITE( *,*) 'time step (d): ',dt_forcesoil273 WRITE(numout,*) 'time step (d): ',dt_forcesoil 225 274 !- 226 275 ! read (and partially write) the restart file … … 258 307 l1d = ALL(vardims(1:varnbdim) == 1) 259 308 !---- 260 ALLOCATE( var_3d( kjpindex,vardims(3)), stat=ier)309 ALLOCATE( var_3d(nbp_glo,vardims(3)), stat=ier) 261 310 IF (ier /= 0) STOP 'ALLOCATION PROBLEM' 262 311 !---- read it … … 267 316 ELSE 268 317 CALL restget & 269 & (rest_id_sto, TRIM(varnames(iv)), kjpindex, vardims(3), &270 & 1, itau_dep, .TRUE., var_3d, "gather", kjpindex, indices)318 & (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), & 319 & 1, itau_dep, .TRUE., var_3d, "gather", nbp_glo, indices_g) 271 320 ENDIF 272 321 !---- write it … … 277 326 ELSE 278 327 CALL restput & 279 & (rest_id_sto, TRIM(varnames(iv)), kjpindex, vardims(3), &280 & 1, itau_dep, var_3d, 'scatter', kjpindex, indices)328 & (rest_id_sto, TRIM(varnames(iv)), nbp_glo, vardims(3), & 329 & 1, itau_dep, var_3d, 'scatter', nbp_glo, indices_g) 281 330 ENDIF 282 331 !---- … … 287 336 ! read soil carbon 288 337 !- 289 ALLOCATE(carbon (kjpindex,ncarb,nvm))290 carbon (:,:,:) = val_exp338 ALLOCATE(carbon_g(nbp_glo,ncarb,nvm)) 339 carbon_g(:,:,:) = val_exp 291 340 DO m = 1, nvm 292 341 WRITE (part_str, '(I2)') m 293 342 IF (m<10) part_str(1:1)='0' 294 343 var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str)) 295 CALL restget _p&296 & (rest_id_sto, var_name, kjpindex, ncarb , 1, itau_dep, &297 & .TRUE., carbon (:,:,m), 'gather', kjpindex, indices)298 IF (ALL(carbon (:,:,m) == val_exp)) carbon(:,:,m) = zero344 CALL restget & 345 & (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, & 346 & .TRUE., carbon_g(:,:,m), 'gather', nbp_glo, indices_g) 347 IF (ALL(carbon_g(:,:,m) == val_exp)) carbon_g(:,:,m) = zero 299 348 !-- do not write this variable: it will be modified. 300 349 ENDDO … … 304 353 WRITE(time_str,'(a)') '10000Y' 305 354 CALL getin('TIME_LENGTH', time_str) 355 write(numout,*) 'Number of years for carbon spinup : ',time_str 306 356 ! transform into itau 307 357 CALL tlen2itau(time_str, dt_forcesoil*one_year, date0, itau_len) 308 write( *,*) 'Number of time steps to do: ',itau_len358 write(numout,*) 'Number of time steps to do: ',itau_len 309 359 !- 310 360 ! read the rest of the forcing file and store forcing in an array. 311 361 ! We read an average year. 312 362 !- 313 ALLOCATE(soilcarbon_input (kjpindex,ncarb,nvm,nparan))314 ALLOCATE(control_temp (kjpindex,nlevs,nparan))315 ALLOCATE(control_moist (kjpindex,nlevs,nparan))363 ALLOCATE(soilcarbon_input_g(nbp_glo,ncarb,nvm,nparan)) 364 ALLOCATE(control_temp_g(nbp_glo,nlevs,nparan)) 365 ALLOCATE(control_moist_g(nbp_glo,nlevs,nparan)) 316 366 !- 317 367 ier = NF90_INQ_VARID (Cforcing_id,'soilcarbon_input',v_id) 318 ier = NF90_GET_VAR (Cforcing_id,v_id,soilcarbon_input )368 ier = NF90_GET_VAR (Cforcing_id,v_id,soilcarbon_input_g) 319 369 ier = NF90_INQ_VARID (Cforcing_id, 'control_moist',v_id) 320 ier = NF90_GET_VAR (Cforcing_id,v_id,control_moist )370 ier = NF90_GET_VAR (Cforcing_id,v_id,control_moist_g) 321 371 ier = NF90_INQ_VARID (Cforcing_id, 'control_temp',v_id) 322 ier = NF90_GET_VAR (Cforcing_id,v_id,control_temp )372 ier = NF90_GET_VAR (Cforcing_id,v_id,control_temp_g) 323 373 !- 324 374 ier = NF90_CLOSE (Cforcing_id) 325 375 !- 326 !MM Problem here with dpu which depends on soil type327 DO iv = 1, nbdl-1328 ! first 2.0 is dpu329 ! second 2.0 is average330 diaglev(iv) = 2.0/(2**(nbdl-1) -1) * ( ( 2**(iv-1) -1) + ( 2**(iv) -1) ) / 2.0331 ENDDO332 diaglev(nbdl) = 2.0333 !-334 ! For sequential use only, we must initialize data_para :376 ENDIF 377 CALL bcast(nparan) 378 CALL bcast(dt_forcesoil) 379 CALL bcast(iim_g) 380 CALL bcast(jjm_g) 381 call bcast(nbp_glo) 382 CALL bcast(itau_len) 383 ! 384 ! We must initialize data_para : 335 385 ! 336 386 ! 337 ENDIF 338 339 CALL bcast(iim) 340 CALL bcast(jjm) 341 call bcast(kjpindex) 342 CALL init_data_para(iim,jjm,kjpindex,indices) 343 387 CALL init_data_para(iim_g,jjm_g,nbp_glo,indices_g) 388 389 kjpindex=nbp_loc 390 jjm=jj_nb 391 iim=iim_g 392 IF (debug) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm 393 394 !--- 395 !--- Create the index table 396 !--- 397 !--- This job return a LOCAL kindex 398 !--- 399 ALLOCATE (indices(kjpindex),stat=ier) 400 CALL scatter(indices_g,indices) 401 indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g 402 WRITE(numout,*) mpi_rank,"indices locaux = ",indices(1:kjpindex) 344 403 !- 345 404 !- 346 405 ! there we go: time loop 347 406 !- 348 CALL bcast(nparan) 349 ALLOCATE(clay_loc(nbp_loc)) 350 ALLOCATE(soilcarbon_input_loc(nbp_loc,ncarb,nvm,nparan)) 351 ALLOCATE(control_temp_loc(nbp_loc,nlevs,nparan)) 352 ALLOCATE(control_moist_loc(nbp_loc,nlevs,nparan)) 353 ALLOCATE(carbon_loc(nbp_loc,ncarb,nvm)) 354 ALLOCATE(resp_hetero_soil(nbp_loc,nvm)) 407 ALLOCATE(clay(kjpindex)) 408 ALLOCATE(soilcarbon_input(kjpindex,ncarb,nvm,nparan)) 409 ALLOCATE(control_temp(kjpindex,nlevs,nparan)) 410 ALLOCATE(control_moist(kjpindex,nlevs,nparan)) 411 ALLOCATE(carbon(kjpindex,ncarb,nvm)) 412 ALLOCATE(resp_hetero_soil(kjpindex,nvm)) 355 413 iatt = 0 356 414 357 CALL bcast(itau_len) 358 CALL bcast(nparan) 359 CALL bcast(dt_forcesoil) 360 CALL Scatter(clay,clay_loc) 361 CALL Scatter(soilcarbon_input,soilcarbon_input_loc) 362 CALL Scatter(control_temp,control_temp_loc) 363 CALL Scatter(control_moist,control_moist_loc) 364 CALL Scatter(carbon,carbon_loc) 365 366 !!$ DS 16/06/2011 : calling the new_values of soilcarbon parameters before loop 367 ! 368 CALL getin_p('FRAC_CARB_AA',frac_carb_aa) 369 CALL getin_p('FRAC_CARB_AP',frac_carb_ap) 370 CALL getin_p('FRAC_CARB_SS',frac_carb_ss) 371 CALL getin_p('FRAC_CARB_SA',frac_carb_sa) 372 CALL getin_p('FRAC_CARB_SP',frac_carb_sp) 373 CALL getin_p('FRAC_CARB_PP',frac_carb_pp) 374 CALL getin_p('FRAC_CARB_PA',frac_carb_pa) 375 CALL getin_p('FRAC_CARB_PS',frac_carb_ps) 376 ! 377 CALL getin_p('ACTIVE_TO_PASS_CLAY_FRAC',active_to_pass_clay_frac) 378 CALL getin_p('CARBON_TAU_IACTIVE',carbon_tau_iactive) 379 CALL getin_p('CARBON_TAU_ISLOW',carbon_tau_islow) 380 CALL getin_p('CARBON_TAU_IPASSIVE',carbon_tau_ipassive) 381 CALL getin_p('FLUX_TOT_COEFF',flux_tot_coeff) 415 CALL Scatter(clay_g,clay) 416 CALL Scatter(soilcarbon_input_g,soilcarbon_input) 417 CALL Scatter(control_temp_g,control_temp) 418 CALL Scatter(control_moist_g,control_moist) 419 CALL Scatter(carbon_g,carbon) 420 421 !!$ DS 16/06/2011 : for externalization 422 ! 423 !Config Key = FRAC_CARB_AA 424 !Config Desc = frac carb coefficients from active pool: depends on clay content 425 !Config if = OK_STOMATE 426 !Config Def = 0.0 427 !Config Help = fraction of the active pool going to the active pool 428 !Config Units = NONE 429 CALL getin_p('FRAC_CARB_AA',frac_carb_aa) 430 ! 431 !Config Key = FRAC_CARB_AP 432 !Config Desc = frac carb coefficients from active pool: depends on clay content 433 !Config if = OK_STOMATE 434 !Config Def = 0.004 435 !Config Help = fraction of the active pool going to the passive pool 436 !Config Units = NONE 437 CALL getin_p('FRAC_CARB_AP',frac_carb_ap) 438 ! 439 !Config Key = FRAC_CARB_SS 440 !Config Desc = frac_carb_coefficients from slow pool 441 !Config if = OK_STOMATE 442 !Config Def = 0.0 443 !Config Help = fraction of the slow pool going to the slow pool 444 !Config Units = NONE 445 CALL getin_p('FRAC_CARB_SS',frac_carb_ss) 446 ! 447 !Config Key = FRAC_CARB_SA 448 !Config Desc = frac_carb_coefficients from slow pool 449 !Config if = OK_STOMATE 450 !Config Def = 0.42 451 !Config Help = fraction of the slow pool going to the active pool 452 !Config Units = NONE 453 CALL getin_p('FRAC_CARB_SA',frac_carb_sa) 454 ! 455 !Config Key = FRAC_CARB_SP 456 !Config Desc = frac_carb_coefficients from slow pool 457 !Config if = OK_STOMATE 458 !Config Def = 0.03 459 !Config Help = fraction of the slow pool going to the passive pool 460 !Config Units = NONE 461 CALL getin_p('FRAC_CARB_SP',frac_carb_sp) 462 ! 463 !Config Key = FRAC_CARB_PP 464 !Config Desc = frac_carb_coefficients from passive pool 465 !Config if = OK_STOMATE 466 !Config Def = 0.0 467 !Config Help = fraction of the passive pool going to the passive pool 468 !Config Units = NONE 469 CALL getin_p('FRAC_CARB_PP',frac_carb_pp) 470 ! 471 !Config Key = FRAC_CARB_PA 472 !Config Desc = frac_carb_coefficients from passive pool 473 !Config if = OK_STOMATE 474 !Config Def = 0.45 475 !Config Help = fraction of the passive pool going to the passive pool 476 !Config Units = NONE 477 CALL getin_p('FRAC_CARB_PA',frac_carb_pa) 478 ! 479 !Config Key = FRAC_CARB_PS 480 !Config Desc = frac_carb_coefficients from passive pool 481 !Config if = OK_STOMATE 482 !Config Def = 0.0 483 !Config Help = fraction of the passive pool going to the passive pool 484 !Config Units = NONE 485 CALL getin_p('FRAC_CARB_PS',frac_carb_ps) 486 ! 487 !Config Key = ACTIVE_TO_PASS_CLAY_FRAC 488 !Config Desc = 489 !Config if = OK_STOMATE 490 !Config Def = .68 491 !Config Help = 492 !Config Units = NONE 493 CALL getin_p('ACTIVE_TO_PASS_CLAY_FRAC',active_to_pass_clay_frac) 494 ! 495 !Config Key = CARBON_TAU_IACTIVE 496 !Config Desc = residence times in carbon pools 497 !Config if = OK_STOMATE 498 !Config Def = 0.149 499 !Config Help = 500 !Config Units = days [d] 501 CALL getin_p('CARBON_TAU_IACTIVE',carbon_tau_iactive) 502 ! 503 !Config Key = CARBON_TAU_ISLOW 504 !Config Desc = residence times in carbon pools 505 !Config if = OK_STOMATE 506 !Config Def = 5.48 507 !Config Help = 508 !Config Units = days [d] 509 CALL getin_p('CARBON_TAU_ISLOW',carbon_tau_islow) 510 ! 511 !Config Key = CARBON_TAU_IPASSIVE 512 !Config Desc = residence times in carbon pools 513 !Config if = OK_STOMATE 514 !Config Def = 241. 515 !Config Help = 516 !Config Units = days [d] 517 CALL getin_p('CARBON_TAU_IPASSIVE',carbon_tau_ipassive) 518 ! 519 !Config Key = FLUX_TOT_COEFF 520 !Config Desc = 521 !Config if = OK_STOMATE 522 !Config Def = 1.2, 1.4,.75 523 !Config Help = 524 !Config Units = days [d] 525 CALL getin_p('FLUX_TOT_COEFF',flux_tot_coeff) 382 526 383 527 DO i=1,itau_len … … 385 529 IF (iatt > nparan) iatt = 1 386 530 CALL soilcarbon & 387 & ( nbp_loc, dt_forcesoil, clay_loc, &388 & soilcarbon_input _loc(:,:,:,iatt), &389 & control_temp _loc(:,:,iatt), control_moist_loc(:,:,iatt), &390 & carbon _loc, resp_hetero_soil)531 & (kjpindex, dt_forcesoil, clay, & 532 & soilcarbon_input(:,:,:,iatt), & 533 & control_temp(:,:,iatt), control_moist(:,:,iatt), & 534 & carbon, resp_hetero_soil) 391 535 ENDDO 392 536 393 CALL Gather(carbon _loc,carbon)537 CALL Gather(carbon,carbon_g) 394 538 !- 395 539 ! write new carbon into restart file … … 400 544 IF (m<10) part_str(1:1)='0' 401 545 var_name = 'carbon_'//part_str(1:LEN_TRIM(part_str)) 402 CALL restput _p&403 & (rest_id_sto, var_name, kjpindex, ncarb , 1, itau_dep, &404 & carbon (:,:,m), 'scatter', kjpindex, indices)546 CALL restput & 547 & (rest_id_sto, var_name, nbp_glo, ncarb , 1, itau_dep, & 548 & carbon_g(:,:,m), 'scatter', nbp_glo, indices_g) 405 549 ENDDO 406 550 !- … … 409 553 ENDIF 410 554 #ifdef CPP_PARA 411 CALL MPI_FINALIZE(ier r)555 CALL MPI_FINALIZE(ier) 412 556 #endif 557 WRITE(numout,*) "End of forcesoil." 413 558 !-------------------- 414 559 END PROGRAM forcesoil -
branches/ORCHIDEE_EXT/ORCHIDEE_OL/teststomate.f90
r327 r367 21 21 USE slowproc 22 22 USE stomate 23 USE intersurf, ONLY: stom_define_history , intsurf_time23 USE intersurf, ONLY: stom_define_history, stom_ipcc_define_history, intsurf_time, l_first_intersurf, check_time 24 24 USE parallel 25 25 !- … … 28 28 ! Declarations 29 29 !- 30 INTEGER(i_std) :: vegax_id31 30 INTEGER(i_std) :: kjpij,kjpindex 32 31 REAL(r_std) :: dtradia,dt_force 32 33 33 INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices 34 34 INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indexveg … … 41 41 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: veget_max_force_x 42 42 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lai_force_x 43 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: lon,lat44 43 REAL(r_std),DIMENSION(:),ALLOCATABLE :: t2m,t2m_min,temp_sol 45 44 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: soiltemp,soilhum … … 53 52 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: qsintmax_x 54 53 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: co2_flux 55 !!$ >> DS add for externalised version 54 REAL(r_std),DIMENSION(:),ALLOCATABLE :: fco2_lu 55 56 INTEGER(i_std),DIMENSION(:),ALLOCATABLE :: indices_g 57 REAL(r_std),DIMENSION(:),ALLOCATABLE :: x_indices_g 58 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: x_neighbours_g 59 60 !!$ >> DS for the externalized version 56 61 REAL(r_std),DIMENSION(:),ALLOCATABLE :: hist_PFTaxis 57 62 !!$ >> DS 58 !!$ >> Add for correct 1.9.5.1 version 01/07/201159 ! TYPE(control_type), SAVE :: control_flags60 REAL(r_std),DIMENSION(:),ALLOCATABLE :: fco2_lu61 !!$>> DS62 63 63 64 64 INTEGER :: ier,iret … … 68 68 & dri_restname_in,dri_restname_out, & 69 69 & sec_restname_in,sec_restname_out, & 70 & sto_restname_in,sto_restname_out, stom_histname 71 INTEGER(i_std) :: iim,jjm 72 INTEGER(i_std),PARAMETER :: llm = 1 73 REAL(r_std),DIMENSION(llm) :: lev 74 REAL(r_std) :: dt_files 70 & sto_restname_in,sto_restname_out, & 71 & stom_histname, stom_ipcc_histname 72 INTEGER(i_std) :: iim,jjm,llm 73 REAL, ALLOCATABLE, DIMENSION(:,:) :: lon, lat 74 REAL, ALLOCATABLE, DIMENSION(:) :: lev 75 LOGICAL :: rectilinear 76 REAL, ALLOCATABLE, DIMENSION(:) :: lon_rect, lat_rect 77 REAL(r_std) :: dt 75 78 INTEGER(i_std) :: itau_dep,itau,itau_len,itau_step 76 79 REAL(r_std) :: date0 77 80 INTEGER(i_std) :: rest_id_sec,rest_id_sto 78 INTEGER(i_std) :: hist_id_sec,hist_id_sec2,hist_id_sto ,hist_id_stom_IPCC81 INTEGER(i_std) :: hist_id_sec,hist_id_sec2,hist_id_stom,hist_id_stom_IPCC 79 82 CHARACTER(LEN=30) :: time_str 80 REAL :: hist_days_stom,hist_dt_stom 83 REAL(r_std) :: dt_slow_ 84 REAL :: hist_days_stom,hist_days_stom_ipcc,hist_dt_stom,hist_dt_stom_ipcc 81 85 REAL(r_std),DIMENSION(10) :: hist_pool_10axis 82 86 REAL(r_std),DIMENSION(100) :: hist_pool_100axis 83 87 REAL(r_std),DIMENSION(11) :: hist_pool_11axis 84 88 REAL(r_std),DIMENSION(101) :: hist_pool_101axis 85 INTEGER :: hist_PFTaxis_id,h ori_id89 INTEGER :: hist_PFTaxis_id,hist_IPCC_PFTaxis_id,hori_id 86 90 INTEGER :: hist_pool_10axis_id 87 91 INTEGER :: hist_pool_100axis_id 88 92 INTEGER :: hist_pool_11axis_id 89 93 INTEGER :: hist_pool_101axis_id 90 INTEGER :: hist_level 91 INTEGER,PARAMETER :: max_hist_level = 10 92 INTEGER(i_std) :: i,j,iv,id 93 CHARACTER*80 :: var_name 94 CHARACTER(LEN=40),DIMENSION(10) :: fluxop 94 INTEGER(i_std) :: i,j,iv 95 95 LOGICAL :: ldrestart_read,ldrestart_write 96 96 LOGICAL :: l1d … … 104 104 REAL(r_std),DIMENSION(1) :: xtmp 105 105 INTEGER :: nsfm,nsft 106 INTEGER :: iisf 107 INTEGER(i_std) :: max_totsize,totsize_1step 108 INTEGER(i_std),DIMENSION(0:2) :: ifirst,ilast 109 INTEGER(i_std) :: iblocks,nblocks 110 INTEGER,PARAMETER :: ndm = 10 111 INTEGER,DIMENSION(ndm) :: start,count 112 INTEGER :: ndim,v_id 113 INTEGER :: force_id 106 INTEGER :: iisf,iiisf 107 INTEGER(i_std) :: max_totsize,totsize_1step,totsize_tmp 108 109 INTEGER :: vid 114 110 CHARACTER(LEN=100) :: forcing_name 115 111 REAL :: x 116 REAL(r_std),DIMENSION(:),ALLOCATABLE :: x_indices 117 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: x_neighbours 112 118 113 REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: var_3d 114 119 115 !- 120 116 REAL(r_std) :: time_sec,time_step_sec … … 122 118 REAL(r_std),DIMENSION(1) :: r1d 123 119 REAL(r_std) :: julian,djulian 124 ! REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: soiltype 125 !- 126 ! the following variables contain the forcing data 127 !- 128 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: clay_fm 129 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: humrel_x_fm 130 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: litterhum_fm 131 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: t2m_fm 132 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: t2m_min_fm 133 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: temp_sol_fm 134 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: soiltemp_fm 135 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: soilhum_fm 136 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:) :: precip_fm 137 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: gpp_x_fm 138 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: veget_force_x_fm 139 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: veget_max_force_x_fm 140 REAL(r_std),ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: lai_force_x_fm 141 INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:) :: isf 142 LOGICAL,ALLOCATABLE,SAVE,DIMENSION(:) :: nf_written 143 INTEGER(i_std),ALLOCATABLE,SAVE,DIMENSION(:) :: nf_cumul 144 INTEGER(i_std) :: ji,jv 120 121 INTEGER(i_std) :: ji,jv,l 122 123 LOGICAL :: debug 124 145 125 !--------------------------------------------------------------------- 146 !- No parallelisation yet in teststomate ! 147 #ifdef CPP_PARA 148 CALL ipslerr (3,'teststomate', & 149 & 'You try to run testsomate compiled with parallelisation. (CPP_PARA key)', & 150 & 'But it wasn''t programmed yet and teststomate will stop.','You must compiled it without CPP_PARA key.') 151 #endif 152 126 127 CALL init_para(.FALSE.) 128 CALL init_timer 153 129 !- 154 130 ! calendar … … 156 132 CALL ioconf_calendar ('noleap') 157 133 CALL ioget_calendar (one_year,one_day) 158 !- 159 ! open STOMATE's forcing file to read some basic info 160 !- 161 forcing_name = 'stomate_forcing.nc' 162 CALL getin ('STOMATE_FORCING_NAME',forcing_name) 163 iret = NF90_OPEN (TRIM(forcing_name),NF90_NOWRITE,force_id) 164 IF (iret /= NF90_NOERR) THEN 165 CALL ipslerr (3,'teststomate', & 166 & 'Could not open file : ', & 167 & forcing_name,'(Do you have forget it ?)') 168 ENDIF 169 ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'dtradia',dtradia) 170 ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'dt_slow',dt_force) 171 ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'nsft',x) 172 nsft = NINT(x) 173 ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'kjpij',x) 174 kjpij = NINT(x) 175 ier = NF90_GET_ATT (force_id,NF90_GLOBAL,'kjpindex',x) 176 kjpindex = NINT(x) 177 CALL init_grid( kjpindex ) 178 !- 179 write(*,*) 'ATTENTION',dtradia,dt_force 134 135 IF (is_root_prc) THEN 136 !- 137 ! open STOMATE's forcing file to read some basic info 138 !- 139 forcing_name = 'stomate_forcing.nc' 140 CALL getin ('STOMATE_FORCING_NAME',forcing_name) 141 iret = NF90_OPEN (TRIM(forcing_name),NF90_NOWRITE,forcing_id) 142 IF (iret /= NF90_NOERR) THEN 143 CALL ipslerr (3,'teststomate', & 144 & 'Could not open file : ', & 145 & forcing_name,'(Do you have forget it ?)') 146 ENDIF 147 ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'dtradia',dtradia) 148 ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'dt_slow',dt_force) 149 ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'nsft',x) 150 nsft = NINT(x) 151 ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'kjpij',x) 152 kjpij = NINT(x) 153 ier = NF90_GET_ATT (forcing_id,NF90_GLOBAL,'kjpindex',x) 154 nbp_glo = NINT(x) 155 ENDIF 156 CALL bcast(dtradia) 157 CALL bcast(dt_force) 158 CALL bcast(nsft) 159 CALL bcast(nbp_glo) 160 !- 161 write(numout,*) 'ATTENTION',dtradia,dt_force 162 !- 163 ! read info about land points 164 !- 165 IF (is_root_prc) THEN 166 a_er=.FALSE. 167 ALLOCATE (indices_g(nbp_glo),stat=ier) 168 a_er = a_er .OR. (ier.NE.0) 169 IF (a_er) THEN 170 CALL ipslerr (3,'teststomate', & 171 & 'PROBLEM WITH ALLOCATION', & 172 & 'for local variables 1','') 173 ENDIF 174 ! 175 ALLOCATE (x_indices_g(nbp_glo),stat=ier) 176 a_er = a_er .OR. (ier.NE.0) 177 IF (a_er) THEN 178 CALL ipslerr (3,'teststomate', & 179 & 'PROBLEM WITH ALLOCATION', & 180 & 'for global variables 1','') 181 ENDIF 182 ier = NF90_INQ_VARID (forcing_id,'index',vid) 183 IF (ier .NE. 0) THEN 184 CALL ipslerr (3,'teststomate', & 185 & 'PROBLEM WITH ALLOCATION', & 186 & 'for global variables 1','') 187 ENDIF 188 ier = NF90_GET_VAR (forcing_id,vid,x_indices_g) 189 IF (iret /= NF90_NOERR) THEN 190 CALL ipslerr (3,'teststomate', & 191 & 'PROBLEM WITH variable "index" in file ', & 192 & forcing_name,'(check this file)') 193 ENDIF 194 indices_g(:) = NINT(x_indices_g(:)) 195 DEALLOCATE (x_indices_g) 196 ELSE 197 ALLOCATE (indices_g(0)) 198 ENDIF 199 !--------------------------------------------------------------------- 200 !- 201 ! set debug to have more information 202 !- 203 !Config Key = DEBUG_INFO 204 !Config Desc = Flag for debug information 205 !Config Def = n 206 !Config Help = This option allows to switch on the output of debug 207 !Config information without recompiling the code. 208 !- 209 debug = .FALSE. 210 CALL getin_p('DEBUG_INFO',debug) 211 ! 212 !Config Key = LONGPRINT 213 !Config Desc = ORCHIDEE will print more messages 214 !Config Def = n 215 !Config Help = This flag permits to print more debug messages in the run. 216 ! 217 long_print = .FALSE. 218 CALL getin_p('LONGPRINT',long_print) 180 219 181 220 !!$ >> DS added for the externalised version … … 183 222 184 223 ! 1. Read the number of PFTs 185 CALL getin('NVM',nvm) 224 ! 225 !Config Key = NVM 226 !Config Desc = number of PFTs 227 !Config if = ANYTIME 228 !Config Def = 13 229 !Config Help = The number of vegetation types define by the user 230 !Config Units = NONE 231 CALL getin_p('NVM',nvm) 232 186 233 ! 2. Allocate and read the pft parameters stomate specific 187 234 CALL pft_parameters_main 235 !!$ DS 236 237 !- 238 ! activate CO2, STOMATE, but not sechiba 239 !- 240 control%river_routing = .FALSE. 241 control%hydrol_cwrr = .FALSE. 242 control%ok_sechiba = .FALSE. 243 ! 244 control%stomate_watchout = .TRUE. 245 control%ok_co2 = .TRUE. 246 control%ok_stomate = .TRUE. 247 !- 248 ! is DGVM activated? 249 !- 250 control%ok_dgvm = .FALSE. 251 CALL getin_p('STOMATE_OK_DGVM',control%ok_dgvm) 252 WRITE(numout,*) 'LPJ is activated: ',control%ok_dgvm 253 254 !!$ >> DS now we search the values for the scalar parameters in some .def files 255 256 ! 07/2011 257 CALL activate_sub_models(control%ok_sechiba,control%river_routing,control%ok_stomate) 258 CALL veget_config 259 ! 188 260 CALL getin_stomate_pft_parameters 189 !!$ DS 190 191 !- 192 ! allocate arrays 193 !- 261 CALL getin_co2_parameters 262 CALL getin_stomate_parameters 263 IF (control%ok_dgvm) THEN 264 CALL getin_dgvm_parameters 265 ENDIF 266 !!$ >> DS 267 268 !- 269 ! restart files 270 !- 271 IF (is_root_prc) THEN 272 ! Sechiba's restart files 273 sec_restname_in = 'sechiba_start.nc' 274 CALL getin('SECHIBA_restart_in',sec_restname_in) 275 WRITE(numout,*) 'SECHIBA INPUT RESTART_FILE: ',TRIM(sec_restname_in) 276 IF ( TRIM(sec_restname_in) .EQ. 'NONE' ) THEN 277 STOP 'Need a restart file for Sechiba' 278 ENDIF 279 sec_restname_out = 'sechiba_rest_out.nc' 280 CALL getin('SECHIBA_rest_out',sec_restname_out) 281 WRITE(numout,*) 'SECHIBA OUTPUT RESTART_FILE: ',TRIM(sec_restname_out) 282 ! Stomate's restart files 283 sto_restname_in = 'stomate_start.nc' 284 CALL getin('STOMATE_RESTART_FILEIN',sto_restname_in) 285 WRITE(numout,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in) 286 sto_restname_out = 'stomate_rest_out.nc' 287 CALL getin('STOMATE_RESTART_FILEOUT',sto_restname_out) 288 WRITE(numout,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out) 289 290 !- 291 ! We need to know iim, jjm. 292 ! Get them from the restart files themselves. 293 !- 294 iret = NF90_OPEN (sec_restname_in,NF90_NOWRITE,ncfid) 295 IF (iret /= NF90_NOERR) THEN 296 CALL ipslerr (3,'teststomate', & 297 & 'Could not open file : ', & 298 & sec_restname_in,'(Do you have forget it ?)') 299 ENDIF 300 iret = NF90_INQUIRE_DIMENSION (ncfid,1,len=iim_g) 301 iret = NF90_INQUIRE_DIMENSION (ncfid,2,len=jjm_g) 302 iret = NF90_CLOSE (ncfid) 303 304 ENDIF 305 CALL bcast(iim_g) 306 CALL bcast(jjm_g) 307 ! 308 ! Parallelization : 309 ! 310 CALL init_data_para(iim_g,jjm_g,nbp_glo,indices_g) 311 kjpindex=nbp_loc 312 jjm=jj_nb 313 iim=iim_g 314 kjpij=iim*jjm 315 IF (debug) WRITE(numout,*) "Local grid : ",kjpindex,iim,jjm 316 !- 317 !- 318 ! read info about grids 319 !- 320 !- 321 llm=1 322 ALLOCATE(lev(llm)) 323 IF (is_root_prc) THEN 324 !- 325 ier = NF90_INQ_VARID (forcing_id,'lalo',vid) 326 ier = NF90_GET_VAR (forcing_id,vid,lalo_g) 327 !- 328 ALLOCATE (x_neighbours_g(nbp_glo,8),stat=ier) 329 ier = NF90_INQ_VARID (forcing_id,'neighbours',vid) 330 ier = NF90_GET_VAR (forcing_id,vid,x_neighbours_g) 331 neighbours_g(:,:) = NINT(x_neighbours_g(:,:)) 332 DEALLOCATE (x_neighbours_g) 333 !- 334 ier = NF90_INQ_VARID (forcing_id,'resolution',vid) 335 ier = NF90_GET_VAR (forcing_id,vid,resolution_g) 336 !- 337 ier = NF90_INQ_VARID (forcing_id,'contfrac',vid) 338 ier = NF90_GET_VAR (forcing_id,vid,contfrac_g) 339 340 lon_g(:,:) = 0.0 341 lat_g(:,:) = 0.0 342 lev(1) = 0.0 343 !- 344 CALL restini & 345 & (sec_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, & 346 & sec_restname_out, itau_dep, date0, dt, rest_id_sec) 347 !- 348 IF ( dt .NE. dtradia ) THEN 349 WRITE(numout,*) 'dt',dt 350 WRITE(numout,*) 'dtradia',dtradia 351 CALL ipslerr (3,'teststomate', & 352 & 'PROBLEM with time steps.', & 353 & sec_restname_in,'(dt .NE. dtradia)') 354 ENDIF 355 !- 356 CALL restini & 357 & (sto_restname_in, iim_g, jjm_g, lon_g, lat_g, llm, lev, & 358 & sto_restname_out, itau_dep, date0, dt, rest_id_sto) 359 !- 360 IF ( dt .NE. dtradia ) THEN 361 WRITE(numout,*) 'dt',dt 362 WRITE(numout,*) 'dtradia',dtradia 363 CALL ipslerr (3,'teststomate', & 364 & 'PROBLEM with time steps.', & 365 & sto_restname_in,'(dt .NE. dtradia)') 366 ENDIF 367 ENDIF 368 CALL bcast(rest_id_sec) 369 CALL bcast(rest_id_sto) 370 CALL bcast(date0) 371 CALL bcast(dt) 372 CALL bcast(lev) 373 !--- 374 !--- Create the index table 375 !--- 376 !--- This job return a LOCAL kindex 377 !--- 378 ALLOCATE (indices(kjpindex),stat=ier) 379 IF (debug .AND. is_root_prc) WRITE(numout,*) "indices_g = ",indices_g(1:nbp_glo) 380 CALL scatter(indices_g,indices) 381 indices(1:kjpindex)=indices(1:kjpindex)-(jj_begin-1)*iim_g 382 IF (debug) WRITE(numout,*) "indices = ",indices(1:kjpindex) 383 384 !--- 385 !--- initialize global grid 386 !--- 387 CALL init_grid( kjpindex ) 388 CALL grid_stuff (nbp_glo, iim_g, jjm_g, lon_g, lat_g, indices_g) 389 390 !--- 391 !--- initialize local grid 392 !--- 393 jlandindex = (((indices(1:kjpindex)-1)/iim) + 1) 394 if (debug) WRITE(numout,*) "jlandindex = ",jlandindex(1:kjpindex) 395 ilandindex = (indices(1:kjpindex) - (jlandindex(1:kjpindex)-1)*iim) 396 IF (debug) WRITE(numout,*) "ilandindex = ",ilandindex(1:kjpindex) 397 ALLOCATE(lon(iim,jjm)) 398 ALLOCATE(lat(iim,jjm)) 399 CALL scatter2D(lon_g,lon) 400 CALL scatter2D(lat_g,lat) 401 IF (debug) THEN 402 IF (is_root_prc) THEN 403 WRITE(numout,*) mpi_rank, lat_g 404 WRITE(numout,*) mpi_rank, lon_g 405 ENDIF 406 WRITE(numout,*) mpi_rank, lat 407 WRITE(numout,*) mpi_rank, lon 408 ENDIF 409 410 DO ji=1,kjpindex 411 412 j = jlandindex(ji) 413 i = ilandindex(ji) 414 415 !- Create the internal coordinate table 416 !- 417 lalo(ji,1) = lat(i,j) 418 lalo(ji,2) = lon(i,j) 419 ENDDO 420 CALL scatter(neighbours_g,neighbours) 421 CALL scatter(resolution_g,resolution) 422 CALL scatter(contfrac_g,contfrac) 423 CALL scatter(area_g,area) 424 !- 425 !- Check if we have by any change a rectilinear grid. This would allow us to 426 !- simplify the output files. 427 ! 428 rectilinear = .FALSE. 429 IF (is_root_prc) THEN 430 IF ( ALL(lon_g(:,:) == SPREAD(lon_g(:,1), 2, SIZE(lon_g,2))) .AND. & 431 & ALL(lat_g(:,:) == SPREAD(lat_g(1,:), 1, SIZE(lat_g,1))) ) THEN 432 rectilinear = .TRUE. 433 ENDIF 434 ENDIF 435 CALL bcast(rectilinear) 436 IF (rectilinear) THEN 437 ALLOCATE(lon_rect(iim),stat=ier) 438 IF (ier .NE. 0) THEN 439 WRITE (numout,*) ' error in lon_rect allocation. We stop. We need iim words = ',iim 440 STOP 'intersurf_history' 441 ENDIF 442 ALLOCATE(lat_rect(jjm),stat=ier) 443 IF (ier .NE. 0) THEN 444 WRITE (numout,*) ' error in lat_rect allocation. We stop. We need jjm words = ',jjm 445 STOP 'intersurf_history' 446 ENDIF 447 lon_rect(:) = lon(:,1) 448 lat_rect(:) = lat(1,:) 449 ENDIF 450 !- 451 ! allocate arrays 452 !- 453 ! 194 454 a_er = .FALSE. 195 455 !!$ >> DS added for the externalised version … … 200 460 a_er = a_er .OR. (ier.NE.0) 201 461 !!$ end add DS 202 ALLOCATE (indices(kjpindex),stat=ier)203 a_er = a_er .OR. (ier.NE.0)204 462 ALLOCATE (indexveg(kjpindex*nvm), stat=ier) 205 463 a_er = a_er .OR. (ier.NE.0) … … 252 510 ALLOCATE (co2_flux(kjpindex,nvm),stat=ier) 253 511 a_er = a_er .OR. (ier.NE.0) 254 IF ( a_er ) STOP 'PROBLEM WITH ALLOCATION' 512 ALLOCATE (fco2_lu(kjpindex),stat=ier) 513 a_er = a_er .OR. (ier.NE.0) 514 IF (a_er) THEN 515 CALL ipslerr (3,'teststomate', & 516 & 'PROBLEM WITH ALLOCATION', & 517 & 'for local variables 1','') 518 ENDIF 255 519 ! 256 520 ! prepare forcing 257 521 ! 258 522 max_totsize = 50 259 CALL getin ('STOMATE_FORCING_MEMSIZE',max_totsize)523 CALL getin_p ('STOMATE_FORCING_MEMSIZE',max_totsize) 260 524 max_totsize = max_totsize * 1000000 525 261 526 totsize_1step = SIZE(soiltype(:,3))*KIND(soiltype(:,3)) + & 262 SIZE(humrel_x)*KIND(humrel_x) + & 263 SIZE(litterhum)*KIND(litterhum) + & 264 SIZE(t2m)*KIND(t2m) + & 265 SIZE(t2m_min)*KIND(t2m_min) + & 266 SIZE(temp_sol)*KIND(temp_sol) + & 267 SIZE(soiltemp)*KIND(soiltemp) + & 268 SIZE(soilhum)*KIND(soilhum) + & 269 SIZE(precip_rain)*KIND(precip_rain) + & 270 SIZE(precip_snow)*KIND(precip_snow) + & 271 SIZE(gpp_x)*KIND(gpp_x) + & 272 SIZE(veget_force_x)*KIND(veget_force_x) + & 273 SIZE(veget_max_force_x)*KIND(veget_max_force_x) + & 274 SIZE(lai_force_x)*KIND(lai_force_x) 527 SIZE(humrel_x)*KIND(humrel_x) + & 528 SIZE(litterhum)*KIND(litterhum) + & 529 SIZE(t2m)*KIND(t2m) + & 530 SIZE(t2m_min)*KIND(t2m_min) + & 531 SIZE(temp_sol)*KIND(temp_sol) + & 532 SIZE(soiltemp)*KIND(soiltemp) + & 533 SIZE(soilhum)*KIND(soilhum) + & 534 SIZE(precip_rain)*KIND(precip_rain) + & 535 SIZE(precip_snow)*KIND(precip_snow) + & 536 SIZE(gpp_x)*KIND(gpp_x) + & 537 SIZE(veget_force_x)*KIND(veget_force_x) + & 538 SIZE(veget_max_force_x)*KIND(veget_max_force_x) + & 539 SIZE(lai_force_x)*KIND(lai_force_x) 540 CALL reduce_sum(totsize_1step,totsize_tmp) 541 CALL bcast(totsize_tmp) 542 totsize_1step=totsize_tmp 543 275 544 ! total number of forcing steps 276 nsft = INT(one_year/(dt_force/one_day)) 545 IF ( nsft .NE. INT(one_year/(dt_force/one_day)) ) THEN 546 CALL ipslerr (3,'teststomate', & 547 & 'stomate: error with total number of forcing steps', & 548 & 'nsft','teststomate computation different with forcing file value.') 549 ENDIF 277 550 ! number of forcing steps in memory 278 nsfm = MIN(nsft,MAX(1,NINT(FLOAT(max_totsize)/FLOAT(totsize_1step)))) 551 nsfm = MIN(nsft, & 552 & MAX(1,NINT( REAL(max_totsize,r_std) & 553 & /REAL(totsize_1step,r_std)))) 279 554 !- 280 555 WRITE(numout,*) 'Offline forcing of Stomate:' … … 282 557 WRITE(numout,*) ' Number of forcing steps in memory:',nsfm 283 558 !- 284 a_er = .FALSE. 285 ALLOCATE (clay_fm(kjpindex,nsfm),stat=ier) 286 a_er = a_er.OR.(ier.NE.0) 287 ALLOCATE (humrel_x_fm(kjpindex,nvm,nsfm),stat=ier) 288 a_er = a_er.OR.(ier.NE.0) 289 ALLOCATE (litterhum_fm(kjpindex,nsfm),stat=ier) 290 a_er = a_er.OR.(ier.NE.0) 291 ALLOCATE (t2m_fm(kjpindex,nsfm),stat=ier) 292 a_er = a_er.OR.(ier.NE.0) 293 ALLOCATE (t2m_min_fm(kjpindex,nsfm),stat=ier) 294 a_er = a_er.OR.(ier.NE.0) 295 ALLOCATE (temp_sol_fm(kjpindex,nsfm),stat=ier) 296 a_er = a_er.OR.(ier.NE.0) 297 ALLOCATE (soiltemp_fm(kjpindex,nbdl,nsfm),stat=ier) 298 a_er = a_er.OR.(ier.NE.0) 299 ALLOCATE (soilhum_fm(kjpindex,nbdl,nsfm),stat=ier) 300 a_er = a_er.OR.(ier.NE.0) 301 ALLOCATE (precip_fm(kjpindex,nsfm),stat=ier) 302 a_er = a_er.OR.(ier.NE.0) 303 ALLOCATE (gpp_x_fm(kjpindex,nvm,nsfm),stat=ier) 304 a_er = a_er.OR.(ier.NE.0) 305 ALLOCATE (veget_force_x_fm(kjpindex,nvm,nsfm),stat=ier) 306 a_er = a_er.OR.(ier.NE.0) 307 ALLOCATE (veget_max_force_x_fm(kjpindex,nvm,nsfm),stat=ier) 308 a_er = a_er.OR.(ier.NE.0) 309 ALLOCATE (lai_force_x_fm(kjpindex,nvm,nsfm),stat=ier) 310 a_er = a_er.OR.(ier.NE.0) 311 ALLOCATE (isf(nsfm),stat=ier) 312 a_er = a_er.OR.(ier.NE.0) 313 ALLOCATE (nf_written(nsft),stat=ier) 314 a_er = a_er.OR.(ier.NE.0) 315 ALLOCATE (nf_cumul(nsft),stat=ier) 316 a_er = a_er.OR.(ier.NE.0) 317 IF (a_er) THEN 318 STOP 'stomate: error in memory allocation for forcing data' 319 ENDIF 559 CALL init_forcing(kjpindex,nsfm,nsft) 560 !- 320 561 ! ensure that we read all new forcing states 321 562 iisf = nsfm … … 323 564 ! of the forcing states that will be in memory 324 565 isf(:) = (/ (i,i=1,nsfm) /) 325 !- 326 ! read info about grids 327 !- 328 contfrac(:) = 1.0 329 !- 330 ALLOCATE (x_indices(kjpindex),stat=ier) 331 ier = NF90_INQ_VARID (force_id,'index',v_id) 332 ier = NF90_GET_VAR (force_id,v_id,x_indices) 333 indices(:) = NINT(x_indices(:)) 334 DEALLOCATE (x_indices) 335 !- 336 DO ji=1,kjpindex 337 DO jv=1,nvm 338 indexveg((jv-1)*kjpindex+ji) = indices(ji)+(jv-1)*kjpij 339 ENDDO 340 ENDDO 341 !- 342 ier = NF90_INQ_VARID (force_id,'lalo',v_id) 343 ier = NF90_GET_VAR (force_id,v_id,lalo) 344 !- 345 ALLOCATE (x_neighbours(kjpindex,8),stat=ier) 346 ier = NF90_INQ_VARID (force_id,'neighbours',v_id) 347 ier = NF90_GET_VAR (force_id,v_id,x_neighbours) 348 neighbours(:,:) = NINT(x_neighbours(:,:)) 349 DEALLOCATE (x_neighbours) 350 !- 351 ier = NF90_INQ_VARID (force_id,'resolution',v_id) 352 ier = NF90_GET_VAR (force_id,v_id,resolution) 353 !- 354 ier = NF90_INQ_VARID (force_id,'contfrac',v_id) 355 ier = NF90_GET_VAR (force_id,v_id,contfrac) 356 !- 357 ! activate CO2, STOMATE, but not sechiba 358 !- 359 control%river_routing = .FALSE. 360 control%hydrol_cwrr = .FALSE. 361 control%ok_sechiba = .FALSE. 362 ! 363 control%stomate_watchout = .TRUE. 364 control%ok_co2 = .TRUE. 365 control%ok_stomate = .TRUE. 366 !!$ >> DS now we search the values for the scalar parameters in some .def files 367 ! 07/2011 368 CALL activate_sub_models(control%ok_sechiba,control%river_routing,control%ok_stomate) 369 CALL veget_config 370 ! 371 CALL getin_co2_parameters 372 CALL getin_stomate_parameters 373 !!$ >> DS 374 !- 375 ! is DGVM activated? 376 !- 377 control%ok_dgvm = .FALSE. 378 CALL getin('STOMATE_OK_DGVM',control%ok_dgvm) 379 WRITE(*,*) 'LPJ is activated: ',control%ok_dgvm 380 !!$ >> DS for the externalisation reading the scalar parameters when ok_dgvm is set to true 381 IF (control%ok_dgvm) THEN 382 CALL getin_dgvm_parameters 383 ENDIF 384 !!$ >> DS 385 !- 386 ! restart files 387 !- 388 ! Sechiba's restart files 389 sec_restname_in = 'sechiba_start.nc' 390 CALL getin('SECHIBA_restart_in',sec_restname_in) 391 WRITE(*,*) 'SECHIBA INPUT RESTART_FILE: ',TRIM(sec_restname_in) 392 IF ( TRIM(sec_restname_in) .EQ. 'NONE' ) THEN 393 STOP 'Need a restart file for Sechiba' 394 ENDIF 395 sec_restname_out = 'sechiba_restart.nc' 396 CALL getin('SECHIBA_rest_out',sec_restname_out) 397 WRITE(*,*) 'SECHIBA OUTPUT RESTART_FILE: ',TRIM(sec_restname_out) 398 ! Stomate's restart files 399 sto_restname_in = 'stomate_start.nc' 400 CALL getin ('STOMATE_RESTART_FILEIN',sto_restname_in) 401 WRITE(*,*) 'STOMATE INPUT RESTART_FILE: ',TRIM(sto_restname_in) 402 sto_restname_out = 'stomate_restart.nc' 403 CALL getin ('STOMATE_RESTART_FILEOUT',sto_restname_out) 404 WRITE(*,*) 'STOMATE OUTPUT RESTART_FILE: ',TRIM(sto_restname_out) 405 !- 406 ! We need to know iim, jjm. 407 ! Get them from the restart files themselves. 408 !- 409 iret = NF90_OPEN (sec_restname_in,NF90_NOWRITE,ncfid) 410 IF (iret /= NF90_NOERR) THEN 411 CALL ipslerr (3,'teststomate', & 412 & 'Could not open file : ', & 413 & sec_restname_in,'(Do you have forget it ?)') 414 ENDIF 415 iret = NF90_INQUIRE_DIMENSION (ncfid,1,len=iim) 416 iret = NF90_INQUIRE_DIMENSION (ncfid,2,len=jjm) 417 iret = NF90_CLOSE (ncfid) 418 ! Allocate longitudes and latitudes 419 ALLOCATE (lon(iim,jjm),stat=ier) 420 a_er = a_er.OR.(ier.NE.0) 421 ALLOCATE (lat(iim,jjm),stat=ier) 422 a_er = a_er.OR.(ier.NE.0) 423 lon(:,:) = 0.0 424 lat(:,:) = 0.0 425 lev(1) = 0.0 426 !- 427 CALL restini & 428 & (sec_restname_in, iim, jjm, lon, lat, llm, lev, & 429 & sec_restname_out, itau_dep, date0, dt_files, rest_id_sec) 430 !- 431 CALL restini & 432 & (sto_restname_in, iim, jjm, lon, lat, llm, lev, & 433 & sto_restname_out, itau_dep, date0, dt_files, rest_id_sto) 434 !- 435 IF ( dt_files .NE. dtradia ) THEN 436 WRITE(*,*) 'dt_files',dt_files 437 WRITE(*,*) 'dtradia',dtradia 438 STOP 'PROBLEM with time steps.' 439 ENDIF 566 567 nf_written(:) = .FALSE. 568 nf_written(isf(:)) = .TRUE. 569 440 570 !- 441 571 ! a time step for STOMATE corresponds to itau_step timesteps in SECHIBA 442 572 !- 443 573 itau_step = NINT(dt_force/dtradia) 574 IF (debug) WRITE(numout,*) "dtradia, dt_rest, dt_force, itau_step",dtradia, dt, dt_force, itau_step 444 575 ! 445 576 CALL ioconf_startdate(date0) … … 451 582 CALL getin ('TIME_LENGTH', time_str) 452 583 ! transform into itau 453 CALL tlen2itau(time_str, dt _files, date0, itau_len)584 CALL tlen2itau(time_str, dt, date0, itau_len) 454 585 ! itau_len*dtradia must be a multiple of dt_force 455 586 itau_len = NINT( MAX(1.,FLOAT(NINT(itau_len*dtradia/dt_force))) & 456 & *dt_force/dtradia)457 !-458 ! set up STOMATE history file459 460 461 462 463 464 465 466 467 468 469 !-587 & *dt_force/dtradia) 588 !- 589 ! set up STOMATE history file 590 !- 591 !Config Key = STOMATE_OUTPUT_FILE 592 !Config Desc = Name of file in which STOMATE's output is going 593 !Config to be written 594 !Config Def = stomate_history.nc 595 !Config Help = This file is going to be created by the model 596 !Config and will contain the output from the model. 597 !Config This file is a truly COADS compliant netCDF file. 598 !Config It will be generated by the hist software from 599 !Config the IOIPSL package. 600 !- 470 601 stom_histname='stomate_history.nc' 471 CALL getin ('STOMATE_OUTPUT_FILE', stom_histname)472 WRITE( *,*) 'STOMATE_OUTPUT_FILE', TRIM(stom_histname)473 474 475 476 477 478 !-602 CALL getin_p ('STOMATE_OUTPUT_FILE', stom_histname) 603 WRITE(numout,*) 'STOMATE_OUTPUT_FILE', TRIM(stom_histname) 604 !- 605 !Config Key = STOMATE_HIST_DT 606 !Config Desc = STOMATE history time step (d) 607 !Config Def = 10. 608 !Config Help = Time step of the STOMATE history file 609 !- 479 610 hist_days_stom = 10. 480 CALL getin ('STOMATE_HIST_DT', hist_days_stom)611 CALL getin_p ('STOMATE_HIST_DT', hist_days_stom) 481 612 IF ( hist_days_stom == -1. ) THEN 482 613 hist_dt_stom = -1. … … 488 619 ENDIF 489 620 !- 490 ! initialize 491 CALL histbeg(stom_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, & 492 & itau_dep, date0, dt_files, hori_id, hist_id_sto) 493 ! define PFT axis 621 !- 622 !- initialize 623 WRITE(numout,*) "before histbeg : ",date0,dt 624 IF ( rectilinear ) THEN 625 #ifdef CPP_PARA 626 CALL histbeg(stom_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, & 627 & itau_dep, date0, dt, hori_id, hist_id_stom, domain_id=orch_domain_id) 628 #else 629 CALL histbeg(stom_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, & 630 & itau_dep, date0, dt, hori_id, hist_id_stom) 631 #endif 632 ELSE 633 #ifdef CPP_PARA 634 CALL histbeg(stom_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, & 635 & itau_dep, date0, dt, hori_id, hist_id_stom, domain_id=orch_domain_id) 636 #else 637 CALL histbeg(stom_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, & 638 & itau_dep, date0, dt, hori_id, hist_id_stom) 639 #endif 640 ENDIF 641 !- define PFT axis 494 642 hist_PFTaxis = (/ ( REAL(i,r_std), i=1,nvm ) /) 495 !declare this axis496 CALL histvert (hist_id_sto , 'PFT', 'Plant functional type', &497 & '-', nvm, hist_PFTaxis, hist_PFTaxis_id)643 !- declare this axis 644 CALL histvert (hist_id_stom, 'PFT', 'Plant functional type', & 645 & '1', nvm, hist_PFTaxis, hist_PFTaxis_id) 498 646 !- define Pool_10 axis 499 647 hist_pool_10axis = (/ ( REAL(i,r_std), i=1,10 ) /) 500 648 !- declare this axis 501 CALL histvert (hist_id_sto, 'P10', 'Pool 10 years', & 502 & '-', 10, hist_pool_10axis, hist_pool_10axis_id) 649 CALL histvert (hist_id_stom, 'P10', 'Pool 10 years', & 650 & '1', 10, hist_pool_10axis, hist_pool_10axis_id) 651 503 652 !- define Pool_100 axis 504 653 hist_pool_100axis = (/ ( REAL(i,r_std), i=1,100 ) /) 505 654 !- declare this axis 506 CALL histvert (hist_id_sto, 'P100', 'Pool 100 years', & 507 & '-', 100, hist_pool_100axis, hist_pool_100axis_id) 655 CALL histvert (hist_id_stom, 'P100', 'Pool 100 years', & 656 & '1', 100, hist_pool_100axis, hist_pool_100axis_id) 657 508 658 !- define Pool_11 axis 509 659 hist_pool_11axis = (/ ( REAL(i,r_std), i=1,11 ) /) 510 660 !- declare this axis 511 CALL histvert (hist_id_sto, 'P11', 'Pool 10 years + 1', & 512 & '-', 11, hist_pool_11axis, hist_pool_11axis_id) 661 CALL histvert (hist_id_stom, 'P11', 'Pool 10 years + 1', & 662 & '1', 11, hist_pool_11axis, hist_pool_11axis_id) 663 513 664 !- define Pool_101 axis 514 665 hist_pool_101axis = (/ ( REAL(i,r_std), i=1,101 ) /) 515 666 !- declare this axis 516 CALL histvert (hist_id_sto, 'P101', 'Pool 100 years + 1', & 517 & '-', 101, hist_pool_101axis, hist_pool_101axis_id) 518 ! define STOMATE history file 519 CALL stom_define_history (hist_id_sto, nvm, iim, jjm, & 520 & dt_files, hist_dt_stom, hori_id, hist_PFTaxis_id, & 667 CALL histvert (hist_id_stom, 'P101', 'Pool 100 years + 1', & 668 & '1', 101, hist_pool_101axis, hist_pool_101axis_id) 669 670 !- define STOMATE history file 671 CALL stom_define_history (hist_id_stom, nvm, iim, jjm, & 672 & dt, hist_dt_stom, hori_id, hist_PFTaxis_id, & 521 673 & hist_pool_10axis_id, hist_pool_100axis_id, & 522 674 & hist_pool_11axis_id, hist_pool_101axis_id) 523 ! end definition 524 CALL histend(hist_id_sto) 675 676 !- end definition 677 CALL histend(hist_id_stom) 678 !- 679 !- 680 ! STOMATE IPCC OUTPUTS IS ACTIVATED 681 !- 682 !Config Key = STOMATE_IPCC_OUTPUT_FILE 683 !Config Desc = Name of file in which STOMATE's output is going 684 !Config to be written 685 !Config Def = stomate_ipcc_history.nc 686 !Config Help = This file is going to be created by the model 687 !Config and will contain the output from the model. 688 !Config This file is a truly COADS compliant netCDF file. 689 !Config It will be generated by the hist software from 690 !Config the IOIPSL package. 691 !- 692 stom_ipcc_histname='stomate_ipcc_history.nc' 693 CALL getin_p('STOMATE_IPCC_OUTPUT_FILE', stom_ipcc_histname) 694 WRITE(numout,*) 'STOMATE_IPCC_OUTPUT_FILE', TRIM(stom_ipcc_histname) 695 !- 696 !Config Key = STOMATE_IPCC_HIST_DT 697 !Config Desc = STOMATE IPCC history time step (d) 698 !Config Def = 0. 699 !Config Help = Time step of the STOMATE IPCC history file 700 !- 701 hist_days_stom_ipcc = zero 702 CALL getin_p('STOMATE_IPCC_HIST_DT', hist_days_stom_ipcc) 703 IF ( hist_days_stom_ipcc == moins_un ) THEN 704 hist_dt_stom_ipcc = moins_un 705 WRITE(numout,*) 'output frequency for STOMATE IPCC history file (d): one month.' 706 ELSE 707 hist_dt_stom_ipcc = NINT( hist_days_stom_ipcc ) * one_day 708 WRITE(numout,*) 'output frequency for STOMATE IPCC history file (d): ', & 709 hist_dt_stom_ipcc/one_day 710 ENDIF 711 712 ! test consistency between STOMATE_IPCC_HIST_DT and DT_SLOW parameters 713 dt_slow_ = one_day 714 CALL getin_p('DT_SLOW', dt_slow_) 715 IF ( hist_days_stom_ipcc > zero ) THEN 716 IF (dt_slow_ > hist_dt_stom_ipcc) THEN 717 WRITE(numout,*) "DT_SLOW = ",dt_slow_," , STOMATE_IPCC_HIST_DT = ",hist_dt_stom_ipcc 718 CALL ipslerr (3,'intsurf_history', & 719 & 'Problem with DT_SLOW > STOMATE_IPCC_HIST_DT','', & 720 & '(must be less or equal)') 721 ENDIF 722 ENDIF 723 724 IF ( hist_dt_stom_ipcc == 0 ) THEN 725 hist_id_stom_ipcc = -1 726 ELSE 727 !- 728 !- initialize 729 IF ( rectilinear ) THEN 730 #ifdef CPP_PARA 731 CALL histbeg(stom_ipcc_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, & 732 & itau_dep, date0, dt, hori_id, hist_id_stom_ipcc,domain_id=orch_domain_id) 733 #else 734 CALL histbeg(stom_ipcc_histname, iim, lon_rect, jjm, lat_rect, 1, iim, 1, jjm, & 735 & itau_dep, date0, dt, hori_id, hist_id_stom_ipcc) 736 #endif 737 ELSE 738 #ifdef CPP_PARA 739 CALL histbeg(stom_ipcc_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, & 740 & itau_dep, date0, dt, hori_id, hist_id_stom_ipcc,domain_id=orch_domain_id) 741 #else 742 CALL histbeg(stom_ipcc_histname, iim, lon, jjm, lat, 1, iim, 1, jjm, & 743 & itau_dep, date0, dt, hori_id, hist_id_stom_ipcc) 744 #endif 745 ENDIF 746 !- declare this axis 747 CALL histvert (hist_id_stom_IPCC, 'PFT', 'Plant functional type', & 748 & '1', nvm, hist_PFTaxis, hist_IPCC_PFTaxis_id) 749 750 !- define STOMATE history file 751 CALL stom_IPCC_define_history (hist_id_stom_IPCC, nvm, iim, jjm, & 752 & dt, hist_dt_stom_ipcc, hori_id, hist_IPCC_PFTaxis_id) 753 754 !- end definition 755 CALL histend(hist_id_stom_IPCC) 756 757 ENDIF 758 ! 759 CALL histwrite(hist_id_stom, 'Areas', itau_dep+itau_step, area, kjpindex, indices) 760 IF ( hist_id_stom_IPCC > 0 ) THEN 761 CALL histwrite(hist_id_stom_IPCC, 'Areas', itau_dep+itau_step, area, kjpindex, indices) 762 ENDIF 763 ! 525 764 hist_id_sec = -1 526 765 hist_id_sec2 = -1 527 hist_id_stom_IPCC = -1 528 !- 529 ! read some variables we need from SECHIBA's restart file 530 !- 766 !- 767 ! first call of slowproc to initialize variables 768 !- 769 itau = itau_dep 770 ! 771 DO ji=1,kjpindex 772 DO jv=1,nvm 773 indexveg((jv-1)*kjpindex + ji) = indices(ji) + (jv-1)*kjpij 774 ENDDO 775 ENDDO 776 !- 777 !MM Problem here with dpu which depends on soil type 778 DO l = 1, nbdl-1 779 ! first 2.0 is dpu 780 ! second 2.0 is average 781 diaglev(l) = dpu_cste/(2**(nbdl-1) -1) * ( ( 2**(l-1) -1) + ( 2**(l) -1) ) / 2.0 782 ENDDO 783 diaglev(nbdl) = dpu_cste 784 ! 785 786 ! >> DS : getin the value of slowproc parameters 787 ! 788 !Config Key = CLAYFRACTION_DEFAULT 789 !Config Desc = 790 !Config If = OK_SECHIBA 791 !Config Def = 0.2 792 !Config Help = 793 !Config Units = NONE 794 CALL getin_p('CLAYFRACTION_DEFAULT',clayfraction_default) 795 ! 796 !Config Key = MIN_VEGFRAC 797 !Config Desc = Minimal fraction of mesh a vegetation type can occupy 798 !Config If = OK_SECHIBA 799 !Config Def = 0.001 800 !Config Help = 801 !Config Units = NONE 802 CALL getin_p('MIN_VEGFRAC',min_vegfrac) 803 ! 804 !Config Key = STEMPDIAG_BID 805 !Config Desc = only needed for an initial LAI if there is no restart file 806 !Config If = OK_SECHIBA 807 !Config Def = 280. 808 !Config Help = 809 !Config Units = 810 CALL getin_p('STEMPDIAG_BID',stempdiag_bid) 811 ! 812 !Config Key = SOILTYPE_DEFAULT 813 !Config Desc = Default soil texture distribution in the following order : sand, loam and clay 814 !Config If = OK_SECHIBA 815 !Config Def = 0.0, 1.0, 0.0 816 !Config Help = 817 !Config Units = NONE 818 CALL getin('SOILTYPE_DEFAULT',soiltype_default) 819 ! >> DS 820 531 821 CALL ioget_expval(val_exp) 532 !-533 ! first call of slowproc to initialize variables534 !-535 itau = itau_dep536 822 ldrestart_read = .TRUE. 537 823 ldrestart_write = .FALSE. 538 !- 539 !MM Problem here with dpu which depends on soil type 540 DO jv = 1, nbdl-1 541 ! first 2.0 is dpu 542 ! second 2.0 is average 543 diaglev(jv) = 2.0/(2**(nbdl-1) -1) * ( ( 2**(jv-1) -1) + ( 2**(jv) -1) ) / 2.0 544 ENDDO 545 diaglev(nbdl) = 2.0 546 !- 547 ! For sequential use only, we must initialize data_para : 548 ! 549 CALL init_para(.FALSE.) 550 ! 551 CALL init_data_para(iim,jjm,kjpindex,indices) 552 ! 553 !- global index index_g is the index_l of root proc 554 IF (is_root_prc) index_g(:)=indices(1:kjpindex) 555 !- 556 557 ! >> DS : getin the value of slowproc parameters 558 CALL getin('CLAYFRACTION_DEFAULT',clayfraction_default) 559 CALL getin('SOILTYPE_DEFAULT',soiltype_default) 560 ! >> DS 561 824 !- 825 ! read some variables we need from SECHIBA's restart file 826 !- 562 827 CALL slowproc_main & 563 828 & (itau, kjpij, kjpindex, dt_force, date0, & … … 568 833 & deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, & 569 834 & frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, & 570 & rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_sto, hist_id_stom_IPCC, co2_flux, & 571 !!>> DS add 01/07/2011 correction due to some modifications in the trunk 572 & fco2_lu) 835 & rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu) 573 836 ! correct date 574 837 day_counter = one_day - dt_force … … 577 840 ! time loop 578 841 !- 842 IF (debug) check_time=.TRUE. 843 CALL intsurf_time( itau_dep+itau_step, date0, dtradia ) 844 l_first_intersurf=.FALSE. 845 ! 579 846 DO itau = itau_dep+itau_step,itau_dep+itau_len,itau_step 847 ! 580 848 !-- next forcing state 581 849 iisf = iisf+1 … … 587 855 !---- determine blocks of forcing states that are contiguous in memory 588 856 !----- 589 nblocks = 0 590 ifirst(:) = 1 591 ilast(:) = 1 592 !----- 593 DO iisf=1,nsfm 594 IF ( (nblocks .NE. 0) ) THEN 595 IF ( (isf(iisf) .EQ. isf(ilast(nblocks))+1) ) THEN 596 !-------- element is contiguous with last element found 597 ilast(nblocks) = iisf 598 ELSE 599 !-------- found first element of new block 600 nblocks = nblocks+1 601 IF (nblocks .GT. nsfm) THEN 602 ! IF (nblocks .GT. 2) THEN 603 STOP 'Problem in teststomate' 604 ENDIF 605 ifirst(nblocks) = iisf 606 ilast(nblocks) = iisf 607 ENDIF 608 ELSE 609 !-------- found first element of new block 610 nblocks = nblocks+1 611 IF (nblocks .GT. nsfm) THEN 612 ! IF (nblocks .GT. 2) THEN 613 STOP 'Problem in teststomate' 614 ENDIF 615 ifirst(nblocks) = iisf 616 ilast(nblocks) = iisf 617 ENDIF 618 ENDDO 619 !----- 620 DO iblocks=1,nblocks 621 IF (ifirst(iblocks) .NE. ilast(iblocks)) THEN 622 ndim = 2; 623 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 624 count(1:ndim) = SHAPE(clay_fm) 625 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 626 ier = NF90_INQ_VARID (force_id,'clay',v_id) 627 a_er = a_er.OR.(ier.NE.0) 628 ier = NF90_GET_VAR (force_id,v_id, & 629 & clay_fm(:,ifirst(iblocks):ilast(iblocks)), & 630 & start=start(1:ndim), count=count(1:ndim)) 631 a_er = a_er.OR.(ier.NE.0) 632 !--------- 633 ndim = 3; 634 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 635 count(1:ndim) = SHAPE(humrel_x_fm) 636 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 637 ier = NF90_INQ_VARID (force_id,'humrel',v_id) 638 a_er = a_er.OR.(ier.NE.0) 639 ier = NF90_GET_VAR (force_id,v_id, & 640 & humrel_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), & 641 & start=start(1:ndim), count=count(1:ndim)) 642 a_er = a_er.OR.(ier.NE.0) 643 !--------- 644 ndim = 2; 645 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 646 count(1:ndim) = SHAPE(litterhum_fm) 647 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 648 ier = NF90_INQ_VARID (force_id,'litterhum',v_id) 649 a_er = a_er.OR.(ier.NE.0) 650 ier = NF90_GET_VAR (force_id,v_id, & 651 & litterhum_fm(:,ifirst(iblocks):ilast(iblocks)), & 652 & start=start(1:ndim), count=count(1:ndim)) 653 a_er = a_er.OR.(ier.NE.0) 654 !--------- 655 ndim = 2; 656 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 657 count(1:ndim) = SHAPE(t2m_fm) 658 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 659 ier = NF90_INQ_VARID (force_id,'t2m',v_id) 660 a_er = a_er.OR.(ier.NE.0) 661 ier = NF90_GET_VAR (force_id,v_id, & 662 & t2m_fm(:,ifirst(iblocks):ilast(iblocks)), & 663 & start=start(1:ndim), count=count(1:ndim)) 664 a_er = a_er.OR.(ier.NE.0) 665 !--------- 666 ndim = 2; 667 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 668 count(1:ndim) = SHAPE(t2m_min_fm) 669 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 670 ier = NF90_INQ_VARID (force_id,'t2m_min',v_id) 671 a_er = a_er.OR.(ier.NE.0) 672 ier = NF90_GET_VAR (force_id,v_id, & 673 & t2m_min_fm(:,ifirst(iblocks):ilast(iblocks)), & 674 & start=start(1:ndim), count=count(1:ndim)) 675 a_er = a_er.OR.(ier.NE.0) 676 !--------- 677 ndim = 2; 678 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 679 count(1:ndim) = SHAPE(temp_sol_fm) 680 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 681 ier = NF90_INQ_VARID (force_id,'tsurf',v_id) 682 a_er = a_er.OR.(ier.NE.0) 683 ier = NF90_GET_VAR (force_id,v_id, & 684 & temp_sol_fm(:,ifirst(iblocks):ilast(iblocks)), & 685 & start=start(1:ndim), count=count(1:ndim)) 686 a_er = a_er.OR.(ier.NE.0) 687 !--------- 688 ndim = 3; 689 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 690 count(1:ndim) = SHAPE(soiltemp_fm) 691 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 692 ier = NF90_INQ_VARID (force_id,'tsoil',v_id) 693 a_er = a_er.OR.(ier.NE.0) 694 ier = NF90_GET_VAR (force_id,v_id, & 695 & soiltemp_fm(:,:,ifirst(iblocks):ilast(iblocks)), & 696 & start=start(1:ndim), count=count(1:ndim)) 697 a_er = a_er.OR.(ier.NE.0) 698 !--------- 699 ndim = 3; 700 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 701 count(1:ndim) = SHAPE(soilhum_fm) 702 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 703 ier = NF90_INQ_VARID (force_id,'soilhum',v_id) 704 a_er = a_er.OR.(ier.NE.0) 705 ier = NF90_GET_VAR (force_id,v_id, & 706 & soilhum_fm(:,:,ifirst(iblocks):ilast(iblocks)), & 707 & start=start(1:ndim), count=count(1:ndim)) 708 a_er = a_er.OR.(ier.NE.0) 709 !--------- 710 ndim = 2; 711 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 712 count(1:ndim) = SHAPE(precip_fm) 713 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 714 ier = NF90_INQ_VARID (force_id,'precip',v_id) 715 a_er = a_er.OR.(ier.NE.0) 716 ier = NF90_GET_VAR (force_id,v_id, & 717 & precip_fm(:,ifirst(iblocks):ilast(iblocks)), & 718 & start=start(1:ndim), count=count(1:ndim)) 719 a_er = a_er.OR.(ier.NE.0) 720 !--------- 721 ndim = 3; 722 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 723 count(1:ndim) = SHAPE(gpp_x_fm) 724 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 725 ier = NF90_INQ_VARID (force_id,'gpp',v_id) 726 a_er = a_er.OR.(ier.NE.0) 727 ier = NF90_GET_VAR (force_id,v_id, & 728 & gpp_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), & 729 & start=start(1:ndim), count=count(1:ndim)) 730 a_er = a_er.OR.(ier.NE.0) 731 !--------- 732 ndim = 3; 733 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 734 count(1:ndim) = SHAPE(veget_force_x_fm) 735 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 736 ier = NF90_INQ_VARID (force_id,'veget',v_id) 737 a_er = a_er.OR.(ier.NE.0) 738 ier = NF90_GET_VAR (force_id,v_id, & 739 & veget_force_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), & 740 & start=start(1:ndim), count=count(1:ndim)) 741 a_er = a_er.OR.(ier.NE.0) 742 !--------- 743 ndim = 3; 744 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 745 count(1:ndim) = SHAPE(veget_max_force_x_fm) 746 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 747 ier = NF90_INQ_VARID (force_id,'veget_max',v_id) 748 a_er = a_er.OR.(ier.NE.0) 749 ier = NF90_GET_VAR (force_id,v_id, & 750 & veget_max_force_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), & 751 & start=start(1:ndim), count=count(1:ndim)) 752 a_er = a_er.OR.(ier.NE.0) 753 !--------- 754 ndim = 3; 755 start(:) = 1; start(ndim) = isf(ifirst(iblocks)); 756 count(1:ndim) = SHAPE(lai_force_x_fm) 757 count(ndim) = isf(ilast(iblocks))-isf(ifirst(iblocks))+1 758 ier = NF90_INQ_VARID (force_id,'lai',v_id) 759 a_er = a_er.OR.(ier.NE.0) 760 ier = NF90_GET_VAR (force_id,v_id, & 761 & lai_force_x_fm(:,:,ifirst(iblocks):ilast(iblocks)), & 762 & start=start(1:ndim), count=count(1:ndim)) 763 a_er = a_er.OR.(ier.NE.0) 764 ENDIF 765 ENDDO 857 CALL forcing_read(forcing_id,nsfm) 858 859 !-------------------------- 860 766 861 !----- 767 862 !---- determine which forcing states must be read next time … … 769 864 isf(1) = isf(nsfm)+1 770 865 IF ( isf(1) .GT. nsft ) isf(1) = 1 771 DO iisf = 2, nsfm 772 isf(iisf) = isf(iisf-1)+1 773 IF ( isf(iisf) .GT. nsft ) isf(iisf) = 1 774 ENDDO 866 DO iiisf = 2, nsfm 867 isf(iiisf) = isf(iiisf-1)+1 868 IF ( isf(iiisf) .GT. nsft ) isf(iiisf) = 1 869 ENDDO 870 nf_written(isf(:)) = .TRUE. 775 871 !---- start again at first forcing state 776 iisf = 1 777 ENDIF 778 soiltype(:,3) = clay_fm(:,iisf) 779 humrel_x(:,:) = humrel_x_fm(:,:,iisf) 780 litterhum(:) = litterhum_fm(:,iisf) 781 t2m(:) = t2m_fm(:,iisf) 782 t2m_min(:) = t2m_min_fm(:,iisf) 783 temp_sol(:) = temp_sol_fm(:,iisf) 784 soiltemp(:,:) = soiltemp_fm(:,:,iisf) 785 soilhum(:,:) = soilhum_fm(:,:,iisf) 786 precip_rain(:) = precip_fm(:,iisf) 787 gpp_x(:,:) = gpp_x_fm(:,:,iisf) 788 veget_force_x(:,:) = veget_force_x_fm(:,:,iisf) 789 veget_max_force_x(:,:) = veget_max_force_x_fm(:,:,iisf) 790 lai_force_x(:,:) = lai_force_x_fm(:,:,iisf) 791 WHERE ( t2m(:) .LT. ZeroCelsius ) 792 precip_snow(:) = precip_rain(:) 793 precip_rain(:) = 0. 794 ELSEWHERE 795 precip_snow(:) = 0. 796 ENDWHERE 872 iisf = 1 873 ENDIF 874 ! Bug here ! soiltype(:,3) != clay 875 ! soiltype(:,3) = clay_fm(:,iisf) 876 humrel_x(:,:) = humrel_daily_fm(:,:,iisf) 877 litterhum(:) = litterhum_daily_fm(:,iisf) 878 t2m(:) = t2m_daily_fm(:,iisf) 879 t2m_min(:) = t2m_min_daily_fm(:,iisf) 880 temp_sol(:) = tsurf_daily_fm(:,iisf) 881 soiltemp(:,:) = tsoil_daily_fm(:,:,iisf) 882 soilhum(:,:) = soilhum_daily_fm(:,:,iisf) 883 precip_rain(:) = precip_fm(:,iisf) 884 gpp_x(:,:) = gpp_daily_fm(:,:,iisf) 885 veget_force_x(:,:) = veget_fm(:,:,iisf) 886 veget_max_force_x(:,:) = veget_max_fm(:,:,iisf) 887 lai_force_x(:,:) = lai_fm(:,:,iisf) 888 WHERE ( t2m(:) .LT. ZeroCelsius ) 889 precip_snow(:) = precip_rain(:) 890 precip_rain(:) = 0. 891 ELSEWHERE 892 precip_snow(:) = 0. 893 ENDWHERE 797 894 !--- 798 895 !-- scale GPP to new lai and veget_max 799 896 !--- 800 WHERE ( lai_x(:,:) .EQ. 0.0 ) gpp_x(:,:) = 0.0897 WHERE ( lai_x(:,:) .EQ. 0.0 ) gpp_x(:,:) = 0.0 801 898 !-- scale GPP to new LAI 802 WHERE (lai_force_x(:,:) .GT. 0.0 )803 gpp_x(:,:) = gpp_x(:,:)*atan(2.*lai_x(:,:)) &804 & / atan( 2.*MAX(lai_force_x(:,:),0.01))805 ENDWHERE899 WHERE (lai_force_x(:,:) .GT. 0.0 ) 900 gpp_x(:,:) = gpp_x(:,:)*ATAN(2.*lai_x(:,:)) & 901 & /ATAN( 2.*MAX(lai_force_x(:,:),0.01)) 902 ENDWHERE 806 903 !-- scale GPP to new veget_max 807 WHERE (veget_max_force_x(:,:) .GT. 0.0 )808 gpp_x(:,:) = gpp_x(:,:)*veget_max_x(:,:)/veget_max_force_x(:,:)809 ENDWHERE904 WHERE (veget_max_force_x(:,:) .GT. 0.0 ) 905 gpp_x(:,:) = gpp_x(:,:)*veget_max_x(:,:)/veget_max_force_x(:,:) 906 ENDWHERE 810 907 !--- 811 908 !-- number crunching 812 909 !--- 813 CALL intsurf_time( itau, date0, dtradia ) 814 ldrestart_read = .FALSE. 815 ldrestart_write = .FALSE. 816 CALL slowproc_main & 910 ldrestart_read = .FALSE. 911 ldrestart_write = .FALSE. 912 CALL slowproc_main & 817 913 & (itau, kjpij, kjpindex, dt_force, date0, & 818 914 & ldrestart_read, ldrestart_write, .FALSE., .TRUE., & … … 822 918 & deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, & 823 919 & frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, & 824 & rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_sto , hist_id_stom_IPCC, co2_flux, &825 !!>> DS add 01/07/2011 correction due to some modifications in the trunk 826 & fco2_lu)827 day_counter = one_day - dt_force920 & rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu) 921 day_counter = one_day - dt_force 922 ! 923 CALL intsurf_time( itau+itau_step, date0, dtradia ) 828 924 ENDDO ! end of the time loop 829 925 !- 830 itau = itau -itau_step926 itau = itau -itau_step 831 927 !- 832 928 ! write restart files 833 929 !- 930 IF (is_root_prc) THEN 834 931 ! first, read and write variables that are not managed otherwise 835 taboo_vars = & 836 & '$lat$ $lon$ $lev$ $veget_year$ '// & 837 & '$height$ $day_counter$ $veget$ $veget_max$ $frac_nobio$ '// & 838 & '$lai$ $soiltype_frac$ $clay_frac$ '// & 839 & '$nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$' 840 !- 841 CALL ioget_vname(rest_id_sec, nbvar, varnames) 842 !!$!- 843 !!$! read and write some special variables (1D or variables that we need) 844 !!$!- 845 !!$ var_name = 'day_counter' 846 !!$ CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) 847 !!$ CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) 848 !!$!- 849 !!$ var_name = 'dt_days' 850 !!$ CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) 851 !!$ CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) 852 !!$!- 853 !!$ var_name = 'date' 854 !!$ CALL restget (rest_id_sto, var_name, 1, 1, 1, itau_dep, .TRUE., xtmp) 855 !!$ CALL restput (rest_id_sto, var_name, 1, 1, 1, itau_dep, xtmp) 856 !- 857 DO iv = 1, nbvar 932 taboo_vars = & 933 & '$lat$ $lon$ $lev$ $veget_year$ '// & 934 & '$height$ $day_counter$ $veget$ $veget_max$ $frac_nobio$ '// & 935 & '$lai$ $soiltype_frac$ $clay_frac$ '// & 936 & '$nav_lon$ $nav_lat$ $nav_lev$ $time$ $time_steps$' 937 !- 938 CALL ioget_vname(rest_id_sec, nbvar, varnames) 939 !- 940 DO iv = 1, nbvar 858 941 !-- check if the variable is to be written here 859 IF (INDEX( taboo_vars,'$'//TRIM(varnames(iv))//'$') .EQ. 0 ) THEN942 IF (INDEX( taboo_vars,'$'//TRIM(varnames(iv))//'$') .EQ. 0 ) THEN 860 943 !---- get variable dimensions, especially 3rd dimension 861 CALL ioget_vdim &862 & (rest_id_sec, varnames(iv), varnbdim_max, varnbdim, vardims)863 l1d = ALL(vardims(1:varnbdim) .EQ. 1)864 ALLOCATE(var_3d(kjpindex,vardims(3)),stat=ier)865 IF (ier .NE. 0) STOP 'ALLOCATION PROBLEM'944 CALL ioget_vdim & 945 & (rest_id_sec, varnames(iv), varnbdim_max, varnbdim, vardims) 946 l1d = ALL(vardims(1:varnbdim) .EQ. 1) 947 ALLOCATE(var_3d(nbp_glo,vardims(3)),stat=ier) 948 IF (ier .NE. 0) STOP 'ALLOCATION PROBLEM' 866 949 !---- read it 867 IF (l1d) THEN868 CALL restget (rest_id_sec,TRIM(varnames(iv)), &869 870 ELSE871 CALL restget (rest_id_sec,TRIM(varnames(iv)), &872 kjpindex,vardims(3),1,itau_dep,.TRUE.,var_3d, &873 "gather",kjpindex,indices)874 ENDIF950 IF (l1d) THEN 951 CALL restget (rest_id_sec,TRIM(varnames(iv)), & 952 1,vardims(3),1,itau_dep,.TRUE.,var_3d) 953 ELSE 954 CALL restget (rest_id_sec,TRIM(varnames(iv)), & 955 nbp_glo,vardims(3),1,itau_dep,.TRUE.,var_3d, & 956 "gather",nbp_glo,indices_g) 957 ENDIF 875 958 !---- write it 876 IF (l1d) THEN 877 CALL restput (rest_id_sec,TRIM(varnames(iv)), & 878 1,vardims(3),1,itau,var_3d) 879 ELSE 880 CALL restput (rest_id_sec,TRIM(varnames(iv)), & 881 kjpindex,vardims(3),1,itau,var_3d, & 882 'scatter',kjpindex,indices) 883 ENDIF 884 DEALLOCATE (var_3d) 885 ENDIF 886 ENDDO 959 IF (l1d) THEN 960 CALL restput (rest_id_sec,TRIM(varnames(iv)), & 961 1,vardims(3),1,itau,var_3d) 962 ELSE 963 CALL restput (rest_id_sec,TRIM(varnames(iv)), & 964 nbp_glo,vardims(3),1,itau,var_3d, & 965 'scatter',nbp_glo,indices_g) 966 ENDIF 967 DEALLOCATE (var_3d) 968 ENDIF 969 ENDDO 970 ENDIF 971 CALL barrier_para() 972 887 973 ! call slowproc to write restart files 888 974 ldrestart_read = .FALSE. … … 897 983 & deadleaf_cover, assim_param_x, lai_x, height_x, veget_x, & 898 984 & frac_nobio, veget_max_x, totfrac_nobio, qsintmax_x, & 899 & rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_sto, hist_id_stom_IPCC, co2_flux, & 900 !!>> DS add 01/07/2011 correction due to some modifications in the trunk 901 & fco2_lu) 902 985 & rest_id_sec, hist_id_sec, hist_id_sec2, rest_id_sto, hist_id_stom, hist_id_stom_IPCC, co2_flux, fco2_lu) 903 986 !- 904 987 ! close files 905 988 !- 906 CALL restclo 989 IF (is_root_prc) & 990 CALL restclo 907 991 CALL histclo 908 ier = NF90_CLOSE (force_id) 909 !- 910 ! write a new driver restart file with correct time step 911 !- 912 write(*,*) 'teststomate: writing driver restart file with correct time step.' 913 dri_restname_in = 'driver_start.nc' 914 CALL getin ('RESTART_FILEIN',dri_restname_in) 915 dri_restname_out = 'driver_restart.nc' 916 CALL getin ('RESTART_FILEOUT',dri_restname_out) 917 CALL SYSTEM & 918 & ('cp '//TRIM(dri_restname_in)//' '//TRIM(dri_restname_out)) 919 !- 920 iret = NF90_OPEN (TRIM(sec_restname_out),NF90_NOWRITE,ncfid) 921 iret = NF90_INQ_VARID (ncfid,'time',v_id) 922 iret = NF90_GET_VAR (ncfid,v_id,r1d) 923 time_sec = r1d(1) 924 iret = NF90_INQ_VARID (ncfid,'time_steps',v_id) 925 iret = NF90_GET_VAR (ncfid,v_id,time_step_sec) 926 iret = NF90_CLOSE (ncfid) 927 !- 928 iret = NF90_OPEN (TRIM(dri_restname_out),NF90_WRITE,ncfid) 929 iret = NF90_INQ_VARID (ncfid,'time',v_id) 930 iret = NF90_GET_VAR (ncfid,v_id,r1d) 931 time_dri = r1d(1) 932 r1d(1) = time_sec 933 iret = NF90_PUT_VAR (ncfid,v_id,r1d) 934 iret = NF90_INQ_VARID (ncfid,'time_steps',v_id) 935 iret = NF90_GET_VAR (ncfid,v_id,time_step_dri) 936 iret = NF90_PUT_VAR (ncfid,v_id,time_step_sec) 937 iret = NF90_INQ_VARID (ncfid,'julian',v_id) 938 iret = NF90_GET_VAR (ncfid,v_id,r1d) 939 julian = r1d(1) 940 djulian = (time_step_sec-time_step_dri)*dtradia/one_day 941 julian = julian & 942 & +djulian-FLOAT(INT((julian+djulian)/one_year))*one_year 943 r1d(1) = julian 944 iret = NF90_PUT_VAR (ncfid,v_id,r1d) 945 iret = NF90_CLOSE (ncfid) 946 947 CALL getin_dump 992 993 IF (is_root_prc) & 994 ier = NF90_CLOSE (forcing_id) 995 996 IF (is_root_prc) THEN 997 CALL getin_dump() 998 IF ( debug ) WRITE(numout,*) 'REST CLOSED' 999 ENDIF 1000 #ifdef CPP_PARA 1001 CALL MPI_FINALIZE(ier) 1002 #endif 1003 WRITE(numout,*) "End of teststomate." 1004 948 1005 !--------------- 949 1006 END PROGRAM teststomate
Note: See TracChangeset
for help on using the changeset viewer.
