Changeset 3605 for LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1DUTILS.h
- Timestamp:
- Nov 21, 2019, 4:43:45 PM (4 years ago)
- Location:
- LMDZ6/branches/Ocean_skin
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
libf/phylmd/dyn1d/1DUTILS.h (modified) (9 diffs, 1 prop)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Ocean_skin
-
LMDZ6/branches/Ocean_skin/libf/phylmd/dyn1d/1DUTILS.h
-
Property
svn:keywords
set to
Id
r3316 r3605 2 2 3 3 ! 4 ! $Id : conf_unicol.F 1279 2010-08-04 17:20:56Z lahellec$4 ! $Id$ 5 5 ! 6 6 ! … … 540 540 CALL getin('nudging_w',nudging_w) 541 541 542 ! RELIQUE ANCIENS FORMAT. ECRASE PAR LE SUIVANT 542 543 !Config Key = nudging_q 543 544 !Config Desc = forcage ou non par nudging sur q 544 545 !Config Def = false 545 546 !Config Help = forcage ou non par nudging sur q 546 nudging_q =0 547 CALL getin('nudging_q',nudging_q) 547 nudging_qv =0 548 CALL getin('nudging_q',nudging_qv) 549 CALL getin('nudging_qv',nudging_qv) 550 551 p_nudging_u=11000. 552 p_nudging_v=11000. 553 p_nudging_t=11000. 554 p_nudging_qv=11000. 555 CALL getin('p_nudging_u',p_nudging_u) 556 CALL getin('p_nudging_v',p_nudging_v) 557 CALL getin('p_nudging_t',p_nudging_t) 558 CALL getin('p_nudging_qv',p_nudging_qv) 548 559 549 560 !Config Key = nudging_t … … 599 610 write(lunout,*)' nudging_v = ', nudging_v 600 611 write(lunout,*)' nudging_t = ', nudging_t 601 write(lunout,*)' nudging_q = ', nudging_q612 write(lunout,*)' nudging_qv = ', nudging_qv 602 613 IF (forcing_type .eq.40) THEN 603 614 write(lunout,*) '--- Forcing type GCSS Old --- with:' … … 814 825 character*80 abort_message 815 826 ! 816 INTEGER nb 817 SAVE nb 818 DATA nb / 0 / 827 INTEGER pass 819 828 820 829 CALL open_restartphy(fichnom) … … 828 837 ENDDO 829 838 830 831 832 833 834 835 839 ! modname = 'dyn1dredem' 840 ! ierr = NF_OPEN(fichnom, NF_WRITE, nid) 841 ! IF (ierr .NE. NF_NOERR) THEN 842 ! abort_message="Pb. d ouverture "//fichnom 843 ! CALL abort_gcm('Modele 1D',abort_message,1) 844 ! ENDIF 836 845 837 846 DO l=1,length … … 885 894 tab_cntrl(31) = FLOAT(itau_dyn + itaufin) 886 895 ! 887 CALL put_var("controle","Param. de controle Dyn1D",tab_cntrl) 896 DO pass=1,2 897 CALL put_var(pass,"controle","Param. de controle Dyn1D",tab_cntrl) 888 898 ! 889 899 890 900 ! Ecriture/extension de la coordonnee temps 891 901 892 nb = nb + 1893 902 894 903 ! Ecriture des champs 895 904 ! 896 CALL put_field( "plev","p interfaces sauf la nulle",plev)897 CALL put_field( "play","",play)898 CALL put_field( "phi","geopotentielle",phi)899 CALL put_field( "phis","geopotentiell de surface",phis)900 CALL put_field( "presnivs","",presnivs)901 CALL put_field( "ucov","",ucov)902 CALL put_field( "vcov","",vcov)903 CALL put_field( "temp","",temp)904 CALL put_field( "omega2","",omega2)905 CALL put_field(pass,"plev","p interfaces sauf la nulle",plev) 906 CALL put_field(pass,"play","",play) 907 CALL put_field(pass,"phi","geopotentielle",phi) 908 CALL put_field(pass,"phis","geopotentiell de surface",phis) 909 CALL put_field(pass,"presnivs","",presnivs) 910 CALL put_field(pass,"ucov","",ucov) 911 CALL put_field(pass,"vcov","",vcov) 912 CALL put_field(pass,"temp","",temp) 913 CALL put_field(pass,"omega2","",omega2) 905 914 906 915 Do iq=1,nqtot 907 CALL put_field( "q"//nmq(iq),"eau vap ou condens et traceurs", &916 CALL put_field(pass,"q"//nmq(iq),"eau vap ou condens et traceurs", & 908 917 & q(:,:,iq)) 909 918 EndDo 910 CALL close_restartphy 919 IF (pass==1) CALL enddef_restartphy 920 IF (pass==2) CALL close_restartphy 921 922 923 ENDDO 911 924 912 925 ! … … 1458 1471 1459 1472 !====================================================================== 1460 SUBROUTINE read_togacoare(fich_toga,nlev_toga,nt_toga &1461 & ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga &1462 & ,ht_toga,vt_toga,hq_toga,vq_toga)1463 implicit none1464 1465 !-------------------------------------------------------------------------1466 ! Read TOGA-COARE forcing data1467 !-------------------------------------------------------------------------1468 1469 integer nlev_toga,nt_toga1470 real ts_toga(nt_toga),plev_toga(nlev_toga,nt_toga)1471 real t_toga(nlev_toga,nt_toga),q_toga(nlev_toga,nt_toga)1472 real u_toga(nlev_toga,nt_toga),v_toga(nlev_toga,nt_toga)1473 real w_toga(nlev_toga,nt_toga)1474 real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)1475 real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)1476 character*80 fich_toga1477 1478 integer k,ip1479 real bid1480 1481 integer iy,im,id,ih1482 1483 real plev_min1484 1485 plev_min = 55. ! pas de tendance de vap. d eau au-dessus de 55 hPa1486 1487 open(21,file=trim(fich_toga),form='formatted')1488 read(21,'(a)')1489 do ip = 1, nt_toga1490 read(21,'(a)')1491 read(21,'(a)')1492 read(21,223) iy, im, id, ih, bid, ts_toga(ip), bid,bid,bid,bid1493 read(21,'(a)')1494 read(21,'(a)')1495 1496 do k = 1, nlev_toga1497 read(21,230) plev_toga(k,ip), t_toga(k,ip), q_toga(k,ip) &1498 & ,u_toga(k,ip), v_toga(k,ip), w_toga(k,ip) &1499 & ,ht_toga(k,ip), vt_toga(k,ip), hq_toga(k,ip), vq_toga(k,ip)1500 1501 ! conversion in SI units:1502 t_toga(k,ip)=t_toga(k,ip)+273.15 ! K1503 q_toga(k,ip)=q_toga(k,ip)*0.001 ! kg/kg1504 w_toga(k,ip)=w_toga(k,ip)*100./3600. ! Pa/s1505 ! no water vapour tendency above 55 hPa1506 if (plev_toga(k,ip) .lt. plev_min) then1507 q_toga(k,ip) = 0.1508 hq_toga(k,ip) = 0.1509 vq_toga(k,ip) =0.1510 endif1511 enddo1512 1513 ts_toga(ip)=ts_toga(ip)+273.15 ! K1514 enddo1515 close(21)1516 1517 223 format(4i3,6f8.2)1518 230 format(6f9.3,4e11.3)1519 1520 return1521 end1522 1523 !-------------------------------------------------------------------------1524 SUBROUTINE read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu)1525 implicit none1526 1527 !-------------------------------------------------------------------------1528 ! Read I.SANDU case forcing data1529 !-------------------------------------------------------------------------1530 1531 integer nlev_sandu,nt_sandu1532 real ts_sandu(nt_sandu)1533 character*80 fich_sandu1534 1535 integer ip1536 integer iy,im,id,ih1537 1538 real plev_min1539 1540 print*,'nlev_sandu',nlev_sandu1541 plev_min = 55000. ! pas de tendance de vap. d eau au-dessus de 55 hPa1542 1543 open(21,file=trim(fich_sandu),form='formatted')1544 read(21,'(a)')1545 do ip = 1, nt_sandu1546 read(21,'(a)')1547 read(21,'(a)')1548 read(21,223) iy, im, id, ih, ts_sandu(ip)1549 print *,'ts=',iy,im,id,ih,ip,ts_sandu(ip)1550 enddo1551 close(21)1552 1553 223 format(4i3,f8.2)1554 1555 return1556 end1557 1558 !=====================================================================1559 !-------------------------------------------------------------------------1560 SUBROUTINE read_astex(fich_astex,nlev_astex,nt_astex,div_astex, &1561 & ts_astex,ug_astex,vg_astex,ufa_astex,vfa_astex)1562 implicit none1563 1564 !-------------------------------------------------------------------------1565 ! Read Astex case forcing data1566 !-------------------------------------------------------------------------1567 1568 integer nlev_astex,nt_astex1569 real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)1570 real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)1571 character*80 fich_astex1572 1573 integer ip1574 integer iy,im,id,ih1575 1576 real plev_min1577 1578 print*,'nlev_astex',nlev_astex1579 plev_min = 55000. ! pas de tendance de vap. d eau au-dessus de 55 hPa1580 1581 open(21,file=trim(fich_astex),form='formatted')1582 read(21,'(a)')1583 read(21,'(a)')1584 do ip = 1, nt_astex1585 read(21,'(a)')1586 read(21,'(a)')1587 read(21,223) iy, im, id, ih, div_astex(ip),ts_astex(ip), &1588 &ug_astex(ip),vg_astex(ip),ufa_astex(ip),vfa_astex(ip)1589 ts_astex(ip)=ts_astex(ip)+273.151590 print *,'ts=',iy,im,id,ih,ip,div_astex(ip),ts_astex(ip), &1591 &ug_astex(ip),vg_astex(ip),ufa_astex(ip),vg_astex(ip)1592 enddo1593 close(21)1594 1595 223 format(4i3,e13.2,f7.2,f7.3,f7.2,f7.3,f7.2)1596 1597 return1598 end1599 !=====================================================================1600 subroutine read_twpice(fich_twpice,nlevel,ntime &1601 & ,T_srf,plev,T,q,u,v,omega &1602 & ,T_adv_h,T_adv_v,q_adv_h,q_adv_v)1603 1604 !program reading forcings of the TWP-ICE experiment1605 1606 ! use netcdf1607 1608 implicit none1609 1610 #include "netcdf.inc"1611 1612 integer ntime,nlevel1613 integer l,k1614 character*80 :: fich_twpice1615 real*8 time(ntime)1616 real*8 lat, lon, alt, phis1617 real*8 lev(nlevel)1618 real*8 plev(nlevel,ntime)1619 1620 real*8 T(nlevel,ntime)1621 real*8 q(nlevel,ntime),u(nlevel,ntime)1622 real*8 v(nlevel,ntime)1623 real*8 omega(nlevel,ntime), div(nlevel,ntime)1624 real*8 T_adv_h(nlevel,ntime)1625 real*8 T_adv_v(nlevel,ntime), q_adv_h(nlevel,ntime)1626 real*8 q_adv_v(nlevel,ntime)1627 real*8 s(nlevel,ntime), s_adv_h(nlevel,ntime)1628 real*8 s_adv_v(nlevel,ntime)1629 real*8 p_srf_aver(ntime), p_srf_center(ntime)1630 real*8 T_srf(ntime)1631 1632 integer nid, ierr1633 integer nbvar3d1634 parameter(nbvar3d=20)1635 integer var3didin(nbvar3d)1636 1637 ierr = NF_OPEN(fich_twpice,NF_NOWRITE,nid)1638 if (ierr.NE.NF_NOERR) then1639 write(*,*) 'ERROR: Pb opening forcings cdf file '1640 write(*,*) NF_STRERROR(ierr)1641 stop ""1642 endif1643 1644 ierr=NF_INQ_VARID(nid,"lat",var3didin(1))1645 if(ierr/=NF_NOERR) then1646 write(*,*) NF_STRERROR(ierr)1647 stop 'lat'1648 endif1649 1650 ierr=NF_INQ_VARID(nid,"lon",var3didin(2))1651 if(ierr/=NF_NOERR) then1652 write(*,*) NF_STRERROR(ierr)1653 stop 'lon'1654 endif1655 1656 ierr=NF_INQ_VARID(nid,"alt",var3didin(3))1657 if(ierr/=NF_NOERR) then1658 write(*,*) NF_STRERROR(ierr)1659 stop 'alt'1660 endif1661 1662 ierr=NF_INQ_VARID(nid,"phis",var3didin(4))1663 if(ierr/=NF_NOERR) then1664 write(*,*) NF_STRERROR(ierr)1665 stop 'phis'1666 endif1667 1668 ierr=NF_INQ_VARID(nid,"T",var3didin(5))1669 if(ierr/=NF_NOERR) then1670 write(*,*) NF_STRERROR(ierr)1671 stop 'T'1672 endif1673 1674 ierr=NF_INQ_VARID(nid,"q",var3didin(6))1675 if(ierr/=NF_NOERR) then1676 write(*,*) NF_STRERROR(ierr)1677 stop 'q'1678 endif1679 1680 ierr=NF_INQ_VARID(nid,"u",var3didin(7))1681 if(ierr/=NF_NOERR) then1682 write(*,*) NF_STRERROR(ierr)1683 stop 'u'1684 endif1685 1686 ierr=NF_INQ_VARID(nid,"v",var3didin(8))1687 if(ierr/=NF_NOERR) then1688 write(*,*) NF_STRERROR(ierr)1689 stop 'v'1690 endif1691 1692 ierr=NF_INQ_VARID(nid,"omega",var3didin(9))1693 if(ierr/=NF_NOERR) then1694 write(*,*) NF_STRERROR(ierr)1695 stop 'omega'1696 endif1697 1698 ierr=NF_INQ_VARID(nid,"div",var3didin(10))1699 if(ierr/=NF_NOERR) then1700 write(*,*) NF_STRERROR(ierr)1701 stop 'div'1702 endif1703 1704 ierr=NF_INQ_VARID(nid,"T_adv_h",var3didin(11))1705 if(ierr/=NF_NOERR) then1706 write(*,*) NF_STRERROR(ierr)1707 stop 'T_adv_h'1708 endif1709 1710 ierr=NF_INQ_VARID(nid,"T_adv_v",var3didin(12))1711 if(ierr/=NF_NOERR) then1712 write(*,*) NF_STRERROR(ierr)1713 stop 'T_adv_v'1714 endif1715 1716 ierr=NF_INQ_VARID(nid,"q_adv_h",var3didin(13))1717 if(ierr/=NF_NOERR) then1718 write(*,*) NF_STRERROR(ierr)1719 stop 'q_adv_h'1720 endif1721 1722 ierr=NF_INQ_VARID(nid,"q_adv_v",var3didin(14))1723 if(ierr/=NF_NOERR) then1724 write(*,*) NF_STRERROR(ierr)1725 stop 'q_adv_v'1726 endif1727 1728 ierr=NF_INQ_VARID(nid,"s",var3didin(15))1729 if(ierr/=NF_NOERR) then1730 write(*,*) NF_STRERROR(ierr)1731 stop 's'1732 endif1733 1734 ierr=NF_INQ_VARID(nid,"s_adv_h",var3didin(16))1735 if(ierr/=NF_NOERR) then1736 write(*,*) NF_STRERROR(ierr)1737 stop 's_adv_h'1738 endif1739 1740 ierr=NF_INQ_VARID(nid,"s_adv_v",var3didin(17))1741 if(ierr/=NF_NOERR) then1742 write(*,*) NF_STRERROR(ierr)1743 stop 's_adv_v'1744 endif1745 1746 ierr=NF_INQ_VARID(nid,"p_srf_aver",var3didin(18))1747 if(ierr/=NF_NOERR) then1748 write(*,*) NF_STRERROR(ierr)1749 stop 'p_srf_aver'1750 endif1751 1752 ierr=NF_INQ_VARID(nid,"p_srf_center",var3didin(19))1753 if(ierr/=NF_NOERR) then1754 write(*,*) NF_STRERROR(ierr)1755 stop 'p_srf_center'1756 endif1757 1758 ierr=NF_INQ_VARID(nid,"T_srf",var3didin(20))1759 if(ierr/=NF_NOERR) then1760 write(*,*) NF_STRERROR(ierr)1761 stop 'T_srf'1762 endif1763 1764 !dimensions lecture1765 call catchaxis(nid,ntime,nlevel,time,lev,ierr)1766 1767 !pressure1768 do l=1,ntime1769 do k=1,nlevel1770 plev(k,l)=lev(k)1771 enddo1772 enddo1773 1774 #ifdef NC_DOUBLE1775 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),lat)1776 #else1777 ierr = NF_GET_VAR_REAL(nid,var3didin(1),lat)1778 #endif1779 if(ierr/=NF_NOERR) then1780 write(*,*) NF_STRERROR(ierr)1781 stop "getvarup"1782 endif1783 ! write(*,*)'lecture lat ok',lat1784 1785 #ifdef NC_DOUBLE1786 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),lon)1787 #else1788 ierr = NF_GET_VAR_REAL(nid,var3didin(2),lon)1789 #endif1790 if(ierr/=NF_NOERR) then1791 write(*,*) NF_STRERROR(ierr)1792 stop "getvarup"1793 endif1794 ! write(*,*)'lecture lon ok',lon1795 1796 #ifdef NC_DOUBLE1797 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),alt)1798 #else1799 ierr = NF_GET_VAR_REAL(nid,var3didin(3),alt)1800 #endif1801 if(ierr/=NF_NOERR) then1802 write(*,*) NF_STRERROR(ierr)1803 stop "getvarup"1804 endif1805 ! write(*,*)'lecture alt ok',alt1806 1807 #ifdef NC_DOUBLE1808 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),phis)1809 #else1810 ierr = NF_GET_VAR_REAL(nid,var3didin(4),phis)1811 #endif1812 if(ierr/=NF_NOERR) then1813 write(*,*) NF_STRERROR(ierr)1814 stop "getvarup"1815 endif1816 ! write(*,*)'lecture phis ok',phis1817 1818 #ifdef NC_DOUBLE1819 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),T)1820 #else1821 ierr = NF_GET_VAR_REAL(nid,var3didin(5),T)1822 #endif1823 if(ierr/=NF_NOERR) then1824 write(*,*) NF_STRERROR(ierr)1825 stop "getvarup"1826 endif1827 ! write(*,*)'lecture T ok'1828 1829 #ifdef NC_DOUBLE1830 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),q)1831 #else1832 ierr = NF_GET_VAR_REAL(nid,var3didin(6),q)1833 #endif1834 if(ierr/=NF_NOERR) then1835 write(*,*) NF_STRERROR(ierr)1836 stop "getvarup"1837 endif1838 ! write(*,*)'lecture q ok'1839 !q in kg/kg1840 do l=1,ntime1841 do k=1,nlevel1842 q(k,l)=q(k,l)/1000.1843 enddo1844 enddo1845 #ifdef NC_DOUBLE1846 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),u)1847 #else1848 ierr = NF_GET_VAR_REAL(nid,var3didin(7),u)1849 #endif1850 if(ierr/=NF_NOERR) then1851 write(*,*) NF_STRERROR(ierr)1852 stop "getvarup"1853 endif1854 ! write(*,*)'lecture u ok'1855 1856 #ifdef NC_DOUBLE1857 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),v)1858 #else1859 ierr = NF_GET_VAR_REAL(nid,var3didin(8),v)1860 #endif1861 if(ierr/=NF_NOERR) then1862 write(*,*) NF_STRERROR(ierr)1863 stop "getvarup"1864 endif1865 ! write(*,*)'lecture v ok'1866 1867 #ifdef NC_DOUBLE1868 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),omega)1869 #else1870 ierr = NF_GET_VAR_REAL(nid,var3didin(9),omega)1871 #endif1872 if(ierr/=NF_NOERR) then1873 write(*,*) NF_STRERROR(ierr)1874 stop "getvarup"1875 endif1876 ! write(*,*)'lecture omega ok'1877 !omega in mb/hour1878 do l=1,ntime1879 do k=1,nlevel1880 omega(k,l)=omega(k,l)*100./3600.1881 enddo1882 enddo1883 1884 #ifdef NC_DOUBLE1885 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),div)1886 #else1887 ierr = NF_GET_VAR_REAL(nid,var3didin(10),div)1888 #endif1889 if(ierr/=NF_NOERR) then1890 write(*,*) NF_STRERROR(ierr)1891 stop "getvarup"1892 endif1893 ! write(*,*)'lecture div ok'1894 1895 #ifdef NC_DOUBLE1896 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),T_adv_h)1897 #else1898 ierr = NF_GET_VAR_REAL(nid,var3didin(11),T_adv_h)1899 #endif1900 if(ierr/=NF_NOERR) then1901 write(*,*) NF_STRERROR(ierr)1902 stop "getvarup"1903 endif1904 ! write(*,*)'lecture T_adv_h ok'1905 !T adv in K/s1906 do l=1,ntime1907 do k=1,nlevel1908 T_adv_h(k,l)=T_adv_h(k,l)/3600.1909 enddo1910 enddo1911 1912 1913 #ifdef NC_DOUBLE1914 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),T_adv_v)1915 #else1916 ierr = NF_GET_VAR_REAL(nid,var3didin(12),T_adv_v)1917 #endif1918 if(ierr/=NF_NOERR) then1919 write(*,*) NF_STRERROR(ierr)1920 stop "getvarup"1921 endif1922 ! write(*,*)'lecture T_adv_v ok'1923 !T adv in K/s1924 do l=1,ntime1925 do k=1,nlevel1926 T_adv_v(k,l)=T_adv_v(k,l)/3600.1927 enddo1928 enddo1929 1930 #ifdef NC_DOUBLE1931 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),q_adv_h)1932 #else1933 ierr = NF_GET_VAR_REAL(nid,var3didin(13),q_adv_h)1934 #endif1935 if(ierr/=NF_NOERR) then1936 write(*,*) NF_STRERROR(ierr)1937 stop "getvarup"1938 endif1939 ! write(*,*)'lecture q_adv_h ok'1940 !q adv in kg/kg/s1941 do l=1,ntime1942 do k=1,nlevel1943 q_adv_h(k,l)=q_adv_h(k,l)/1000./3600.1944 enddo1945 enddo1946 1947 1948 #ifdef NC_DOUBLE1949 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),q_adv_v)1950 #else1951 ierr = NF_GET_VAR_REAL(nid,var3didin(14),q_adv_v)1952 #endif1953 if(ierr/=NF_NOERR) then1954 write(*,*) NF_STRERROR(ierr)1955 stop "getvarup"1956 endif1957 ! write(*,*)'lecture q_adv_v ok'1958 !q adv in kg/kg/s1959 do l=1,ntime1960 do k=1,nlevel1961 q_adv_v(k,l)=q_adv_v(k,l)/1000./3600.1962 enddo1963 enddo1964 1965 1966 #ifdef NC_DOUBLE1967 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),s)1968 #else1969 ierr = NF_GET_VAR_REAL(nid,var3didin(15),s)1970 #endif1971 if(ierr/=NF_NOERR) then1972 write(*,*) NF_STRERROR(ierr)1973 stop "getvarup"1974 endif1975 1976 #ifdef NC_DOUBLE1977 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),s_adv_h)1978 #else1979 ierr = NF_GET_VAR_REAL(nid,var3didin(16),s_adv_h)1980 #endif1981 if(ierr/=NF_NOERR) then1982 write(*,*) NF_STRERROR(ierr)1983 stop "getvarup"1984 endif1985 1986 #ifdef NC_DOUBLE1987 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),s_adv_v)1988 #else1989 ierr = NF_GET_VAR_REAL(nid,var3didin(17),s_adv_v)1990 #endif1991 if(ierr/=NF_NOERR) then1992 write(*,*) NF_STRERROR(ierr)1993 stop "getvarup"1994 endif1995 1996 #ifdef NC_DOUBLE1997 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),p_srf_aver)1998 #else1999 ierr = NF_GET_VAR_REAL(nid,var3didin(18),p_srf_aver)2000 #endif2001 if(ierr/=NF_NOERR) then2002 write(*,*) NF_STRERROR(ierr)2003 stop "getvarup"2004 endif2005 2006 #ifdef NC_DOUBLE2007 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),p_srf_center)2008 #else2009 ierr = NF_GET_VAR_REAL(nid,var3didin(19),p_srf_center)2010 #endif2011 if(ierr/=NF_NOERR) then2012 write(*,*) NF_STRERROR(ierr)2013 stop "getvarup"2014 endif2015 2016 #ifdef NC_DOUBLE2017 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),T_srf)2018 #else2019 ierr = NF_GET_VAR_REAL(nid,var3didin(20),T_srf)2020 #endif2021 if(ierr/=NF_NOERR) then2022 write(*,*) NF_STRERROR(ierr)2023 stop "getvarup"2024 endif2025 ! write(*,*)'lecture T_srf ok', T_srf2026 2027 return2028 end subroutine read_twpice2029 !=====================================================================2030 subroutine catchaxis(nid,ttm,llm,time,lev,ierr)2031 2032 ! use netcdf2033 2034 implicit none2035 #include "netcdf.inc"2036 integer nid,ttm,llm2037 real*8 time(ttm)2038 real*8 lev(llm)2039 integer ierr2040 2041 integer timevar,levvar2042 integer timelen,levlen2043 integer timedimin,levdimin2044 2045 ! Control & lecture on dimensions2046 ! ===============================2047 ierr=NF_INQ_DIMID(nid,"time",timedimin)2048 ierr=NF_INQ_VARID(nid,"time",timevar)2049 if (ierr.NE.NF_NOERR) then2050 write(*,*) 'ERROR: Field <time> is missing'2051 stop ""2052 endif2053 ierr=NF_INQ_DIMLEN(nid,timedimin,timelen)2054 2055 ierr=NF_INQ_DIMID(nid,"lev",levdimin)2056 ierr=NF_INQ_VARID(nid,"lev",levvar)2057 if (ierr.NE.NF_NOERR) then2058 write(*,*) 'ERROR: Field <lev> is lacking'2059 stop ""2060 endif2061 ierr=NF_INQ_DIMLEN(nid,levdimin,levlen)2062 2063 if((timelen/=ttm).or.(levlen/=llm)) then2064 write(*,*) 'ERROR: Not the good lenght for axis'2065 write(*,*) 'longitude: ',timelen,ttm+12066 write(*,*) 'latitude: ',levlen,llm2067 stop ""2068 endif2069 2070 !#ifdef NC_DOUBLE2071 ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)2072 ierr = NF_GET_VAR_DOUBLE(nid,levvar,lev)2073 !#else2074 ! ierr = NF_GET_VAR_REAL(nid,timevar,time)2075 ! ierr = NF_GET_VAR_REAL(nid,levvar,lev)2076 !#endif2077 2078 return2079 end2080 !=====================================================================2081 2082 SUBROUTINE interp_sandu_vertical(play,nlev_sandu,plev_prof &2083 & ,t_prof,thl_prof,q_prof,u_prof,v_prof,w_prof &2084 & ,omega_prof,o3mmr_prof &2085 & ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod &2086 & ,omega_mod,o3mmr_mod,mxcalc)2087 2088 implicit none2089 2090 #include "dimensions.h"2091 2092 !-------------------------------------------------------------------------2093 ! Vertical interpolation of SANDUREF forcing data onto model levels2094 !-------------------------------------------------------------------------2095 2096 integer nlevmax2097 parameter (nlevmax=41)2098 integer nlev_sandu,mxcalc2099 ! real play(llm), plev_prof(nlevmax)2100 ! real t_prof(nlevmax),q_prof(nlevmax)2101 ! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)2102 ! real ht_prof(nlevmax),vt_prof(nlevmax)2103 ! real hq_prof(nlevmax),vq_prof(nlevmax)2104 2105 real play(llm), plev_prof(nlev_sandu)2106 real t_prof(nlev_sandu),thl_prof(nlev_sandu),q_prof(nlev_sandu)2107 real u_prof(nlev_sandu),v_prof(nlev_sandu), w_prof(nlev_sandu)2108 real omega_prof(nlev_sandu),o3mmr_prof(nlev_sandu)2109 2110 real t_mod(llm),thl_mod(llm),q_mod(llm)2111 real u_mod(llm),v_mod(llm), w_mod(llm)2112 real omega_mod(llm),o3mmr_mod(llm)2113 2114 integer l,k,k1,k22115 real frac,frac1,frac2,fact2116 2117 do l = 1, llm2118 2119 if (play(l).ge.plev_prof(nlev_sandu)) then2120 2121 mxcalc=l2122 k1=02123 k2=02124 2125 if (play(l).le.plev_prof(1)) then2126 2127 do k = 1, nlev_sandu-12128 if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then2129 k1=k2130 k2=k+12131 endif2132 enddo2133 2134 if (k1.eq.0 .or. k2.eq.0) then2135 write(*,*) 'PB! k1, k2 = ',k1,k22136 write(*,*) 'l,play(l) = ',l,play(l)/1002137 do k = 1, nlev_sandu-12138 write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/1002139 enddo2140 endif2141 2142 frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))2143 t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))2144 thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))2145 q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))2146 u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))2147 v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))2148 w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))2149 omega_mod(l)=omega_prof(k2)-frac*(omega_prof(k2)-omega_prof(k1))2150 o3mmr_mod(l)=o3mmr_prof(k2)-frac*(o3mmr_prof(k2)-o3mmr_prof(k1))2151 2152 else !play>plev_prof(1)2153 2154 k1=12155 k2=22156 frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))2157 frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))2158 t_mod(l)= frac1*t_prof(k1) - frac2*t_prof(k2)2159 thl_mod(l)= frac1*thl_prof(k1) - frac2*thl_prof(k2)2160 q_mod(l)= frac1*q_prof(k1) - frac2*q_prof(k2)2161 u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)2162 v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)2163 w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)2164 omega_mod(l)= frac1*omega_prof(k1) - frac2*omega_prof(k2)2165 o3mmr_mod(l)= frac1*o3mmr_prof(k1) - frac2*o3mmr_prof(k2)2166 2167 endif ! play.le.plev_prof(1)2168 2169 else ! above max altitude of forcing file2170 2171 !jyg2172 fact=20.*(plev_prof(nlev_sandu)-play(l))/plev_prof(nlev_sandu) !jyg2173 fact = max(fact,0.) !jyg2174 fact = exp(-fact) !jyg2175 t_mod(l)= t_prof(nlev_sandu) !jyg2176 thl_mod(l)= thl_prof(nlev_sandu) !jyg2177 q_mod(l)= q_prof(nlev_sandu)*fact !jyg2178 u_mod(l)= u_prof(nlev_sandu)*fact !jyg2179 v_mod(l)= v_prof(nlev_sandu)*fact !jyg2180 w_mod(l)= w_prof(nlev_sandu)*fact !jyg2181 omega_mod(l)= omega_prof(nlev_sandu)*fact !jyg2182 o3mmr_mod(l)= o3mmr_prof(nlev_sandu)*fact !jyg2183 2184 endif ! play2185 2186 enddo ! l2187 2188 do l = 1,llm2189 ! print *,'t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l) ',2190 ! $ l,t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l)2191 enddo2192 2193 return2194 end2195 !=====================================================================2196 SUBROUTINE interp_astex_vertical(play,nlev_astex,plev_prof &2197 & ,t_prof,thl_prof,qv_prof,ql_prof,qt_prof,u_prof,v_prof &2198 & ,w_prof,tke_prof,o3mmr_prof &2199 & ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod &2200 & ,tke_mod,o3mmr_mod,mxcalc)2201 2202 implicit none2203 2204 #include "dimensions.h"2205 2206 !-------------------------------------------------------------------------2207 ! Vertical interpolation of Astex forcing data onto model levels2208 !-------------------------------------------------------------------------2209 2210 integer nlevmax2211 parameter (nlevmax=41)2212 integer nlev_astex,mxcalc2213 ! real play(llm), plev_prof(nlevmax)2214 ! real t_prof(nlevmax),qv_prof(nlevmax)2215 ! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)2216 ! real ht_prof(nlevmax),vt_prof(nlevmax)2217 ! real hq_prof(nlevmax),vq_prof(nlevmax)2218 2219 real play(llm), plev_prof(nlev_astex)2220 real t_prof(nlev_astex),thl_prof(nlev_astex),qv_prof(nlev_astex)2221 real u_prof(nlev_astex),v_prof(nlev_astex), w_prof(nlev_astex)2222 real o3mmr_prof(nlev_astex),ql_prof(nlev_astex)2223 real qt_prof(nlev_astex),tke_prof(nlev_astex)2224 2225 real t_mod(llm),thl_mod(llm),qv_mod(llm)2226 real u_mod(llm),v_mod(llm), w_mod(llm),tke_mod(llm)2227 real o3mmr_mod(llm),ql_mod(llm),qt_mod(llm)2228 2229 integer l,k,k1,k22230 real frac,frac1,frac2,fact2231 2232 do l = 1, llm2233 2234 if (play(l).ge.plev_prof(nlev_astex)) then2235 2236 mxcalc=l2237 k1=02238 k2=02239 2240 if (play(l).le.plev_prof(1)) then2241 2242 do k = 1, nlev_astex-12243 if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then2244 k1=k2245 k2=k+12246 endif2247 enddo2248 2249 if (k1.eq.0 .or. k2.eq.0) then2250 write(*,*) 'PB! k1, k2 = ',k1,k22251 write(*,*) 'l,play(l) = ',l,play(l)/1002252 do k = 1, nlev_astex-12253 write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/1002254 enddo2255 endif2256 2257 frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))2258 t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))2259 thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))2260 qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1))2261 ql_mod(l)= ql_prof(k2) - frac*(ql_prof(k2)-ql_prof(k1))2262 qt_mod(l)= qt_prof(k2) - frac*(qt_prof(k2)-qt_prof(k1))2263 u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))2264 v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))2265 w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))2266 tke_mod(l)= tke_prof(k2) - frac*(tke_prof(k2)-tke_prof(k1))2267 o3mmr_mod(l)=o3mmr_prof(k2)-frac*(o3mmr_prof(k2)-o3mmr_prof(k1))2268 2269 else !play>plev_prof(1)2270 2271 k1=12272 k2=22273 frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))2274 frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))2275 t_mod(l)= frac1*t_prof(k1) - frac2*t_prof(k2)2276 thl_mod(l)= frac1*thl_prof(k1) - frac2*thl_prof(k2)2277 qv_mod(l)= frac1*qv_prof(k1) - frac2*qv_prof(k2)2278 ql_mod(l)= frac1*ql_prof(k1) - frac2*ql_prof(k2)2279 qt_mod(l)= frac1*qt_prof(k1) - frac2*qt_prof(k2)2280 u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)2281 v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)2282 w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)2283 tke_mod(l)= frac1*tke_prof(k1) - frac2*tke_prof(k2)2284 o3mmr_mod(l)= frac1*o3mmr_prof(k1) - frac2*o3mmr_prof(k2)2285 2286 endif ! play.le.plev_prof(1)2287 2288 else ! above max altitude of forcing file2289 2290 !jyg2291 fact=20.*(plev_prof(nlev_astex)-play(l))/plev_prof(nlev_astex) !jyg2292 fact = max(fact,0.) !jyg2293 fact = exp(-fact) !jyg2294 t_mod(l)= t_prof(nlev_astex) !jyg2295 thl_mod(l)= thl_prof(nlev_astex) !jyg2296 qv_mod(l)= qv_prof(nlev_astex)*fact !jyg2297 ql_mod(l)= ql_prof(nlev_astex)*fact !jyg2298 qt_mod(l)= qt_prof(nlev_astex)*fact !jyg2299 u_mod(l)= u_prof(nlev_astex)*fact !jyg2300 v_mod(l)= v_prof(nlev_astex)*fact !jyg2301 w_mod(l)= w_prof(nlev_astex)*fact !jyg2302 tke_mod(l)= tke_prof(nlev_astex)*fact !jyg2303 o3mmr_mod(l)= o3mmr_prof(nlev_astex)*fact !jyg2304 2305 endif ! play2306 2307 enddo ! l2308 2309 do l = 1,llm2310 ! print *,'t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l) ',2311 ! $ l,t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l)2312 enddo2313 2314 return2315 end2316 2317 !======================================================================2318 SUBROUTINE read_rico(fich_rico,nlev_rico,ps_rico,play &2319 & ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico &2320 & ,dth_dyn,dqh_dyn)2321 implicit none2322 2323 !-------------------------------------------------------------------------2324 ! Read RICO forcing data2325 !-------------------------------------------------------------------------2326 #include "dimensions.h"2327 2328 2329 integer nlev_rico2330 real ts_rico,ps_rico2331 real t_rico(llm),q_rico(llm)2332 real u_rico(llm),v_rico(llm)2333 real w_rico(llm)2334 real dth_dyn(llm)2335 real dqh_dyn(llm)2336 2337 2338 real play(llm),zlay(llm)2339 2340 2341 real prico(nlev_rico),zrico(nlev_rico)2342 2343 character*80 fich_rico2344 2345 integer k,l2346 2347 2348 print*,fich_rico2349 open(21,file=trim(fich_rico),form='formatted')2350 do k=1,llm2351 zlay(k)=0.2352 enddo2353 2354 read(21,*) ps_rico,ts_rico2355 prico(1)=ps_rico2356 zrico(1)=0.02357 do l=2,nlev_rico2358 read(21,*) k,prico(l),zrico(l)2359 enddo2360 close(21)2361 2362 do k=1,llm2363 do l=1,802364 if(prico(l)>play(k)) then2365 if(play(k)>prico(l+1)) then2366 zlay(k)=zrico(l)+(play(k)-prico(l)) * &2367 & (zrico(l+1)-zrico(l))/(prico(l+1)-prico(l))2368 else2369 zlay(k)=zrico(l)+(play(k)-prico(80))* &2370 & (zrico(81)-zrico(80))/(prico(81)-prico(80))2371 endif2372 endif2373 enddo2374 print*,k,zlay(k)2375 ! U2376 if(0 < zlay(k) .and. zlay(k) < 4000) then2377 u_rico(k)=-9.9 + (-1.9 + 9.9)*zlay(k)/40002378 elseif(4000 < zlay(k) .and. zlay(k) < 12000) then2379 u_rico(k)= -1.9 + (30.0 + 1.9) / &2380 & (12000 - 4000) * (zlay(k) - 4000)2381 elseif(12000 < zlay(k) .and. zlay(k) < 13000) then2382 u_rico(k)=30.02383 elseif(13000 < zlay(k) .and. zlay(k) < 20000) then2384 u_rico(k)=30.0 - (30.0) / &2385 & (20000 - 13000) * (zlay(k) - 13000)2386 else2387 u_rico(k)=0.02388 endif2389 2390 !Q_v2391 if(0 < zlay(k) .and. zlay(k) < 740) then2392 q_rico(k)=16.0 + (13.8 - 16.0) / (740) * zlay(k)2393 elseif(740 < zlay(k) .and. zlay(k) < 3260) then2394 q_rico(k)=13.8 + (2.4 - 13.8) / &2395 & (3260 - 740) * (zlay(k) - 740)2396 elseif(3260 < zlay(k) .and. zlay(k) < 4000) then2397 q_rico(k)=2.4 + (1.8 - 2.4) / &2398 & (4000 - 3260) * (zlay(k) - 3260)2399 elseif(4000 < zlay(k) .and. zlay(k) < 9000) then2400 q_rico(k)=1.8 + (0 - 1.8) / &2401 & (9000 - 4000) * (zlay(k) - 4000)2402 else2403 q_rico(k)=0.02404 endif2405 2406 !T2407 if(0 < zlay(k) .and. zlay(k) < 740) then2408 t_rico(k)=299.2 + (292.0 - 299.2) / (740) * zlay(k)2409 elseif(740 < zlay(k) .and. zlay(k) < 4000) then2410 t_rico(k)=292.0 + (278.0 - 292.0) / &2411 & (4000 - 740) * (zlay(k) - 740)2412 elseif(4000 < zlay(k) .and. zlay(k) < 15000) then2413 t_rico(k)=278.0 + (203.0 - 278.0) / &2414 & (15000 - 4000) * (zlay(k) - 4000)2415 elseif(15000 < zlay(k) .and. zlay(k) < 17500) then2416 t_rico(k)=203.0 + (194.0 - 203.0) / &2417 & (17500 - 15000)* (zlay(k) - 15000)2418 elseif(17500 < zlay(k) .and. zlay(k) < 20000) then2419 t_rico(k)=194.0 + (206.0 - 194.0) / &2420 & (20000 - 17500)* (zlay(k) - 17500)2421 elseif(20000 < zlay(k) .and. zlay(k) < 60000) then2422 t_rico(k)=206.0 + (270.0 - 206.0) / &2423 & (60000 - 20000)* (zlay(k) - 20000)2424 endif2425 2426 ! W2427 if(0 < zlay(k) .and. zlay(k) < 2260 ) then2428 w_rico(k)=- (0.005/2260) * zlay(k)2429 elseif(2260 < zlay(k) .and. zlay(k) < 4000 ) then2430 w_rico(k)=- 0.0052431 elseif(4000 < zlay(k) .and. zlay(k) < 5000 ) then2432 w_rico(k)=- 0.005 + (0.005/ (5000 - 4000)) * (zlay(k) - 4000)2433 else2434 w_rico(k)=0.02435 endif2436 2437 ! dThrz+dTsw0+dTlw02438 if(0 < zlay(k) .and. zlay(k) < 4000) then2439 dth_dyn(k)=- 2.51 / 86400 + (-2.18 + 2.51 )/ &2440 & (86400*4000) * zlay(k)2441 elseif(4000 < zlay(k) .and. zlay(k) < 5000) then2442 dth_dyn(k)=- 2.18 / 86400 + ( 2.18 ) / &2443 & (86400*(5000 - 4000)) * (zlay(k) - 4000)2444 else2445 dth_dyn(k)=0.02446 endif2447 ! dQhrz2448 if(0 < zlay(k) .and. zlay(k) < 3000) then2449 dqh_dyn(k)=-1.0 / 86400 + (0.345 + 1.0)/ &2450 & (86400*3000) * (zlay(k))2451 elseif(3000 < zlay(k) .and. zlay(k) < 4000) then2452 dqh_dyn(k)=0.345 / 864002453 elseif(4000 < zlay(k) .and. zlay(k) < 5000) then2454 dqh_dyn(k)=0.345 / 86400 + &2455 & (-0.345)/(86400 * (5000 - 4000)) * (zlay(k)-4000)2456 else2457 dqh_dyn(k)=0.02458 endif2459 2460 !? if(play(k)>6e4) then2461 !? ratqs0(1,k)=ratqsbas*(plev(1)-play(k))/(plev(1)-6e4)2462 !? elseif((play(k)>3e4).and.(play(k)<6e4)) then2463 !? ratqs0(1,k)=ratqsbas+(ratqshaut-ratqsbas)&2464 !? *(6e4-play(k))/(6e4-3e4)2465 !? else2466 !? ratqs0(1,k)=ratqshaut2467 !? endif2468 2469 enddo2470 2471 do k=1,llm2472 q_rico(k)=q_rico(k)/1e32473 dqh_dyn(k)=dqh_dyn(k)/1e32474 v_rico(k)=-3.82475 enddo2476 2477 return2478 end2479 2480 !======================================================================2481 SUBROUTINE interp_sandu_time(day,day1,annee_ref &2482 & ,year_ini_sandu,day_ini_sandu,nt_sandu,dt_sandu &2483 & ,nlev_sandu,ts_sandu,ts_prof)2484 implicit none2485 2486 !---------------------------------------------------------------------------------------2487 ! Time interpolation of a 2D field to the timestep corresponding to day2488 !2489 ! day: current julian day (e.g. 717538.2)2490 ! day1: first day of the simulation2491 ! nt_sandu: total nb of data in the forcing (e.g. 13 for Sanduref)2492 ! dt_sandu: total time interval (in sec) between 2 forcing data (e.g. 6h for Sanduref)2493 !---------------------------------------------------------------------------------------2494 ! inputs:2495 integer annee_ref2496 integer nt_sandu,nlev_sandu2497 integer year_ini_sandu2498 real day, day1,day_ini_sandu,dt_sandu2499 real ts_sandu(nt_sandu)2500 ! outputs:2501 real ts_prof2502 ! local:2503 integer it_sandu1, it_sandu22504 real timeit,time_sandu1,time_sandu2,frac2505 ! Check that initial day of the simulation consistent with SANDU period:2506 if (annee_ref.ne.2006 ) then2507 print*,'Pour SANDUREF, annee_ref doit etre 2006 '2508 print*,'Changer annee_ref dans run.def'2509 stop2510 endif2511 ! if (annee_ref.eq.2006 .and. day1.lt.day_ini_sandu) then2512 ! print*,'SANDUREF debute le 15 Juillet 2006 (jour julien=196)'2513 ! print*,'Changer dayref dans run.def'2514 ! stop2515 ! endif2516 2517 ! Determine timestep relative to the 1st day of TOGA-COARE:2518 ! timeit=(day-day1)*86400.2519 ! if (annee_ref.eq.1992) then2520 ! timeit=(day-day_ini_sandu)*86400.2521 ! else2522 ! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 19922523 ! endif2524 timeit=(day-day_ini_sandu)*864002525 2526 ! Determine the closest observation times:2527 it_sandu1=INT(timeit/dt_sandu)+12528 it_sandu2=it_sandu1 + 12529 time_sandu1=(it_sandu1-1)*dt_sandu2530 time_sandu2=(it_sandu2-1)*dt_sandu2531 print *,'timeit day day_ini_sandu',timeit,day,day_ini_sandu2532 print *,'it_sandu1,it_sandu2,time_sandu1,time_sandu2', &2533 & it_sandu1,it_sandu2,time_sandu1,time_sandu22534 2535 if (it_sandu1 .ge. nt_sandu) then2536 write(*,*) 'PB-stop: day, it_sandu1, it_sandu2, timeit: ' &2537 & ,day,it_sandu1,it_sandu2,timeit/86400.2538 stop2539 endif2540 2541 ! time interpolation:2542 frac=(time_sandu2-timeit)/(time_sandu2-time_sandu1)2543 frac=max(frac,0.0)2544 2545 ts_prof = ts_sandu(it_sandu2) &2546 & -frac*(ts_sandu(it_sandu2)-ts_sandu(it_sandu1))2547 2548 print*, &2549 &'day,annee_ref,day_ini_sandu,timeit,it_sandu1,it_sandu2,SST:', &2550 &day,annee_ref,day_ini_sandu,timeit/86400.,it_sandu1, &2551 &it_sandu2,ts_prof2552 2553 return2554 END2555 !=====================================================================2556 !-------------------------------------------------------------------------2557 SUBROUTINE read_armcu(fich_armcu,nlev_armcu,nt_armcu, &2558 & sens,flat,adv_theta,rad_theta,adv_qt)2559 implicit none2560 2561 !-------------------------------------------------------------------------2562 ! Read ARM_CU case forcing data2563 !-------------------------------------------------------------------------2564 2565 integer nlev_armcu,nt_armcu2566 real sens(nt_armcu),flat(nt_armcu)2567 real adv_theta(nt_armcu),rad_theta(nt_armcu),adv_qt(nt_armcu)2568 character*80 fich_armcu2569 2570 integer ip2571 2572 integer iy,im,id,ih,in2573 2574 print*,'nlev_armcu',nlev_armcu2575 2576 open(21,file=trim(fich_armcu),form='formatted')2577 read(21,'(a)')2578 do ip = 1, nt_armcu2579 read(21,'(a)')2580 read(21,'(a)')2581 read(21,223) iy, im, id, ih, in, sens(ip),flat(ip), &2582 & adv_theta(ip),rad_theta(ip),adv_qt(ip)2583 print *,'forcages=',iy,im,id,ih,in, sens(ip),flat(ip), &2584 & adv_theta(ip),rad_theta(ip),adv_qt(ip)2585 enddo2586 close(21)2587 2588 223 format(5i3,5f8.3)2589 2590 return2591 end2592 2593 !=====================================================================2594 SUBROUTINE interp_toga_vertical(play,nlev_toga,plev_prof &2595 & ,t_prof,q_prof,u_prof,v_prof,w_prof &2596 & ,ht_prof,vt_prof,hq_prof,vq_prof &2597 & ,t_mod,q_mod,u_mod,v_mod,w_mod &2598 & ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)2599 2600 implicit none2601 2602 #include "dimensions.h"2603 2604 !-------------------------------------------------------------------------2605 ! Vertical interpolation of TOGA-COARE forcing data onto model levels2606 !-------------------------------------------------------------------------2607 2608 integer nlevmax2609 parameter (nlevmax=41)2610 integer nlev_toga,mxcalc2611 ! real play(llm), plev_prof(nlevmax)2612 ! real t_prof(nlevmax),q_prof(nlevmax)2613 ! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)2614 ! real ht_prof(nlevmax),vt_prof(nlevmax)2615 ! real hq_prof(nlevmax),vq_prof(nlevmax)2616 2617 real play(llm), plev_prof(nlev_toga)2618 real t_prof(nlev_toga),q_prof(nlev_toga)2619 real u_prof(nlev_toga),v_prof(nlev_toga), w_prof(nlev_toga)2620 real ht_prof(nlev_toga),vt_prof(nlev_toga)2621 real hq_prof(nlev_toga),vq_prof(nlev_toga)2622 2623 real t_mod(llm),q_mod(llm)2624 real u_mod(llm),v_mod(llm), w_mod(llm)2625 real ht_mod(llm),vt_mod(llm)2626 real hq_mod(llm),vq_mod(llm)2627 2628 integer l,k,k1,k22629 real frac,frac1,frac2,fact2630 2631 do l = 1, llm2632 2633 if (play(l).ge.plev_prof(nlev_toga)) then2634 2635 mxcalc=l2636 k1=02637 k2=02638 2639 if (play(l).le.plev_prof(1)) then2640 2641 do k = 1, nlev_toga-12642 if (play(l).le.plev_prof(k).and. play(l).gt.plev_prof(k+1)) then2643 k1=k2644 k2=k+12645 endif2646 enddo2647 2648 if (k1.eq.0 .or. k2.eq.0) then2649 write(*,*) 'PB! k1, k2 = ',k1,k22650 write(*,*) 'l,play(l) = ',l,play(l)/1002651 do k = 1, nlev_toga-12652 write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/1002653 enddo2654 endif2655 2656 frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))2657 t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))2658 q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_prof(k1))2659 u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))2660 v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))2661 w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))2662 ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))2663 vt_mod(l)= vt_prof(k2) - frac*(vt_prof(k2)-vt_prof(k1))2664 hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))2665 vq_mod(l)= vq_prof(k2) - frac*(vq_prof(k2)-vq_prof(k1))2666 2667 else !play>plev_prof(1)2668 2669 k1=12670 k2=22671 frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))2672 frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))2673 t_mod(l)= frac1*t_prof(k1) - frac2*t_prof(k2)2674 q_mod(l)= frac1*q_prof(k1) - frac2*q_prof(k2)2675 u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)2676 v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)2677 w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)2678 ht_mod(l)= frac1*ht_prof(k1) - frac2*ht_prof(k2)2679 vt_mod(l)= frac1*vt_prof(k1) - frac2*vt_prof(k2)2680 hq_mod(l)= frac1*hq_prof(k1) - frac2*hq_prof(k2)2681 vq_mod(l)= frac1*vq_prof(k1) - frac2*vq_prof(k2)2682 2683 endif ! play.le.plev_prof(1)2684 2685 else ! above max altitude of forcing file2686 2687 !jyg2688 fact=20.*(plev_prof(nlev_toga)-play(l))/plev_prof(nlev_toga) !jyg2689 fact = max(fact,0.) !jyg2690 fact = exp(-fact) !jyg2691 t_mod(l)= t_prof(nlev_toga) !jyg2692 q_mod(l)= q_prof(nlev_toga)*fact !jyg2693 u_mod(l)= u_prof(nlev_toga)*fact !jyg2694 v_mod(l)= v_prof(nlev_toga)*fact !jyg2695 w_mod(l)= 0.0 !jyg2696 ht_mod(l)= ht_prof(nlev_toga) !jyg2697 vt_mod(l)= vt_prof(nlev_toga) !jyg2698 hq_mod(l)= hq_prof(nlev_toga)*fact !jyg2699 vq_mod(l)= vq_prof(nlev_toga)*fact !jyg2700 2701 endif ! play2702 2703 enddo ! l2704 2705 ! do l = 1,llm2706 ! print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',2707 ! $ l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)2708 ! enddo2709 2710 return2711 end2712 2713 !=====================================================================2714 SUBROUTINE interp_case_vertical(play,nlev_cas,plev_prof_cas &2715 & ,t_prof_cas,q_prof_cas,u_prof_cas,v_prof_cas,ug_prof_cas,vg_prof_cas,vitw_prof_cas &2716 & ,du_prof_cas,hu_prof_cas,vu_prof_cas,dv_prof_cas,hv_prof_cas,vv_prof_cas &2717 & ,dt_prof_cas,ht_prof_cas,vt_prof_cas,dtrad_prof_cas,dq_prof_cas,hq_prof_cas,vq_prof_cas &2718 & ,t_mod_cas,q_mod_cas,u_mod_cas,v_mod_cas,ug_mod_cas,vg_mod_cas,w_mod_cas &2719 & ,du_mod_cas,hu_mod_cas,vu_mod_cas,dv_mod_cas,hv_mod_cas,vv_mod_cas &2720 & ,dt_mod_cas,ht_mod_cas,vt_mod_cas,dtrad_mod_cas,dq_mod_cas,hq_mod_cas,vq_mod_cas,mxcalc)2721 2722 implicit none2723 2724 #include "dimensions.h"2725 2726 !-------------------------------------------------------------------------2727 ! Vertical interpolation of TOGA-COARE forcing data onto mod_casel levels2728 !-------------------------------------------------------------------------2729 2730 integer nlevmax2731 parameter (nlevmax=41)2732 integer nlev_cas,mxcalc2733 ! real play(llm), plev_prof(nlevmax)2734 ! real t_prof(nlevmax),q_prof(nlevmax)2735 ! real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)2736 ! real ht_prof(nlevmax),vt_prof(nlevmax)2737 ! real hq_prof(nlevmax),vq_prof(nlevmax)2738 2739 real play(llm), plev_prof_cas(nlev_cas)2740 real t_prof_cas(nlev_cas),q_prof_cas(nlev_cas)2741 real u_prof_cas(nlev_cas),v_prof_cas(nlev_cas)2742 real ug_prof_cas(nlev_cas),vg_prof_cas(nlev_cas), vitw_prof_cas(nlev_cas)2743 real du_prof_cas(nlev_cas),hu_prof_cas(nlev_cas),vu_prof_cas(nlev_cas)2744 real dv_prof_cas(nlev_cas),hv_prof_cas(nlev_cas),vv_prof_cas(nlev_cas)2745 real dt_prof_cas(nlev_cas),ht_prof_cas(nlev_cas),vt_prof_cas(nlev_cas),dtrad_prof_cas(nlev_cas)2746 real dq_prof_cas(nlev_cas),hq_prof_cas(nlev_cas),vq_prof_cas(nlev_cas)2747 2748 real t_mod_cas(llm),q_mod_cas(llm)2749 real u_mod_cas(llm),v_mod_cas(llm)2750 real ug_mod_cas(llm),vg_mod_cas(llm), w_mod_cas(llm)2751 real du_mod_cas(llm),hu_mod_cas(llm),vu_mod_cas(llm)2752 real dv_mod_cas(llm),hv_mod_cas(llm),vv_mod_cas(llm)2753 real dt_mod_cas(llm),ht_mod_cas(llm),vt_mod_cas(llm),dtrad_mod_cas(llm)2754 real dq_mod_cas(llm),hq_mod_cas(llm),vq_mod_cas(llm)2755 2756 integer l,k,k1,k22757 real frac,frac1,frac2,fact2758 2759 do l = 1, llm2760 2761 if (play(l).ge.plev_prof_cas(nlev_cas)) then2762 2763 mxcalc=l2764 k1=02765 k2=02766 2767 if (play(l).le.plev_prof_cas(1)) then2768 2769 do k = 1, nlev_cas-12770 if (play(l).le.plev_prof_cas(k).and. play(l).gt.plev_prof_cas(k+1)) then2771 k1=k2772 k2=k+12773 endif2774 enddo2775 2776 if (k1.eq.0 .or. k2.eq.0) then2777 write(*,*) 'PB! k1, k2 = ',k1,k22778 write(*,*) 'l,play(l) = ',l,play(l)/1002779 do k = 1, nlev_cas-12780 write(*,*) 'k,plev_prof_cas(k) = ',k,plev_prof_cas(k)/1002781 enddo2782 endif2783 2784 frac = (plev_prof_cas(k2)-play(l))/(plev_prof_cas(k2)-plev_prof_cas(k1))2785 t_mod_cas(l)= t_prof_cas(k2) - frac*(t_prof_cas(k2)-t_prof_cas(k1))2786 q_mod_cas(l)= q_prof_cas(k2) - frac*(q_prof_cas(k2)-q_prof_cas(k1))2787 u_mod_cas(l)= u_prof_cas(k2) - frac*(u_prof_cas(k2)-u_prof_cas(k1))2788 v_mod_cas(l)= v_prof_cas(k2) - frac*(v_prof_cas(k2)-v_prof_cas(k1))2789 ug_mod_cas(l)= ug_prof_cas(k2) - frac*(ug_prof_cas(k2)-ug_prof_cas(k1))2790 vg_mod_cas(l)= vg_prof_cas(k2) - frac*(vg_prof_cas(k2)-vg_prof_cas(k1))2791 w_mod_cas(l)= vitw_prof_cas(k2) - frac*(vitw_prof_cas(k2)-vitw_prof_cas(k1))2792 du_mod_cas(l)= du_prof_cas(k2) - frac*(du_prof_cas(k2)-du_prof_cas(k1))2793 hu_mod_cas(l)= hu_prof_cas(k2) - frac*(hu_prof_cas(k2)-hu_prof_cas(k1))2794 vu_mod_cas(l)= vu_prof_cas(k2) - frac*(vu_prof_cas(k2)-vu_prof_cas(k1))2795 dv_mod_cas(l)= dv_prof_cas(k2) - frac*(dv_prof_cas(k2)-dv_prof_cas(k1))2796 hv_mod_cas(l)= hv_prof_cas(k2) - frac*(hv_prof_cas(k2)-hv_prof_cas(k1))2797 vv_mod_cas(l)= vv_prof_cas(k2) - frac*(vv_prof_cas(k2)-vv_prof_cas(k1))2798 dt_mod_cas(l)= dt_prof_cas(k2) - frac*(dt_prof_cas(k2)-dt_prof_cas(k1))2799 ht_mod_cas(l)= ht_prof_cas(k2) - frac*(ht_prof_cas(k2)-ht_prof_cas(k1))2800 vt_mod_cas(l)= vt_prof_cas(k2) - frac*(vt_prof_cas(k2)-vt_prof_cas(k1))2801 dq_mod_cas(l)= dq_prof_cas(k2) - frac*(dq_prof_cas(k2)-dq_prof_cas(k1))2802 hq_mod_cas(l)= hq_prof_cas(k2) - frac*(hq_prof_cas(k2)-hq_prof_cas(k1))2803 vq_mod_cas(l)= vq_prof_cas(k2) - frac*(vq_prof_cas(k2)-vq_prof_cas(k1))2804 dtrad_mod_cas(l)= dtrad_prof_cas(k2) - frac*(dtrad_prof_cas(k2)-dtrad_prof_cas(k1))2805 2806 else !play>plev_prof_cas(1)2807 2808 k1=12809 k2=22810 frac1 = (play(l)-plev_prof_cas(k2))/(plev_prof_cas(k1)-plev_prof_cas(k2))2811 frac2 = (play(l)-plev_prof_cas(k1))/(plev_prof_cas(k1)-plev_prof_cas(k2))2812 t_mod_cas(l)= frac1*t_prof_cas(k1) - frac2*t_prof_cas(k2)2813 q_mod_cas(l)= frac1*q_prof_cas(k1) - frac2*q_prof_cas(k2)2814 u_mod_cas(l)= frac1*u_prof_cas(k1) - frac2*u_prof_cas(k2)2815 v_mod_cas(l)= frac1*v_prof_cas(k1) - frac2*v_prof_cas(k2)2816 ug_mod_cas(l)= frac1*ug_prof_cas(k1) - frac2*ug_prof_cas(k2)2817 vg_mod_cas(l)= frac1*vg_prof_cas(k1) - frac2*vg_prof_cas(k2)2818 w_mod_cas(l)= frac1*vitw_prof_cas(k1) - frac2*vitw_prof_cas(k2)2819 du_mod_cas(l)= frac1*du_prof_cas(k1) - frac2*du_prof_cas(k2)2820 hu_mod_cas(l)= frac1*hu_prof_cas(k1) - frac2*hu_prof_cas(k2)2821 vu_mod_cas(l)= frac1*vu_prof_cas(k1) - frac2*vu_prof_cas(k2)2822 dv_mod_cas(l)= frac1*dv_prof_cas(k1) - frac2*dv_prof_cas(k2)2823 hv_mod_cas(l)= frac1*hv_prof_cas(k1) - frac2*hv_prof_cas(k2)2824 vv_mod_cas(l)= frac1*vv_prof_cas(k1) - frac2*vv_prof_cas(k2)2825 dt_mod_cas(l)= frac1*dt_prof_cas(k1) - frac2*dt_prof_cas(k2)2826 ht_mod_cas(l)= frac1*ht_prof_cas(k1) - frac2*ht_prof_cas(k2)2827 vt_mod_cas(l)= frac1*vt_prof_cas(k1) - frac2*vt_prof_cas(k2)2828 dq_mod_cas(l)= frac1*dq_prof_cas(k1) - frac2*dq_prof_cas(k2)2829 hq_mod_cas(l)= frac1*hq_prof_cas(k1) - frac2*hq_prof_cas(k2)2830 vq_mod_cas(l)= frac1*vq_prof_cas(k1) - frac2*vq_prof_cas(k2)2831 dtrad_mod_cas(l)= frac1*dtrad_prof_cas(k1) - frac2*dtrad_prof_cas(k2)2832 2833 endif ! play.le.plev_prof_cas(1)2834 2835 else ! above max altitude of forcing file2836 2837 !jyg2838 fact=20.*(plev_prof_cas(nlev_cas)-play(l))/plev_prof_cas(nlev_cas) !jyg2839 fact = max(fact,0.) !jyg2840 fact = exp(-fact) !jyg2841 t_mod_cas(l)= t_prof_cas(nlev_cas) !jyg2842 q_mod_cas(l)= q_prof_cas(nlev_cas)*fact !jyg2843 u_mod_cas(l)= u_prof_cas(nlev_cas)*fact !jyg2844 v_mod_cas(l)= v_prof_cas(nlev_cas)*fact !jyg2845 ug_mod_cas(l)= ug_prof_cas(nlev_cas)*fact !jyg2846 vg_mod_cas(l)= vg_prof_cas(nlev_cas)*fact !jyg2847 w_mod_cas(l)= 0.0 !jyg2848 du_mod_cas(l)= du_prof_cas(nlev_cas)*fact2849 hu_mod_cas(l)= hu_prof_cas(nlev_cas)*fact !jyg2850 vu_mod_cas(l)= vu_prof_cas(nlev_cas)*fact !jyg2851 dv_mod_cas(l)= dv_prof_cas(nlev_cas)*fact2852 hv_mod_cas(l)= hv_prof_cas(nlev_cas)*fact !jyg2853 vv_mod_cas(l)= vv_prof_cas(nlev_cas)*fact !jyg2854 dt_mod_cas(l)= dt_prof_cas(nlev_cas)2855 ht_mod_cas(l)= ht_prof_cas(nlev_cas) !jyg2856 vt_mod_cas(l)= vt_prof_cas(nlev_cas) !jyg2857 dq_mod_cas(l)= dq_prof_cas(nlev_cas)*fact2858 hq_mod_cas(l)= hq_prof_cas(nlev_cas)*fact !jyg2859 vq_mod_cas(l)= vq_prof_cas(nlev_cas)*fact !jyg2860 dtrad_mod_cas(l)= dtrad_prof_cas(nlev_cas)*fact !jyg2861 2862 endif ! play2863 2864 enddo ! l2865 2866 ! do l = 1,llm2867 ! print *,'t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l) ',2868 ! $ l,t_mod_cas(l),q_mod_cas(l),ht_mod_cas(l),hq_mod_cas(l)2869 ! enddo2870 2871 return2872 end2873 !*****************************************************************************2874 !=====================================================================2875 SUBROUTINE interp_dice_vertical(play,nlev_dice,nt_dice,plev_prof &2876 & ,th_prof,qv_prof,u_prof,v_prof,o3_prof &2877 & ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof &2878 & ,th_mod,qv_mod,u_mod,v_mod,o3_mod &2879 & ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc)2880 2881 implicit none2882 2883 #include "dimensions.h"2884 2885 !-------------------------------------------------------------------------2886 ! Vertical interpolation of Dice forcing data onto model levels2887 !-------------------------------------------------------------------------2888 2889 integer nlevmax2890 parameter (nlevmax=41)2891 integer nlev_dice,mxcalc,nt_dice2892 2893 real play(llm), plev_prof(nlev_dice)2894 real th_prof(nlev_dice),qv_prof(nlev_dice)2895 real u_prof(nlev_dice),v_prof(nlev_dice)2896 real o3_prof(nlev_dice)2897 real ht_prof(nlev_dice),hq_prof(nlev_dice)2898 real hu_prof(nlev_dice),hv_prof(nlev_dice)2899 real w_prof(nlev_dice),omega_prof(nlev_dice)2900 2901 real th_mod(llm),qv_mod(llm)2902 real u_mod(llm),v_mod(llm), o3_mod(llm)2903 real ht_mod(llm),hq_mod(llm)2904 real hu_mod(llm),hv_mod(llm),w_mod(llm),omega_mod(llm)2905 2906 integer l,k,k1,k2,kp2907 real aa,frac,frac1,frac2,fact2908 2909 do l = 1, llm2910 2911 if (play(l).ge.plev_prof(nlev_dice)) then2912 2913 mxcalc=l2914 k1=02915 k2=02916 2917 if (play(l).le.plev_prof(1)) then2918 2919 do k = 1, nlev_dice-12920 if (play(l).le.plev_prof(k) .and. play(l).gt.plev_prof(k+1)) then2921 k1=k2922 k2=k+12923 endif2924 enddo2925 2926 if (k1.eq.0 .or. k2.eq.0) then2927 write(*,*) 'PB! k1, k2 = ',k1,k22928 write(*,*) 'l,play(l) = ',l,play(l)/1002929 do k = 1, nlev_dice-12930 write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/1002931 enddo2932 endif2933 2934 frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))2935 th_mod(l)= th_prof(k2) - frac*(th_prof(k2)-th_prof(k1))2936 qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1))2937 u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))2938 v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))2939 o3_mod(l)= o3_prof(k2) - frac*(o3_prof(k2)-o3_prof(k1))2940 ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))2941 hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))2942 hu_mod(l)= hu_prof(k2) - frac*(hu_prof(k2)-hu_prof(k1))2943 hv_mod(l)= hv_prof(k2) - frac*(hv_prof(k2)-hv_prof(k1))2944 w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))2945 omega_mod(l)= omega_prof(k2) - frac*(omega_prof(k2)-omega_prof(k1))2946 2947 else !play>plev_prof(1)2948 2949 k1=12950 k2=22951 frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))2952 frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))2953 th_mod(l)= frac1*th_prof(k1) - frac2*th_prof(k2)2954 qv_mod(l)= frac1*qv_prof(k1) - frac2*qv_prof(k2)2955 u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)2956 v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)2957 o3_mod(l)= frac1*o3_prof(k1) - frac2*o3_prof(k2)2958 ht_mod(l)= frac1*ht_prof(k1) - frac2*ht_prof(k2)2959 hq_mod(l)= frac1*hq_prof(k1) - frac2*hq_prof(k2)2960 hu_mod(l)= frac1*hu_prof(k1) - frac2*hu_prof(k2)2961 hv_mod(l)= frac1*hv_prof(k1) - frac2*hv_prof(k2)2962 w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)2963 omega_mod(l)= frac1*omega_prof(k1) - frac2*omega_prof(k2)2964 2965 endif ! play.le.plev_prof(1)2966 2967 else ! above max altitude of forcing file2968 2969 !jyg2970 fact=20.*(plev_prof(nlev_dice)-play(l))/plev_prof(nlev_dice) !jyg2971 fact = max(fact,0.) !jyg2972 fact = exp(-fact) !jyg2973 th_mod(l)= th_prof(nlev_dice) !jyg2974 qv_mod(l)= qv_prof(nlev_dice)*fact !jyg2975 u_mod(l)= u_prof(nlev_dice)*fact !jyg2976 v_mod(l)= v_prof(nlev_dice)*fact !jyg2977 o3_mod(l)= o3_prof(nlev_dice)*fact !jyg2978 ht_mod(l)= ht_prof(nlev_dice) !jyg2979 hq_mod(l)= hq_prof(nlev_dice)*fact !jyg2980 hu_mod(l)= hu_prof(nlev_dice) !jyg2981 hv_mod(l)= hv_prof(nlev_dice) !jyg2982 w_mod(l)= 0. !jyg2983 omega_mod(l)= 0. !jyg2984 2985 endif ! play2986 2987 enddo ! l2988 2989 ! do l = 1,llm2990 ! print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',2991 ! $ l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)2992 ! enddo2993 2994 return2995 end2996 2997 !======================================================================2998 SUBROUTINE interp_astex_time(day,day1,annee_ref &2999 & ,year_ini_astex,day_ini_astex,nt_astex,dt_astex &3000 & ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex &3001 & ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof &3002 & ,ufa_prof,vfa_prof)3003 implicit none3004 3005 !---------------------------------------------------------------------------------------3006 ! Time interpolation of a 2D field to the timestep corresponding to day3007 !3008 ! day: current julian day (e.g. 717538.2)3009 ! day1: first day of the simulation3010 ! nt_astex: total nb of data in the forcing (e.g. 41 for Astex)3011 ! dt_astex: total time interval (in sec) between 2 forcing data (e.g. 1h for Astex)3012 !---------------------------------------------------------------------------------------3013 3014 ! inputs:3015 integer annee_ref3016 integer nt_astex,nlev_astex3017 integer year_ini_astex3018 real day, day1,day_ini_astex,dt_astex3019 real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)3020 real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)3021 ! outputs:3022 real div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof3023 ! local:3024 integer it_astex1, it_astex23025 real timeit,time_astex1,time_astex2,frac3026 3027 ! Check that initial day of the simulation consistent with ASTEX period:3028 if (annee_ref.ne.1992 ) then3029 print*,'Pour Astex, annee_ref doit etre 1992 '3030 print*,'Changer annee_ref dans run.def'3031 stop3032 endif3033 if (annee_ref.eq.1992 .and. day1.lt.day_ini_astex) then3034 print*,'Astex debute le 13 Juin 1992 (jour julien=165)'3035 print*,'Changer dayref dans run.def'3036 stop3037 endif3038 3039 ! Determine timestep relative to the 1st day of TOGA-COARE:3040 ! timeit=(day-day1)*86400.3041 ! if (annee_ref.eq.1992) then3042 ! timeit=(day-day_ini_astex)*86400.3043 ! else3044 ! timeit=(day+2.-1.)*86400. ! 2 days between Jun13 and Jun15 19923045 ! endif3046 timeit=(day-day_ini_astex)*864003047 3048 ! Determine the closest observation times:3049 it_astex1=INT(timeit/dt_astex)+13050 it_astex2=it_astex1 + 13051 time_astex1=(it_astex1-1)*dt_astex3052 time_astex2=(it_astex2-1)*dt_astex3053 print *,'timeit day day_ini_astex',timeit,day,day_ini_astex3054 print *,'it_astex1,it_astex2,time_astex1,time_astex2', &3055 & it_astex1,it_astex2,time_astex1,time_astex23056 3057 if (it_astex1 .ge. nt_astex) then3058 write(*,*) 'PB-stop: day, it_astex1, it_astex2, timeit: ' &3059 & ,day,it_astex1,it_astex2,timeit/86400.3060 stop3061 endif3062 3063 ! time interpolation:3064 frac=(time_astex2-timeit)/(time_astex2-time_astex1)3065 frac=max(frac,0.0)3066 3067 div_prof = div_astex(it_astex2) &3068 & -frac*(div_astex(it_astex2)-div_astex(it_astex1))3069 ts_prof = ts_astex(it_astex2) &3070 & -frac*(ts_astex(it_astex2)-ts_astex(it_astex1))3071 ug_prof = ug_astex(it_astex2) &3072 & -frac*(ug_astex(it_astex2)-ug_astex(it_astex1))3073 vg_prof = vg_astex(it_astex2) &3074 & -frac*(vg_astex(it_astex2)-vg_astex(it_astex1))3075 ufa_prof = ufa_astex(it_astex2) &3076 & -frac*(ufa_astex(it_astex2)-ufa_astex(it_astex1))3077 vfa_prof = vfa_astex(it_astex2) &3078 & -frac*(vfa_astex(it_astex2)-vfa_astex(it_astex1))3079 3080 print*, &3081 &'day,annee_ref,day_ini_astex,timeit,it_astex1,it_astex2,SST:', &3082 &day,annee_ref,day_ini_astex,timeit/86400.,it_astex1, &3083 &it_astex2,div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof3084 3085 return3086 END3087 3088 !======================================================================3089 SUBROUTINE interp_toga_time(day,day1,annee_ref &3090 & ,year_ini_toga,day_ini_toga,nt_toga,dt_toga,nlev_toga &3091 & ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga &3092 & ,ht_toga,vt_toga,hq_toga,vq_toga &3093 & ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof &3094 & ,ht_prof,vt_prof,hq_prof,vq_prof)3095 implicit none3096 3097 !---------------------------------------------------------------------------------------3098 ! Time interpolation of a 2D field to the timestep corresponding to day3099 !3100 ! day: current julian day (e.g. 717538.2)3101 ! day1: first day of the simulation3102 ! nt_toga: total nb of data in the forcing (e.g. 480 for TOGA-COARE)3103 ! dt_toga: total time interval (in sec) between 2 forcing data (e.g. 6h for TOGA-COARE)3104 !---------------------------------------------------------------------------------------3105 3106 #include "compar1d.h"3107 3108 ! inputs:3109 integer annee_ref3110 integer nt_toga,nlev_toga3111 integer year_ini_toga3112 real day, day1,day_ini_toga,dt_toga3113 real ts_toga(nt_toga)3114 real plev_toga(nlev_toga,nt_toga),t_toga(nlev_toga,nt_toga)3115 real q_toga(nlev_toga,nt_toga),u_toga(nlev_toga,nt_toga)3116 real v_toga(nlev_toga,nt_toga),w_toga(nlev_toga,nt_toga)3117 real ht_toga(nlev_toga,nt_toga),vt_toga(nlev_toga,nt_toga)3118 real hq_toga(nlev_toga,nt_toga),vq_toga(nlev_toga,nt_toga)3119 ! outputs:3120 real ts_prof3121 real plev_prof(nlev_toga),t_prof(nlev_toga)3122 real q_prof(nlev_toga),u_prof(nlev_toga)3123 real v_prof(nlev_toga),w_prof(nlev_toga)3124 real ht_prof(nlev_toga),vt_prof(nlev_toga)3125 real hq_prof(nlev_toga),vq_prof(nlev_toga)3126 ! local:3127 integer it_toga1, it_toga2,k3128 real timeit,time_toga1,time_toga2,frac3129 3130 3131 if (forcing_type.eq.2) then3132 ! Check that initial day of the simulation consistent with TOGA-COARE period:3133 if (annee_ref.ne.1992 .and. annee_ref.ne.1993) then3134 print*,'Pour TOGA-COARE, annee_ref doit etre 1992 ou 1993'3135 print*,'Changer annee_ref dans run.def'3136 stop3137 endif3138 if (annee_ref.eq.1992 .and. day1.lt.day_ini_toga) then3139 print*,'TOGA-COARE a debute le 1er Nov 1992 (jour julien=306)'3140 print*,'Changer dayref dans run.def'3141 stop3142 endif3143 if (annee_ref.eq.1993 .and. day1.gt.day_ini_toga+119) then3144 print*,'TOGA-COARE a fini le 28 Feb 1993 (jour julien=59)'3145 print*,'Changer dayref ou nday dans run.def'3146 stop3147 endif3148 3149 else if (forcing_type.eq.4) then3150 3151 ! Check that initial day of the simulation consistent with TWP-ICE period:3152 if (annee_ref.ne.2006) then3153 print*,'Pour TWP-ICE, annee_ref doit etre 2006'3154 print*,'Changer annee_ref dans run.def'3155 stop3156 endif3157 if (annee_ref.eq.2006 .and. day1.lt.day_ini_toga) then3158 print*,'TWP-ICE a debute le 17 Jan 2006 (jour julien=17)'3159 print*,'Changer dayref dans run.def'3160 stop3161 endif3162 if (annee_ref.eq.2006 .and. day1.gt.day_ini_toga+26) then3163 print*,'TWP-ICE a fini le 12 Feb 2006 (jour julien=43)'3164 print*,'Changer dayref ou nday dans run.def'3165 stop3166 endif3167 3168 endif3169 3170 ! Determine timestep relative to the 1st day of TOGA-COARE:3171 ! timeit=(day-day1)*86400.3172 ! if (annee_ref.eq.1992) then3173 ! timeit=(day-day_ini_toga)*86400.3174 ! else3175 ! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 19923176 ! endif3177 timeit=(day-day_ini_toga)*864003178 3179 ! Determine the closest observation times:3180 it_toga1=INT(timeit/dt_toga)+13181 it_toga2=it_toga1 + 13182 time_toga1=(it_toga1-1)*dt_toga3183 time_toga2=(it_toga2-1)*dt_toga3184 3185 if (it_toga1 .ge. nt_toga) then3186 write(*,*) 'PB-stop: day, it_toga1, it_toga2, timeit: ' &3187 & ,day,it_toga1,it_toga2,timeit/86400.3188 stop3189 endif3190 3191 ! time interpolation:3192 frac=(time_toga2-timeit)/(time_toga2-time_toga1)3193 frac=max(frac,0.0)3194 3195 ts_prof = ts_toga(it_toga2) &3196 & -frac*(ts_toga(it_toga2)-ts_toga(it_toga1))3197 3198 ! print*,3199 ! :'day,annee_ref,day_ini_toga,timeit,it_toga1,it_toga2,SST:',3200 ! :day,annee_ref,day_ini_toga,timeit/86400.,it_toga1,it_toga2,ts_prof3201 3202 do k=1,nlev_toga3203 plev_prof(k) = 100.*(plev_toga(k,it_toga2) &3204 & -frac*(plev_toga(k,it_toga2)-plev_toga(k,it_toga1)))3205 t_prof(k) = t_toga(k,it_toga2) &3206 & -frac*(t_toga(k,it_toga2)-t_toga(k,it_toga1))3207 q_prof(k) = q_toga(k,it_toga2) &3208 & -frac*(q_toga(k,it_toga2)-q_toga(k,it_toga1))3209 u_prof(k) = u_toga(k,it_toga2) &3210 & -frac*(u_toga(k,it_toga2)-u_toga(k,it_toga1))3211 v_prof(k) = v_toga(k,it_toga2) &3212 & -frac*(v_toga(k,it_toga2)-v_toga(k,it_toga1))3213 w_prof(k) = w_toga(k,it_toga2) &3214 & -frac*(w_toga(k,it_toga2)-w_toga(k,it_toga1))3215 ht_prof(k) = ht_toga(k,it_toga2) &3216 & -frac*(ht_toga(k,it_toga2)-ht_toga(k,it_toga1))3217 vt_prof(k) = vt_toga(k,it_toga2) &3218 & -frac*(vt_toga(k,it_toga2)-vt_toga(k,it_toga1))3219 hq_prof(k) = hq_toga(k,it_toga2) &3220 & -frac*(hq_toga(k,it_toga2)-hq_toga(k,it_toga1))3221 vq_prof(k) = vq_toga(k,it_toga2) &3222 & -frac*(vq_toga(k,it_toga2)-vq_toga(k,it_toga1))3223 enddo3224 3225 return3226 END3227 3228 !======================================================================3229 SUBROUTINE interp_dice_time(day,day1,annee_ref &3230 & ,year_ini_dice,day_ini_dice,nt_dice,dt_dice &3231 & ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice &3232 & ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice &3233 & ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice &3234 & ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof &3235 & ,ustar_prof,psurf_prof,ug_prof,vg_prof &3236 & ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof)3237 implicit none3238 3239 !---------------------------------------------------------------------------------------3240 ! Time interpolation of a 2D field to the timestep corresponding to day3241 !3242 ! day: current julian day (e.g. 717538.2)3243 ! day1: first day of the simulation3244 ! nt_dice: total nb of data in the forcing (e.g. 145 for Dice)3245 ! dt_dice: total time interval (in sec) between 2 forcing data (e.g. 30min. for Dice)3246 !---------------------------------------------------------------------------------------3247 3248 #include "compar1d.h"3249 3250 ! inputs:3251 integer annee_ref3252 integer nt_dice,nlev_dice3253 integer year_ini_dice3254 real day, day1,day_ini_dice,dt_dice3255 real shf_dice(nt_dice),lhf_dice(nt_dice),lwup_dice(nt_dice)3256 real swup_dice(nt_dice),tg_dice(nt_dice),ustar_dice(nt_dice)3257 real psurf_dice(nt_dice),ug_dice(nt_dice),vg_dice(nt_dice)3258 real ht_dice(nlev_dice,nt_dice),hq_dice(nlev_dice,nt_dice)3259 real hu_dice(nlev_dice,nt_dice),hv_dice(nlev_dice,nt_dice)3260 real w_dice(nlev_dice,nt_dice),omega_dice(nlev_dice,nt_dice)3261 ! outputs:3262 real tg_prof,shf_prof,lhf_prof,lwup_prof,swup_prof3263 real ustar_prof,psurf_prof,ug_prof,vg_prof3264 real ht_prof(nlev_dice),hq_prof(nlev_dice)3265 real hu_prof(nlev_dice),hv_prof(nlev_dice)3266 real w_prof(nlev_dice),omega_prof(nlev_dice)3267 ! local:3268 integer it_dice1, it_dice2,k3269 real timeit,time_dice1,time_dice2,frac3270 3271 3272 if (forcing_type.eq.7) then3273 ! Check that initial day of the simulation consistent with Dice period:3274 print *,'annee_ref=',annee_ref3275 print *,'day1=',day13276 print *,'day_ini_dice=',day_ini_dice3277 if (annee_ref.ne.1999) then3278 print*,'Pour Dice, annee_ref doit etre 1999'3279 print*,'Changer annee_ref dans run.def'3280 stop3281 endif3282 if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice) then3283 print*,'Dice a debute le 23 Oct 1999 (jour julien=296)'3284 print*,'Changer dayref dans run.def',day1,day_ini_dice3285 stop3286 endif3287 if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice+2) then3288 print*,'Dice a fini le 25 Oct 1999 (jour julien=298)'3289 print*,'Changer dayref ou nday dans run.def',day1,day_ini_dice3290 stop3291 endif3292 3293 endif3294 3295 ! Determine timestep relative to the 1st day of TOGA-COARE:3296 ! timeit=(day-day1)*86400.3297 ! if (annee_ref.eq.1992) then3298 ! timeit=(day-day_ini_dice)*86400.3299 ! else3300 ! timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 19923301 ! endif3302 timeit=(day-day_ini_dice)*864003303 3304 ! Determine the closest observation times:3305 it_dice1=INT(timeit/dt_dice)+13306 it_dice2=it_dice1 + 13307 time_dice1=(it_dice1-1)*dt_dice3308 time_dice2=(it_dice2-1)*dt_dice3309 3310 if (it_dice1 .ge. nt_dice) then3311 write(*,*) 'PB-stop: day, it_dice1, it_dice2, timeit: ',day,it_dice1,it_dice2,timeit/86400.3312 stop3313 endif3314 3315 ! time interpolation:3316 frac=(time_dice2-timeit)/(time_dice2-time_dice1)3317 frac=max(frac,0.0)3318 3319 shf_prof = shf_dice(it_dice2)-frac*(shf_dice(it_dice2)-shf_dice(it_dice1))3320 lhf_prof = lhf_dice(it_dice2)-frac*(lhf_dice(it_dice2)-lhf_dice(it_dice1))3321 lwup_prof = lwup_dice(it_dice2)-frac*(lwup_dice(it_dice2)-lwup_dice(it_dice1))3322 swup_prof = swup_dice(it_dice2)-frac*(swup_dice(it_dice2)-swup_dice(it_dice1))3323 tg_prof = tg_dice(it_dice2)-frac*(tg_dice(it_dice2)-tg_dice(it_dice1))3324 ustar_prof = ustar_dice(it_dice2)-frac*(ustar_dice(it_dice2)-ustar_dice(it_dice1))3325 psurf_prof = psurf_dice(it_dice2)-frac*(psurf_dice(it_dice2)-psurf_dice(it_dice1))3326 ug_prof = ug_dice(it_dice2)-frac*(ug_dice(it_dice2)-ug_dice(it_dice1))3327 vg_prof = vg_dice(it_dice2)-frac*(vg_dice(it_dice2)-vg_dice(it_dice1))3328 3329 ! print*,3330 ! :'day,annee_ref,day_ini_dice,timeit,it_dice1,it_dice2,SST:',3331 ! :day,annee_ref,day_ini_dice,timeit/86400.,it_dice1,it_dice2,ts_prof3332 3333 do k=1,nlev_dice3334 ht_prof(k) = ht_dice(k,it_dice2)-frac*(ht_dice(k,it_dice2)-ht_dice(k,it_dice1))3335 hq_prof(k) = hq_dice(k,it_dice2)-frac*(hq_dice(k,it_dice2)-hq_dice(k,it_dice1))3336 hu_prof(k) = hu_dice(k,it_dice2)-frac*(hu_dice(k,it_dice2)-hu_dice(k,it_dice1))3337 hv_prof(k) = hv_dice(k,it_dice2)-frac*(hv_dice(k,it_dice2)-hv_dice(k,it_dice1))3338 w_prof(k) = w_dice(k,it_dice2)-frac*(w_dice(k,it_dice2)-w_dice(k,it_dice1))3339 omega_prof(k) = omega_dice(k,it_dice2)-frac*(omega_dice(k,it_dice2)-omega_dice(k,it_dice1))3340 enddo3341 3342 return3343 END3344 3345 !======================================================================3346 SUBROUTINE interp_gabls4_time(day,day1,annee_ref &3347 & ,year_ini_gabls4,day_ini_gabls4,nt_gabls4,dt_gabls4,nlev_gabls4 &3348 & ,ug_gabls4,vg_gabls4,ht_gabls4,hq_gabls4,tg_gabls4 &3349 & ,ug_prof,vg_prof,ht_prof,hq_prof,tg_prof)3350 implicit none3351 3352 !---------------------------------------------------------------------------------------3353 ! Time interpolation of a 2D field to the timestep corresponding to day3354 !3355 ! day: current julian day3356 ! day1: first day of the simulation3357 ! nt_gabls4: total nb of data in the forcing (e.g. 37 for gabls4)3358 ! dt_gabls4: total time interval (in sec) between 2 forcing data (e.g. 60min. for gabls4)3359 !---------------------------------------------------------------------------------------3360 3361 #include "compar1d.h"3362 3363 ! inputs:3364 integer annee_ref3365 integer nt_gabls4,nlev_gabls43366 integer year_ini_gabls43367 real day, day1,day_ini_gabls4,dt_gabls43368 real ug_gabls4(nlev_gabls4,nt_gabls4),vg_gabls4(nlev_gabls4,nt_gabls4)3369 real ht_gabls4(nlev_gabls4,nt_gabls4),hq_gabls4(nlev_gabls4,nt_gabls4)3370 real tg_gabls4(nt_gabls4), tg_prof3371 ! outputs:3372 real ug_prof(nlev_gabls4),vg_prof(nlev_gabls4)3373 real ht_prof(nlev_gabls4),hq_prof(nlev_gabls4)3374 ! local:3375 integer it_gabls41, it_gabls42,k3376 real timeit,time_gabls41,time_gabls42,frac3377 3378 3379 3380 ! Check that initial day of the simulation consistent with gabls4 period:3381 if (forcing_type.eq.8 ) then3382 print *,'annee_ref=',annee_ref3383 print *,'day1=',day13384 print *,'day_ini_gabls4=',day_ini_gabls43385 if (annee_ref.ne.2009) then3386 print*,'Pour gabls4, annee_ref doit etre 2009'3387 print*,'Changer annee_ref dans run.def'3388 stop3389 endif3390 if (annee_ref.eq.2009 .and. day1.gt.day_ini_gabls4) then3391 print*,'gabls4 a debute le 11 dec 2009 (jour julien=345)'3392 print*,'Changer dayref dans run.def',day1,day_ini_gabls43393 stop3394 endif3395 if (annee_ref.eq.2009 .and. day1.gt.day_ini_gabls4+2) then3396 print*,'gabls4 a fini le 12 dec 2009 (jour julien=346)'3397 print*,'Changer dayref ou nday dans run.def',day1,day_ini_gabls43398 stop3399 endif3400 endif3401 3402 timeit=(day-day_ini_gabls4)*864003403 print *,'day,day_ini_gabls4=',day,day_ini_gabls43404 print *,'nt_gabls4,dt,timeit=',nt_gabls4,dt_gabls4,timeit3405 3406 ! Determine the closest observation times:3407 it_gabls41=INT(timeit/dt_gabls4)+13408 it_gabls42=it_gabls41 + 13409 time_gabls41=(it_gabls41-1)*dt_gabls43410 time_gabls42=(it_gabls42-1)*dt_gabls43411 3412 if (it_gabls41 .ge. nt_gabls4) then3413 write(*,*) 'PB-stop: day, it_gabls41, it_gabls42, timeit: ',day,it_gabls41,it_gabls42,timeit/86400.3414 stop3415 endif3416 3417 ! time interpolation:3418 frac=(time_gabls42-timeit)/(time_gabls42-time_gabls41)3419 frac=max(frac,0.0)3420 3421 3422 do k=1,nlev_gabls43423 ug_prof(k) = ug_gabls4(k,it_gabls42)-frac*(ug_gabls4(k,it_gabls42)-ug_gabls4(k,it_gabls41))3424 vg_prof(k) = vg_gabls4(k,it_gabls42)-frac*(vg_gabls4(k,it_gabls42)-vg_gabls4(k,it_gabls41))3425 ht_prof(k) = ht_gabls4(k,it_gabls42)-frac*(ht_gabls4(k,it_gabls42)-ht_gabls4(k,it_gabls41))3426 hq_prof(k) = hq_gabls4(k,it_gabls42)-frac*(hq_gabls4(k,it_gabls42)-hq_gabls4(k,it_gabls41))3427 enddo3428 tg_prof=tg_gabls4(it_gabls42)-frac*(tg_gabls4(it_gabls42)-tg_gabls4(it_gabls41))3429 return3430 END3431 3432 !======================================================================3433 SUBROUTINE interp_armcu_time(day,day1,annee_ref &3434 & ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu &3435 & ,nlev_armcu,fs_armcu,fl_armcu,at_armcu,rt_armcu &3436 & ,aqt_armcu,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof)3437 implicit none3438 3439 !---------------------------------------------------------------------------------------3440 ! Time interpolation of a 2D field to the timestep corresponding to day3441 !3442 ! day: current julian day (e.g. 717538.2)3443 ! day1: first day of the simulation3444 ! nt_armcu: total nb of data in the forcing (e.g. 31 for armcu)3445 ! dt_armcu: total time interval (in sec) between 2 forcing data (e.g. 1/2h for armcu)3446 ! fs= sensible flux3447 ! fl= latent flux3448 ! at,rt,aqt= advective and radiative tendencies3449 !---------------------------------------------------------------------------------------3450 3451 ! inputs:3452 integer annee_ref3453 integer nt_armcu,nlev_armcu3454 integer year_ini_armcu3455 real day, day1,day_ini_armcu,dt_armcu3456 real fs_armcu(nt_armcu),fl_armcu(nt_armcu),at_armcu(nt_armcu)3457 real rt_armcu(nt_armcu),aqt_armcu(nt_armcu)3458 ! outputs:3459 real fs_prof,fl_prof,at_prof,rt_prof,aqt_prof3460 ! local:3461 integer it_armcu1, it_armcu2,k3462 real timeit,time_armcu1,time_armcu2,frac3463 3464 ! Check that initial day of the simulation consistent with ARMCU period:3465 if (annee_ref.ne.1997 ) then3466 print*,'Pour ARMCU, annee_ref doit etre 1997 '3467 print*,'Changer annee_ref dans run.def'3468 stop3469 endif3470 3471 timeit=(day-day_ini_armcu)*864003472 3473 ! Determine the closest observation times:3474 it_armcu1=INT(timeit/dt_armcu)+13475 it_armcu2=it_armcu1 + 13476 time_armcu1=(it_armcu1-1)*dt_armcu3477 time_armcu2=(it_armcu2-1)*dt_armcu3478 print *,'timeit day day_ini_armcu',timeit,day,day_ini_armcu3479 print *,'it_armcu1,it_armcu2,time_armcu1,time_armcu2', &3480 & it_armcu1,it_armcu2,time_armcu1,time_armcu23481 3482 if (it_armcu1 .ge. nt_armcu) then3483 write(*,*) 'PB-stop: day, it_armcu1, it_armcu2, timeit: ' &3484 & ,day,it_armcu1,it_armcu2,timeit/86400.3485 stop3486 endif3487 3488 ! time interpolation:3489 frac=(time_armcu2-timeit)/(time_armcu2-time_armcu1)3490 frac=max(frac,0.0)3491 3492 fs_prof = fs_armcu(it_armcu2) &3493 & -frac*(fs_armcu(it_armcu2)-fs_armcu(it_armcu1))3494 fl_prof = fl_armcu(it_armcu2) &3495 & -frac*(fl_armcu(it_armcu2)-fl_armcu(it_armcu1))3496 at_prof = at_armcu(it_armcu2) &3497 & -frac*(at_armcu(it_armcu2)-at_armcu(it_armcu1))3498 rt_prof = rt_armcu(it_armcu2) &3499 & -frac*(rt_armcu(it_armcu2)-rt_armcu(it_armcu1))3500 aqt_prof = aqt_armcu(it_armcu2) &3501 & -frac*(aqt_armcu(it_armcu2)-aqt_armcu(it_armcu1))3502 3503 print*, &3504 &'day,annee_ref,day_ini_armcu,timeit,it_armcu1,it_armcu2,SST:', &3505 &day,annee_ref,day_ini_armcu,timeit/86400.,it_armcu1, &3506 &it_armcu2,fs_prof,fl_prof,at_prof,rt_prof,aqt_prof3507 3508 return3509 END3510 3511 !=====================================================================3512 subroutine readprofiles(nlev_max,kmax,ntrac,height, &3513 & thlprof,qtprof,uprof, &3514 & vprof,e12prof,ugprof,vgprof, &3515 & wfls,dqtdxls,dqtdyls,dqtdtls, &3516 & thlpcar,tracer,nt1,nt2)3517 implicit none3518 3519 integer nlev_max,kmax,kmax2,ntrac3520 logical :: llesread = .true.3521 3522 real height(nlev_max),thlprof(nlev_max),qtprof(nlev_max), &3523 & uprof(nlev_max),vprof(nlev_max),e12prof(nlev_max), &3524 & ugprof(nlev_max),vgprof(nlev_max),wfls(nlev_max), &3525 & dqtdxls(nlev_max),dqtdyls(nlev_max),dqtdtls(nlev_max), &3526 & thlpcar(nlev_max),tracer(nlev_max,ntrac)3527 3528 real height1(nlev_max)3529 3530 integer, parameter :: ilesfile=13531 integer :: ierr,k,itrac,nt1,nt23532 3533 if(.not.(llesread)) return3534 3535 open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)3536 if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'3537 read (ilesfile,*) kmax3538 do k=1,kmax3539 read (ilesfile,*) height1(k),thlprof(k),qtprof (k), &3540 & uprof (k),vprof (k),e12prof(k)3541 enddo3542 close(ilesfile)3543 3544 open(ilesfile,file='lscale.inp.001',status='old',iostat=ierr)3545 if (ierr /= 0) stop 'ERROR:Lscale.inp does not exist'3546 read (ilesfile,*) kmax23547 if (kmax .ne. kmax2) then3548 print *, 'fichiers prof.inp et lscale.inp incompatibles :'3549 print *, 'nbre de niveaux : ',kmax,' et ',kmax23550 stop 'lecture profiles'3551 endif3552 do k=1,kmax3553 read (ilesfile,*) height(k),ugprof(k),vgprof(k),wfls(k), &3554 & dqtdxls(k),dqtdyls(k),dqtdtls(k),thlpcar(k)3555 end do3556 do k=1,kmax3557 if (height(k) .ne. height1(k)) then3558 print *, 'fichiers prof.inp et lscale.inp incompatibles :'3559 print *, 'les niveaux different : ',k,height1(k), height(k)3560 stop3561 endif3562 end do3563 close(ilesfile)3564 3565 open(ilesfile,file='trac.inp.001',status='old',iostat=ierr)3566 if (ierr /= 0) then3567 print*,'WARNING : trac.inp does not exist'3568 else3569 read (ilesfile,*) kmax2,nt1,nt23570 if (nt2>ntrac) then3571 stop 'Augmenter le nombre de traceurs dans traceur.def'3572 endif3573 if (kmax .ne. kmax2) then3574 print *, 'fichiers prof.inp et lscale.inp incompatibles :'3575 print *, 'nbre de niveaux : ',kmax,' et ',kmax23576 stop 'lecture profiles'3577 endif3578 do k=1,kmax3579 read (ilesfile,*) height(k),(tracer(k,itrac),itrac=nt1,nt2)3580 end do3581 close(ilesfile)3582 endif3583 3584 return3585 end3586 !======================================================================3587 subroutine readprofile_sandu(nlev_max,kmax,height,pprof,tprof, &3588 & thlprof,qprof,uprof,vprof,wprof,omega,o3mmr)3589 !======================================================================3590 implicit none3591 3592 integer nlev_max,kmax3593 logical :: llesread = .true.3594 3595 real height(nlev_max),pprof(nlev_max),tprof(nlev_max)3596 real thlprof(nlev_max)3597 real qprof(nlev_max),uprof(nlev_max),vprof(nlev_max)3598 real wprof(nlev_max),omega(nlev_max),o3mmr(nlev_max)3599 3600 integer, parameter :: ilesfile=13601 integer :: k,ierr3602 3603 if(.not.(llesread)) return3604 3605 open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)3606 if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'3607 read (ilesfile,*) kmax3608 do k=1,kmax3609 read (ilesfile,*) height(k),pprof(k), tprof(k),thlprof(k), &3610 & qprof (k),uprof(k), vprof(k), wprof(k), &3611 & omega (k),o3mmr(k)3612 enddo3613 close(ilesfile)3614 3615 return3616 end3617 3618 !======================================================================3619 subroutine readprofile_astex(nlev_max,kmax,height,pprof,tprof, &3620 & thlprof,qvprof,qlprof,qtprof,uprof,vprof,wprof,tkeprof,o3mmr)3621 !======================================================================3622 implicit none3623 3624 integer nlev_max,kmax3625 logical :: llesread = .true.3626 3627 real height(nlev_max),pprof(nlev_max),tprof(nlev_max), &3628 & thlprof(nlev_max),qlprof(nlev_max),qtprof(nlev_max), &3629 & qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max), &3630 & wprof(nlev_max),tkeprof(nlev_max),o3mmr(nlev_max)3631 3632 integer, parameter :: ilesfile=13633 integer :: ierr,k3634 3635 if(.not.(llesread)) return3636 3637 open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)3638 if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'3639 read (ilesfile,*) kmax3640 do k=1,kmax3641 read (ilesfile,*) height(k),pprof(k), tprof(k),thlprof(k), &3642 & qvprof (k),qlprof (k),qtprof (k), &3643 & uprof(k), vprof(k), wprof(k),tkeprof(k),o3mmr(k)3644 enddo3645 close(ilesfile)3646 3647 return3648 end3649 3650 3651 3652 !======================================================================3653 subroutine readprofile_armcu(nlev_max,kmax,height,pprof,uprof, &3654 & vprof,thetaprof,tprof,qvprof,rvprof,aprof,bprof)3655 !======================================================================3656 implicit none3657 3658 integer nlev_max,kmax3659 logical :: llesread = .true.3660 3661 real height(nlev_max),pprof(nlev_max),tprof(nlev_max)3662 real thetaprof(nlev_max),rvprof(nlev_max)3663 real qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max)3664 real aprof(nlev_max+1),bprof(nlev_max+1)3665 3666 integer, parameter :: ilesfile=13667 integer, parameter :: ifile=23668 integer :: ierr,jtot,k3669 3670 if(.not.(llesread)) return3671 3672 ! Read profiles at full levels3673 IF(nlev_max.EQ.19) THEN3674 open (ilesfile,file='prof.inp.19',status='old',iostat=ierr)3675 print *,'On ouvre prof.inp.19'3676 ELSE3677 open (ilesfile,file='prof.inp.40',status='old',iostat=ierr)3678 print *,'On ouvre prof.inp.40'3679 ENDIF3680 if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'3681 read (ilesfile,*) kmax3682 do k=1,kmax3683 read (ilesfile,*) height(k) ,pprof(k), uprof(k), vprof(k), &3684 & thetaprof(k) ,tprof(k), qvprof(k),rvprof(k)3685 enddo3686 close(ilesfile)3687 3688 ! Vertical coordinates half levels for eta-coordinates (plev = alpha + beta * psurf)3689 IF(nlev_max.EQ.19) THEN3690 open (ifile,file='proh.inp.19',status='old',iostat=ierr)3691 print *,'On ouvre proh.inp.19'3692 if (ierr /= 0) stop 'ERROR:Proh.inp.19 does not exist'3693 ELSE3694 open (ifile,file='proh.inp.40',status='old',iostat=ierr)3695 print *,'On ouvre proh.inp.40'3696 if (ierr /= 0) stop 'ERROR:Proh.inp.40 does not exist'3697 ENDIF3698 read (ifile,*) kmax3699 do k=1,kmax3700 read (ifile,*) jtot,aprof(k),bprof(k)3701 enddo3702 close(ifile)3703 3704 return3705 end3706 3707 !=====================================================================3708 subroutine read_fire(fich_fire,nlevel,ntime &3709 & ,zz,thl,qt,u,v,tke &3710 & ,ug,vg,wls,dqtdx,dqtdy,dqtdt,thl_rad)3711 3712 !program reading forcings of the FIRE case study3713 3714 3715 implicit none3716 3717 #include "netcdf.inc"3718 3719 integer ntime,nlevel3720 character*80 :: fich_fire3721 real*8 zz(nlevel)3722 3723 real*8 thl(nlevel)3724 real*8 qt(nlevel),u(nlevel)3725 real*8 v(nlevel),tke(nlevel)3726 real*8 ug(nlevel,ntime),vg(nlevel,ntime),wls(nlevel,ntime)3727 real*8 dqtdx(nlevel,ntime),dqtdy(nlevel,ntime)3728 real*8 dqtdt(nlevel,ntime),thl_rad(nlevel,ntime)3729 3730 integer nid, ierr3731 integer nbvar3d3732 parameter(nbvar3d=30)3733 integer var3didin(nbvar3d)3734 3735 ierr = NF_OPEN(fich_fire,NF_NOWRITE,nid)3736 if (ierr.NE.NF_NOERR) then3737 write(*,*) 'ERROR: Pb opening forcings nc file '3738 write(*,*) NF_STRERROR(ierr)3739 stop ""3740 endif3741 3742 3743 ierr=NF_INQ_VARID(nid,"zz",var3didin(1))3744 if(ierr/=NF_NOERR) then3745 write(*,*) NF_STRERROR(ierr)3746 stop 'lev'3747 endif3748 3749 3750 ierr=NF_INQ_VARID(nid,"thetal",var3didin(2))3751 if(ierr/=NF_NOERR) then3752 write(*,*) NF_STRERROR(ierr)3753 stop 'temp'3754 endif3755 3756 ierr=NF_INQ_VARID(nid,"qt",var3didin(3))3757 if(ierr/=NF_NOERR) then3758 write(*,*) NF_STRERROR(ierr)3759 stop 'qv'3760 endif3761 3762 ierr=NF_INQ_VARID(nid,"u",var3didin(4))3763 if(ierr/=NF_NOERR) then3764 write(*,*) NF_STRERROR(ierr)3765 stop 'u'3766 endif3767 3768 ierr=NF_INQ_VARID(nid,"v",var3didin(5))3769 if(ierr/=NF_NOERR) then3770 write(*,*) NF_STRERROR(ierr)3771 stop 'v'3772 endif3773 3774 ierr=NF_INQ_VARID(nid,"tke",var3didin(6))3775 if(ierr/=NF_NOERR) then3776 write(*,*) NF_STRERROR(ierr)3777 stop 'tke'3778 endif3779 3780 ierr=NF_INQ_VARID(nid,"ugeo",var3didin(7))3781 if(ierr/=NF_NOERR) then3782 write(*,*) NF_STRERROR(ierr)3783 stop 'ug'3784 endif3785 3786 ierr=NF_INQ_VARID(nid,"vgeo",var3didin(8))3787 if(ierr/=NF_NOERR) then3788 write(*,*) NF_STRERROR(ierr)3789 stop 'vg'3790 endif3791 3792 ierr=NF_INQ_VARID(nid,"wls",var3didin(9))3793 if(ierr/=NF_NOERR) then3794 write(*,*) NF_STRERROR(ierr)3795 stop 'wls'3796 endif3797 3798 ierr=NF_INQ_VARID(nid,"dqtdx",var3didin(10))3799 if(ierr/=NF_NOERR) then3800 write(*,*) NF_STRERROR(ierr)3801 stop 'dqtdx'3802 endif3803 3804 ierr=NF_INQ_VARID(nid,"dqtdy",var3didin(11))3805 if(ierr/=NF_NOERR) then3806 write(*,*) NF_STRERROR(ierr)3807 stop 'dqtdy'3808 endif3809 3810 ierr=NF_INQ_VARID(nid,"dqtdt",var3didin(12))3811 if(ierr/=NF_NOERR) then3812 write(*,*) NF_STRERROR(ierr)3813 stop 'dqtdt'3814 endif3815 3816 ierr=NF_INQ_VARID(nid,"thl_rad",var3didin(13))3817 if(ierr/=NF_NOERR) then3818 write(*,*) NF_STRERROR(ierr)3819 stop 'thl_rad'3820 endif3821 !dimensions lecture3822 ! call catchaxis(nid,ntime,nlevel,time,z,ierr)3823 3824 #ifdef NC_DOUBLE3825 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)3826 #else3827 ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)3828 #endif3829 if(ierr/=NF_NOERR) then3830 write(*,*) NF_STRERROR(ierr)3831 stop "getvarup"3832 endif3833 ! write(*,*)'lecture z ok',zz3834 3835 #ifdef NC_DOUBLE3836 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),thl)3837 #else3838 ierr = NF_GET_VAR_REAL(nid,var3didin(2),thl)3839 #endif3840 if(ierr/=NF_NOERR) then3841 write(*,*) NF_STRERROR(ierr)3842 stop "getvarup"3843 endif3844 ! write(*,*)'lecture thl ok',thl3845 3846 #ifdef NC_DOUBLE3847 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qt)3848 #else3849 ierr = NF_GET_VAR_REAL(nid,var3didin(3),qt)3850 #endif3851 if(ierr/=NF_NOERR) then3852 write(*,*) NF_STRERROR(ierr)3853 stop "getvarup"3854 endif3855 ! write(*,*)'lecture qt ok',qt3856 3857 #ifdef NC_DOUBLE3858 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)3859 #else3860 ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)3861 #endif3862 if(ierr/=NF_NOERR) then3863 write(*,*) NF_STRERROR(ierr)3864 stop "getvarup"3865 endif3866 ! write(*,*)'lecture u ok',u3867 3868 #ifdef NC_DOUBLE3869 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)3870 #else3871 ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)3872 #endif3873 if(ierr/=NF_NOERR) then3874 write(*,*) NF_STRERROR(ierr)3875 stop "getvarup"3876 endif3877 ! write(*,*)'lecture v ok',v3878 3879 #ifdef NC_DOUBLE3880 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tke)3881 #else3882 ierr = NF_GET_VAR_REAL(nid,var3didin(6),tke)3883 #endif3884 if(ierr/=NF_NOERR) then3885 write(*,*) NF_STRERROR(ierr)3886 stop "getvarup"3887 endif3888 ! write(*,*)'lecture tke ok',tke3889 3890 #ifdef NC_DOUBLE3891 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ug)3892 #else3893 ierr = NF_GET_VAR_REAL(nid,var3didin(7),ug)3894 #endif3895 if(ierr/=NF_NOERR) then3896 write(*,*) NF_STRERROR(ierr)3897 stop "getvarup"3898 endif3899 ! write(*,*)'lecture ug ok',ug3900 3901 #ifdef NC_DOUBLE3902 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),vg)3903 #else3904 ierr = NF_GET_VAR_REAL(nid,var3didin(8),vg)3905 #endif3906 if(ierr/=NF_NOERR) then3907 write(*,*) NF_STRERROR(ierr)3908 stop "getvarup"3909 endif3910 ! write(*,*)'lecture vg ok',vg3911 3912 #ifdef NC_DOUBLE3913 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),wls)3914 #else3915 ierr = NF_GET_VAR_REAL(nid,var3didin(9),wls)3916 #endif3917 if(ierr/=NF_NOERR) then3918 write(*,*) NF_STRERROR(ierr)3919 stop "getvarup"3920 endif3921 ! write(*,*)'lecture wls ok',wls3922 3923 #ifdef NC_DOUBLE3924 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),dqtdx)3925 #else3926 ierr = NF_GET_VAR_REAL(nid,var3didin(10),dqtdx)3927 #endif3928 if(ierr/=NF_NOERR) then3929 write(*,*) NF_STRERROR(ierr)3930 stop "getvarup"3931 endif3932 ! write(*,*)'lecture dqtdx ok',dqtdx3933 3934 #ifdef NC_DOUBLE3935 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),dqtdy)3936 #else3937 ierr = NF_GET_VAR_REAL(nid,var3didin(11),dqtdy)3938 #endif3939 if(ierr/=NF_NOERR) then3940 write(*,*) NF_STRERROR(ierr)3941 stop "getvarup"3942 endif3943 ! write(*,*)'lecture dqtdy ok',dqtdy3944 3945 #ifdef NC_DOUBLE3946 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),dqtdt)3947 #else3948 ierr = NF_GET_VAR_REAL(nid,var3didin(12),dqtdt)3949 #endif3950 if(ierr/=NF_NOERR) then3951 write(*,*) NF_STRERROR(ierr)3952 stop "getvarup"3953 endif3954 ! write(*,*)'lecture dqtdt ok',dqtdt3955 3956 #ifdef NC_DOUBLE3957 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),thl_rad)3958 #else3959 ierr = NF_GET_VAR_REAL(nid,var3didin(13),thl_rad)3960 #endif3961 if(ierr/=NF_NOERR) then3962 write(*,*) NF_STRERROR(ierr)3963 stop "getvarup"3964 endif3965 ! write(*,*)'lecture thl_rad ok',thl_rad3966 3967 return3968 end subroutine read_fire3969 !=====================================================================3970 subroutine read_dice(fich_dice,nlevel,ntime &3971 & ,zz,pres,t,qv,u,v,o3 &3972 & ,shf,lhf,lwup,swup,tg,ustar,psurf,ug,vg &3973 & ,hadvt,hadvq,hadvu,hadvv,w,omega)3974 3975 !program reading initial profils and forcings of the Dice case study3976 3977 3978 implicit none3979 3980 #include "netcdf.inc"3981 #include "YOMCST.h"3982 3983 integer ntime,nlevel3984 integer l,k3985 character*80 :: fich_dice3986 real*8 time(ntime)3987 real*8 zz(nlevel)3988 3989 real*8 th(nlevel),pres(nlevel),t(nlevel)3990 real*8 qv(nlevel),u(nlevel),v(nlevel),o3(nlevel)3991 real*8 shf(ntime),lhf(ntime),lwup(ntime),swup(ntime),tg(ntime)3992 real*8 ustar(ntime),psurf(ntime),ug(ntime),vg(ntime)3993 real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime),hadvu(nlevel,ntime)3994 real*8 hadvv(nlevel,ntime),w(nlevel,ntime),omega(nlevel,ntime)3995 real*8 pzero3996 3997 integer nid, ierr3998 integer nbvar3d3999 parameter(nbvar3d=30)4000 integer var3didin(nbvar3d)4001 4002 pzero=100000.4003 ierr = NF_OPEN(fich_dice,NF_NOWRITE,nid)4004 if (ierr.NE.NF_NOERR) then4005 write(*,*) 'ERROR: Pb opening forcings nc file '4006 write(*,*) NF_STRERROR(ierr)4007 stop ""4008 endif4009 4010 4011 ierr=NF_INQ_VARID(nid,"height",var3didin(1))4012 if(ierr/=NF_NOERR) then4013 write(*,*) NF_STRERROR(ierr)4014 stop 'height'4015 endif4016 4017 ierr=NF_INQ_VARID(nid,"pf",var3didin(11))4018 if(ierr/=NF_NOERR) then4019 write(*,*) NF_STRERROR(ierr)4020 stop 'pf'4021 endif4022 4023 ierr=NF_INQ_VARID(nid,"theta",var3didin(12))4024 if(ierr/=NF_NOERR) then4025 write(*,*) NF_STRERROR(ierr)4026 stop 'theta'4027 endif4028 4029 ierr=NF_INQ_VARID(nid,"qv",var3didin(13))4030 if(ierr/=NF_NOERR) then4031 write(*,*) NF_STRERROR(ierr)4032 stop 'qv'4033 endif4034 4035 ierr=NF_INQ_VARID(nid,"u",var3didin(14))4036 if(ierr/=NF_NOERR) then4037 write(*,*) NF_STRERROR(ierr)4038 stop 'u'4039 endif4040 4041 ierr=NF_INQ_VARID(nid,"v",var3didin(15))4042 if(ierr/=NF_NOERR) then4043 write(*,*) NF_STRERROR(ierr)4044 stop 'v'4045 endif4046 4047 ierr=NF_INQ_VARID(nid,"o3mmr",var3didin(16))4048 if(ierr/=NF_NOERR) then4049 write(*,*) NF_STRERROR(ierr)4050 stop 'o3'4051 endif4052 4053 ierr=NF_INQ_VARID(nid,"shf",var3didin(2))4054 if(ierr/=NF_NOERR) then4055 write(*,*) NF_STRERROR(ierr)4056 stop 'shf'4057 endif4058 4059 ierr=NF_INQ_VARID(nid,"lhf",var3didin(3))4060 if(ierr/=NF_NOERR) then4061 write(*,*) NF_STRERROR(ierr)4062 stop 'lhf'4063 endif4064 4065 ierr=NF_INQ_VARID(nid,"lwup",var3didin(4))4066 if(ierr/=NF_NOERR) then4067 write(*,*) NF_STRERROR(ierr)4068 stop 'lwup'4069 endif4070 4071 ierr=NF_INQ_VARID(nid,"swup",var3didin(5))4072 if(ierr/=NF_NOERR) then4073 write(*,*) NF_STRERROR(ierr)4074 stop 'dqtdx'4075 endif4076 4077 ierr=NF_INQ_VARID(nid,"Tg",var3didin(6))4078 if(ierr/=NF_NOERR) then4079 write(*,*) NF_STRERROR(ierr)4080 stop 'Tg'4081 endif4082 4083 ierr=NF_INQ_VARID(nid,"ustar",var3didin(7))4084 if(ierr/=NF_NOERR) then4085 write(*,*) NF_STRERROR(ierr)4086 stop 'ustar'4087 endif4088 4089 ierr=NF_INQ_VARID(nid,"psurf",var3didin(8))4090 if(ierr/=NF_NOERR) then4091 write(*,*) NF_STRERROR(ierr)4092 stop 'psurf'4093 endif4094 4095 ierr=NF_INQ_VARID(nid,"Ug",var3didin(9))4096 if(ierr/=NF_NOERR) then4097 write(*,*) NF_STRERROR(ierr)4098 stop 'Ug'4099 endif4100 4101 ierr=NF_INQ_VARID(nid,"Vg",var3didin(10))4102 if(ierr/=NF_NOERR) then4103 write(*,*) NF_STRERROR(ierr)4104 stop 'Vg'4105 endif4106 4107 ierr=NF_INQ_VARID(nid,"hadvT",var3didin(17))4108 if(ierr/=NF_NOERR) then4109 write(*,*) NF_STRERROR(ierr)4110 stop 'hadvT'4111 endif4112 4113 ierr=NF_INQ_VARID(nid,"hadvq",var3didin(18))4114 if(ierr/=NF_NOERR) then4115 write(*,*) NF_STRERROR(ierr)4116 stop 'hadvq'4117 endif4118 4119 ierr=NF_INQ_VARID(nid,"hadvu",var3didin(19))4120 if(ierr/=NF_NOERR) then4121 write(*,*) NF_STRERROR(ierr)4122 stop 'hadvu'4123 endif4124 4125 ierr=NF_INQ_VARID(nid,"hadvv",var3didin(20))4126 if(ierr/=NF_NOERR) then4127 write(*,*) NF_STRERROR(ierr)4128 stop 'hadvv'4129 endif4130 4131 ierr=NF_INQ_VARID(nid,"w",var3didin(21))4132 if(ierr/=NF_NOERR) then4133 write(*,*) NF_STRERROR(ierr)4134 stop 'w'4135 endif4136 4137 ierr=NF_INQ_VARID(nid,"omega",var3didin(22))4138 if(ierr/=NF_NOERR) then4139 write(*,*) NF_STRERROR(ierr)4140 stop 'omega'4141 endif4142 !dimensions lecture4143 ! call catchaxis(nid,ntime,nlevel,time,z,ierr)4144 4145 #ifdef NC_DOUBLE4146 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)4147 #else4148 ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)4149 #endif4150 if(ierr/=NF_NOERR) then4151 write(*,*) NF_STRERROR(ierr)4152 stop "getvarup"4153 endif4154 ! write(*,*)'lecture zz ok',zz4155 4156 #ifdef NC_DOUBLE4157 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pres)4158 #else4159 ierr = NF_GET_VAR_REAL(nid,var3didin(11),pres)4160 #endif4161 if(ierr/=NF_NOERR) then4162 write(*,*) NF_STRERROR(ierr)4163 stop "getvarup"4164 endif4165 ! write(*,*)'lecture pres ok',pres4166 4167 #ifdef NC_DOUBLE4168 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),th)4169 #else4170 ierr = NF_GET_VAR_REAL(nid,var3didin(12),th)4171 #endif4172 if(ierr/=NF_NOERR) then4173 write(*,*) NF_STRERROR(ierr)4174 stop "getvarup"4175 endif4176 ! write(*,*)'lecture th ok',th4177 do k=1,nlevel4178 t(k)=th(k)*(pres(k)/pzero)**rkappa4179 enddo4180 4181 #ifdef NC_DOUBLE4182 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),qv)4183 #else4184 ierr = NF_GET_VAR_REAL(nid,var3didin(13),qv)4185 #endif4186 if(ierr/=NF_NOERR) then4187 write(*,*) NF_STRERROR(ierr)4188 stop "getvarup"4189 endif4190 ! write(*,*)'lecture qv ok',qv4191 4192 #ifdef NC_DOUBLE4193 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),u)4194 #else4195 ierr = NF_GET_VAR_REAL(nid,var3didin(14),u)4196 #endif4197 if(ierr/=NF_NOERR) then4198 write(*,*) NF_STRERROR(ierr)4199 stop "getvarup"4200 endif4201 ! write(*,*)'lecture u ok',u4202 4203 #ifdef NC_DOUBLE4204 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),v)4205 #else4206 ierr = NF_GET_VAR_REAL(nid,var3didin(15),v)4207 #endif4208 if(ierr/=NF_NOERR) then4209 write(*,*) NF_STRERROR(ierr)4210 stop "getvarup"4211 endif4212 ! write(*,*)'lecture v ok',v4213 4214 #ifdef NC_DOUBLE4215 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),o3)4216 #else4217 ierr = NF_GET_VAR_REAL(nid,var3didin(16),o3)4218 #endif4219 if(ierr/=NF_NOERR) then4220 write(*,*) NF_STRERROR(ierr)4221 stop "getvarup"4222 endif4223 ! write(*,*)'lecture o3 ok',o34224 4225 #ifdef NC_DOUBLE4226 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),shf)4227 #else4228 ierr = NF_GET_VAR_REAL(nid,var3didin(2),shf)4229 #endif4230 if(ierr/=NF_NOERR) then4231 write(*,*) NF_STRERROR(ierr)4232 stop "getvarup"4233 endif4234 ! write(*,*)'lecture shf ok',shf4235 4236 #ifdef NC_DOUBLE4237 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),lhf)4238 #else4239 ierr = NF_GET_VAR_REAL(nid,var3didin(3),lhf)4240 #endif4241 if(ierr/=NF_NOERR) then4242 write(*,*) NF_STRERROR(ierr)4243 stop "getvarup"4244 endif4245 ! write(*,*)'lecture lhf ok',lhf4246 4247 #ifdef NC_DOUBLE4248 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),lwup)4249 #else4250 ierr = NF_GET_VAR_REAL(nid,var3didin(4),lwup)4251 #endif4252 if(ierr/=NF_NOERR) then4253 write(*,*) NF_STRERROR(ierr)4254 stop "getvarup"4255 endif4256 ! write(*,*)'lecture lwup ok',lwup4257 4258 #ifdef NC_DOUBLE4259 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),swup)4260 #else4261 ierr = NF_GET_VAR_REAL(nid,var3didin(5),swup)4262 #endif4263 if(ierr/=NF_NOERR) then4264 write(*,*) NF_STRERROR(ierr)4265 stop "getvarup"4266 endif4267 ! write(*,*)'lecture swup ok',swup4268 4269 #ifdef NC_DOUBLE4270 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tg)4271 #else4272 ierr = NF_GET_VAR_REAL(nid,var3didin(6),tg)4273 #endif4274 if(ierr/=NF_NOERR) then4275 write(*,*) NF_STRERROR(ierr)4276 stop "getvarup"4277 endif4278 ! write(*,*)'lecture tg ok',tg4279 4280 #ifdef NC_DOUBLE4281 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ustar)4282 #else4283 ierr = NF_GET_VAR_REAL(nid,var3didin(7),ustar)4284 #endif4285 if(ierr/=NF_NOERR) then4286 write(*,*) NF_STRERROR(ierr)4287 stop "getvarup"4288 endif4289 ! write(*,*)'lecture ustar ok',ustar4290 4291 #ifdef NC_DOUBLE4292 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),psurf)4293 #else4294 ierr = NF_GET_VAR_REAL(nid,var3didin(8),psurf)4295 #endif4296 if(ierr/=NF_NOERR) then4297 write(*,*) NF_STRERROR(ierr)4298 stop "getvarup"4299 endif4300 ! write(*,*)'lecture psurf ok',psurf4301 4302 #ifdef NC_DOUBLE4303 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),ug)4304 #else4305 ierr = NF_GET_VAR_REAL(nid,var3didin(9),ug)4306 #endif4307 if(ierr/=NF_NOERR) then4308 write(*,*) NF_STRERROR(ierr)4309 stop "getvarup"4310 endif4311 ! write(*,*)'lecture ug ok',ug4312 4313 #ifdef NC_DOUBLE4314 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),vg)4315 #else4316 ierr = NF_GET_VAR_REAL(nid,var3didin(10),vg)4317 #endif4318 if(ierr/=NF_NOERR) then4319 write(*,*) NF_STRERROR(ierr)4320 stop "getvarup"4321 endif4322 ! write(*,*)'lecture vg ok',vg4323 4324 #ifdef NC_DOUBLE4325 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),hadvt)4326 #else4327 ierr = NF_GET_VAR_REAL(nid,var3didin(17),hadvt)4328 #endif4329 if(ierr/=NF_NOERR) then4330 write(*,*) NF_STRERROR(ierr)4331 stop "getvarup"4332 endif4333 ! write(*,*)'lecture hadvt ok',hadvt4334 4335 #ifdef NC_DOUBLE4336 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),hadvq)4337 #else4338 ierr = NF_GET_VAR_REAL(nid,var3didin(18),hadvq)4339 #endif4340 if(ierr/=NF_NOERR) then4341 write(*,*) NF_STRERROR(ierr)4342 stop "getvarup"4343 endif4344 ! write(*,*)'lecture hadvq ok',hadvq4345 4346 #ifdef NC_DOUBLE4347 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),hadvu)4348 #else4349 ierr = NF_GET_VAR_REAL(nid,var3didin(19),hadvu)4350 #endif4351 if(ierr/=NF_NOERR) then4352 write(*,*) NF_STRERROR(ierr)4353 stop "getvarup"4354 endif4355 ! write(*,*)'lecture hadvu ok',hadvu4356 4357 #ifdef NC_DOUBLE4358 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),hadvv)4359 #else4360 ierr = NF_GET_VAR_REAL(nid,var3didin(20),hadvv)4361 #endif4362 if(ierr/=NF_NOERR) then4363 write(*,*) NF_STRERROR(ierr)4364 stop "getvarup"4365 endif4366 ! write(*,*)'lecture hadvv ok',hadvv4367 4368 #ifdef NC_DOUBLE4369 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(21),w)4370 #else4371 ierr = NF_GET_VAR_REAL(nid,var3didin(21),w)4372 #endif4373 if(ierr/=NF_NOERR) then4374 write(*,*) NF_STRERROR(ierr)4375 stop "getvarup"4376 endif4377 ! write(*,*)'lecture w ok',w4378 4379 #ifdef NC_DOUBLE4380 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(22),omega)4381 #else4382 ierr = NF_GET_VAR_REAL(nid,var3didin(22),omega)4383 #endif4384 if(ierr/=NF_NOERR) then4385 write(*,*) NF_STRERROR(ierr)4386 stop "getvarup"4387 endif4388 ! write(*,*)'lecture omega ok',omega4389 4390 return4391 end subroutine read_dice4392 !=====================================================================4393 subroutine read_gabls4(fich_gabls4,nlevel,ntime,nsol &4394 & ,zz,depth_sn,ug,vg,pf,th,t,qv,u,v,hadvt,hadvq,tg,tsnow,snow_dens)4395 4396 !program reading initial profils and forcings of the Gabls4 case study4397 4398 4399 implicit none4400 4401 #include "netcdf.inc"4402 4403 integer ntime,nlevel,nsol4404 integer l,k4405 character*80 :: fich_gabls44406 real*8 time(ntime)4407 4408 ! ATTENTION: visiblement quand on lit gabls4_driver.nc on recupere les donnees4409 ! dans un ordre inverse par rapport a la convention LMDZ4410 ! ==> il faut tout inverser (MPL 20141024)4411 ! les variables indexees "_i" sont celles qui sont lues dans gabls4_driver.nc4412 real*8 zz_i(nlevel),th_i(nlevel),pf_i(nlevel),t_i(nlevel)4413 real*8 qv_i(nlevel),u_i(nlevel),v_i(nlevel),ug_i(nlevel,ntime),vg_i(nlevel,ntime)4414 real*8 hadvt_i(nlevel,ntime),hadvq_i(nlevel,ntime)4415 4416 real*8 zz(nlevel),th(nlevel),pf(nlevel),t(nlevel)4417 real*8 qv(nlevel),u(nlevel),v(nlevel),ug(nlevel,ntime),vg(nlevel,ntime)4418 real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime)4419 4420 real*8 depth_sn(nsol),tsnow(nsol),snow_dens(nsol)4421 real*8 tg(ntime)4422 integer nid, ierr4423 integer nbvar3d4424 parameter(nbvar3d=30)4425 integer var3didin(nbvar3d)4426 4427 ierr = NF_OPEN(fich_gabls4,NF_NOWRITE,nid)4428 if (ierr.NE.NF_NOERR) then4429 write(*,*) 'ERROR: Pb opening forcings nc file '4430 write(*,*) NF_STRERROR(ierr)4431 stop ""4432 endif4433 4434 4435 ierr=NF_INQ_VARID(nid,"height",var3didin(1))4436 if(ierr/=NF_NOERR) then4437 write(*,*) NF_STRERROR(ierr)4438 stop 'height'4439 endif4440 4441 ierr=NF_INQ_VARID(nid,"depth_sn",var3didin(2))4442 if(ierr/=NF_NOERR) then4443 write(*,*) NF_STRERROR(ierr)4444 stop 'depth_sn'4445 endif4446 4447 ierr=NF_INQ_VARID(nid,"Ug",var3didin(3))4448 if(ierr/=NF_NOERR) then4449 write(*,*) NF_STRERROR(ierr)4450 stop 'Ug'4451 endif4452 4453 ierr=NF_INQ_VARID(nid,"Vg",var3didin(4))4454 if(ierr/=NF_NOERR) then4455 write(*,*) NF_STRERROR(ierr)4456 stop 'Vg'4457 endif4458 ierr=NF_INQ_VARID(nid,"pf",var3didin(5))4459 if(ierr/=NF_NOERR) then4460 write(*,*) NF_STRERROR(ierr)4461 stop 'pf'4462 endif4463 4464 ierr=NF_INQ_VARID(nid,"theta",var3didin(6))4465 if(ierr/=NF_NOERR) then4466 write(*,*) NF_STRERROR(ierr)4467 stop 'theta'4468 endif4469 4470 ierr=NF_INQ_VARID(nid,"tempe",var3didin(7))4471 if(ierr/=NF_NOERR) then4472 write(*,*) NF_STRERROR(ierr)4473 stop 'tempe'4474 endif4475 4476 ierr=NF_INQ_VARID(nid,"qv",var3didin(8))4477 if(ierr/=NF_NOERR) then4478 write(*,*) NF_STRERROR(ierr)4479 stop 'qv'4480 endif4481 4482 ierr=NF_INQ_VARID(nid,"u",var3didin(9))4483 if(ierr/=NF_NOERR) then4484 write(*,*) NF_STRERROR(ierr)4485 stop 'u'4486 endif4487 4488 ierr=NF_INQ_VARID(nid,"v",var3didin(10))4489 if(ierr/=NF_NOERR) then4490 write(*,*) NF_STRERROR(ierr)4491 stop 'v'4492 endif4493 4494 ierr=NF_INQ_VARID(nid,"hadvT",var3didin(11))4495 if(ierr/=NF_NOERR) then4496 write(*,*) NF_STRERROR(ierr)4497 stop 'hadvt'4498 endif4499 4500 ierr=NF_INQ_VARID(nid,"hadvQ",var3didin(12))4501 if(ierr/=NF_NOERR) then4502 write(*,*) NF_STRERROR(ierr)4503 stop 'hadvq'4504 endif4505 4506 ierr=NF_INQ_VARID(nid,"Tsnow",var3didin(14))4507 if(ierr/=NF_NOERR) then4508 write(*,*) NF_STRERROR(ierr)4509 stop 'tsnow'4510 endif4511 4512 ierr=NF_INQ_VARID(nid,"snow_density",var3didin(15))4513 if(ierr/=NF_NOERR) then4514 write(*,*) NF_STRERROR(ierr)4515 stop 'snow_density'4516 endif4517 4518 ierr=NF_INQ_VARID(nid,"Tg",var3didin(16))4519 if(ierr/=NF_NOERR) then4520 write(*,*) NF_STRERROR(ierr)4521 stop 'Tg'4522 endif4523 4524 4525 !dimensions lecture4526 ! call catchaxis(nid,ntime,nlevel,time,z,ierr)4527 4528 #ifdef NC_DOUBLE4529 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz_i)4530 #else4531 ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz_i)4532 #endif4533 if(ierr/=NF_NOERR) then4534 write(*,*) NF_STRERROR(ierr)4535 stop "getvarup"4536 endif4537 4538 #ifdef NC_DOUBLE4539 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),depth_sn)4540 #else4541 ierr = NF_GET_VAR_REAL(nid,var3didin(2),depth_sn)4542 #endif4543 if(ierr/=NF_NOERR) then4544 write(*,*) NF_STRERROR(ierr)4545 stop "getvarup"4546 endif4547 4548 #ifdef NC_DOUBLE4549 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),ug_i)4550 #else4551 ierr = NF_GET_VAR_REAL(nid,var3didin(3),ug_i)4552 #endif4553 if(ierr/=NF_NOERR) then4554 write(*,*) NF_STRERROR(ierr)4555 stop "getvarup"4556 endif4557 4558 #ifdef NC_DOUBLE4559 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),vg_i)4560 #else4561 ierr = NF_GET_VAR_REAL(nid,var3didin(4),vg_i)4562 #endif4563 if(ierr/=NF_NOERR) then4564 write(*,*) NF_STRERROR(ierr)4565 stop "getvarup"4566 endif4567 4568 #ifdef NC_DOUBLE4569 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),pf_i)4570 #else4571 ierr = NF_GET_VAR_REAL(nid,var3didin(5),pf_i)4572 #endif4573 if(ierr/=NF_NOERR) then4574 write(*,*) NF_STRERROR(ierr)4575 stop "getvarup"4576 endif4577 4578 #ifdef NC_DOUBLE4579 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),th_i)4580 #else4581 ierr = NF_GET_VAR_REAL(nid,var3didin(6),th_i)4582 #endif4583 if(ierr/=NF_NOERR) then4584 write(*,*) NF_STRERROR(ierr)4585 stop "getvarup"4586 endif4587 4588 #ifdef NC_DOUBLE4589 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),t_i)4590 #else4591 ierr = NF_GET_VAR_REAL(nid,var3didin(7),t_i)4592 #endif4593 if(ierr/=NF_NOERR) then4594 write(*,*) NF_STRERROR(ierr)4595 stop "getvarup"4596 endif4597 4598 #ifdef NC_DOUBLE4599 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),qv_i)4600 #else4601 ierr = NF_GET_VAR_REAL(nid,var3didin(8),qv_i)4602 #endif4603 if(ierr/=NF_NOERR) then4604 write(*,*) NF_STRERROR(ierr)4605 stop "getvarup"4606 endif4607 4608 #ifdef NC_DOUBLE4609 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),u_i)4610 #else4611 ierr = NF_GET_VAR_REAL(nid,var3didin(9),u_i)4612 #endif4613 if(ierr/=NF_NOERR) then4614 write(*,*) NF_STRERROR(ierr)4615 stop "getvarup"4616 endif4617 4618 #ifdef NC_DOUBLE4619 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),v_i)4620 #else4621 ierr = NF_GET_VAR_REAL(nid,var3didin(10),v_i)4622 #endif4623 if(ierr/=NF_NOERR) then4624 write(*,*) NF_STRERROR(ierr)4625 stop "getvarup"4626 endif4627 4628 #ifdef NC_DOUBLE4629 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),hadvt_i)4630 #else4631 ierr = NF_GET_VAR_REAL(nid,var3didin(11),hadvt_i)4632 #endif4633 if(ierr/=NF_NOERR) then4634 write(*,*) NF_STRERROR(ierr)4635 stop "getvarup"4636 endif4637 4638 #ifdef NC_DOUBLE4639 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),hadvq_i)4640 #else4641 ierr = NF_GET_VAR_REAL(nid,var3didin(12),hadvq_i)4642 #endif4643 if(ierr/=NF_NOERR) then4644 write(*,*) NF_STRERROR(ierr)4645 stop "getvarup"4646 endif4647 4648 #ifdef NC_DOUBLE4649 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),tsnow)4650 #else4651 ierr = NF_GET_VAR_REAL(nid,var3didin(14),tsnow)4652 #endif4653 if(ierr/=NF_NOERR) then4654 write(*,*) NF_STRERROR(ierr)4655 stop "getvarup"4656 endif4657 4658 #ifdef NC_DOUBLE4659 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),snow_dens)4660 #else4661 ierr = NF_GET_VAR_REAL(nid,var3didin(15),snow_dens)4662 #endif4663 if(ierr/=NF_NOERR) then4664 write(*,*) NF_STRERROR(ierr)4665 stop "getvarup"4666 endif4667 4668 #ifdef NC_DOUBLE4669 ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),tg)4670 #else4671 ierr = NF_GET_VAR_REAL(nid,var3didin(16),tg)4672 #endif4673 if(ierr/=NF_NOERR) then4674 write(*,*) NF_STRERROR(ierr)4675 stop "getvarup"4676 endif4677 4678 ! On remet les variables lues dans le bon ordre des niveaux (MPL 20141024)4679 do k=1,nlevel4680 zz(k)=zz_i(nlevel+1-k)4681 ug(k,:)=ug_i(nlevel+1-k,:)4682 vg(k,:)=vg_i(nlevel+1-k,:)4683 pf(k)=pf_i(nlevel+1-k)4684 print *,'pf=',pf(k)4685 th(k)=th_i(nlevel+1-k)4686 t(k)=t_i(nlevel+1-k)4687 qv(k)=qv_i(nlevel+1-k)4688 u(k)=u_i(nlevel+1-k)4689 v(k)=v_i(nlevel+1-k)4690 hadvt(k,:)=hadvt_i(nlevel+1-k,:)4691 hadvq(k,:)=hadvq_i(nlevel+1-k,:)4692 enddo4693 return4694 end subroutine read_gabls44695 !=====================================================================4696 4697 ! Reads CIRC input files4698 4699 SUBROUTINE read_circ(nlev_circ,cf,lwp,iwp,reliq,reice,t,z,p,pm,h2o,o3,sza)4700 4701 parameter (ncm_1=49180)4702 #include "YOMCST.h"4703 4704 real albsfc(ncm_1), albsfc_w(ncm_1)4705 real cf(nlev_circ), icefra(nlev_circ), deice(nlev_circ), &4706 reliq(nlev_circ), reice(nlev_circ), lwp(nlev_circ), iwp(nlev_circ)4707 real t(nlev_circ+1), z(nlev_circ+1), dz(nlev_circ), p(nlev_circ+1)4708 real aer_beta(nlev_circ), waer(nlev_circ), gaer(nlev_circ)4709 real pm(nlev_circ), tm(nlev_circ), h2o(nlev_circ), o3(nlev_circ)4710 real co2(nlev_circ), n2o(nlev_circ), co(nlev_circ), ch4(nlev_circ), &4711 o2(nlev_circ), ccl4(nlev_circ), f11(nlev_circ), f12(nlev_circ)4712 ! za= zenital angle4713 ! sza= cosinus angle zenital4714 real wavn(ncm_1), ssf(ncm_1),za,sza4715 integer nlev4716 4717 4718 ! Open the files4719 4720 open (11, file='Tsfc_sza_nlev_case.txt', status='old')4721 open (12, file='level_input_case.txt', status='old')4722 open (13, file='layer_input_case.txt', status='old')4723 open (14, file='aerosol_input_case.txt', status='old')4724 open (15, file='cloud_input_case.txt', status='old')4725 open (16, file='sfcalbedo_input_case.txt', status='old')4726 4727 ! Read scalar information4728 do iskip=1,54729 read (11, *)4730 enddo4731 read (11, '(i8)') nlev4732 read (11, '(f10.2)') tsfc4733 read (11, '(f10.2)') za4734 read (11, '(f10.4)') sw_dn_toa4735 sza=cos(za/180.*RPI)4736 print *,'nlev,tsfc,sza,sw_dn_toa,RPI',nlev,tsfc,sza,sw_dn_toa,RPI4737 close(11)4738 4739 ! Read level information4740 read (12, *)4741 do il=1,nlev4742 read (12, 302) ilev, z(il), p(il), t(il)4743 z(il)=z(il)*1000. ! z donne en km4744 p(il)=p(il)*100. ! p donne en mb4745 enddo4746 302 format (i8, f8.3, 2f9.2)4747 close(12)4748 4749 ! Read layer information (midpoint values)4750 do iskip=1,34751 read (13, *)4752 enddo4753 do il=1,nlev-14754 read (13, 303) ilev,pm(il),tm(il),h2o(il),co2(il),o3(il), &4755 n2o(il),co(il),ch4(il),o2(il),ccl4(il), &4756 f11(il),f12(il)4757 pm(il)=pm(il)*100.4758 enddo4759 303 format (i8, 2f9.2, 10(2x,e13.7))4760 close(13)4761 4762 ! Read aerosol layer information4763 do iskip=1,34764 read (14, *)4765 enddo4766 read (14, '(f10.2)') aer_alpha4767 read (14, *)4768 read (14, *)4769 do il=1,nlev-14770 read (14, 304) ilev, aer_beta(il), waer(il), gaer(il)4771 enddo4772 304 format (i8, f9.5, 2f8.3)4773 close(14)4774 4775 ! Read cloud information4776 do iskip=1,34777 read (15, *)4778 enddo4779 do il=1,nlev-14780 read (15, 305) ilev, cf(il), lwp(il), iwp(il), reliq(il), reice(il)4781 lwp(il)=lwp(il)/1000. ! lwp donne en g/kg4782 iwp(il)=iwp(il)/1000. ! iwp donne en g/kg4783 reliq(il)=reliq(il)/1000000. ! reliq donne en microns4784 reice(il)=reice(il)/1000000. ! reice donne en microns4785 enddo4786 305 format (i8, f8.3, 4f9.2)4787 close(15)4788 4789 ! Read surface albedo (weighted & unweighted) and spectral solar irradiance4790 do iskip=1,64791 read (16, *)4792 enddo4793 do icm_1=1,ncm_14794 read (16, 306) wavn(icm_1), albsfc(icm_1), albsfc_w(icm_1), ssf(icm_1)4795 enddo4796 306 format(f10.1, 2f12.5, f14.8)4797 close(16)4798 4799 return4800 end subroutine read_circ4801 !=====================================================================4802 ! Reads RTMIP input files4803 4804 SUBROUTINE read_rtmip(nlev_rtmip,play,plev,t,h2o,o3)4805 4806 #include "YOMCST.h"4807 4808 real t(nlev_rtmip), pt(nlev_rtmip),pb(nlev_rtmip),h2o(nlev_rtmip), o3(nlev_rtmip)4809 real temp(nlev_rtmip), play(nlev_rtmip),ovap(nlev_rtmip), oz(nlev_rtmip),plev(nlev_rtmip+1)4810 integer nlev4811 4812 4813 ! Open the files4814 4815 open (11, file='low_resolution_profile.txt', status='old')4816 4817 ! Read level information4818 read (11, *)4819 do il=1,nlev_rtmip4820 read (11, 302) pt(il), pb(il), t(il),h2o(il),o3(il)4821 enddo4822 do il=1,nlev_rtmip4823 play(il)=pt(nlev_rtmip-il+1)*100. ! p donne en mb4824 temp(il)=t(nlev_rtmip-il+1)4825 ovap(il)=h2o(nlev_rtmip-il+1)4826 oz(il)=o3(nlev_rtmip-il+1)4827 enddo4828 do il=1,394829 plev(il)=play(il)+(play(il+1)-play(il))/2.4830 print *,'il p t ovap oz=',il,plev(il),temp(il),ovap(il),oz(il)4831 enddo4832 plev(41)=101300.4833 302 format (e16.10,3x,e16.10,3x,e16.10,3x,e12.6,3x,e12.6)4834 close(12)4835 4836 return4837 end subroutine read_rtmip4838 !=====================================================================4839 1473 4840 1474 ! Subroutines for nudging … … 5125 1759 real frac,frac1,frac2,fact 5126 1760 5127 do l = 1, llm5128 print *,'debut interp2, play=',l,play(l)5129 enddo1761 ! do l = 1, llm 1762 ! print *,'debut interp2, play=',l,play(l) 1763 ! enddo 5130 1764 ! do l = 1, nlev_cas 5131 1765 ! print *,'debut interp2, plev_prof_cas=',l,play(l),plev_prof_cas(l) … … 5137 1771 5138 1772 mxcalc=l 5139 print *,'debut interp2, mxcalc=',mxcalc1773 ! print *,'debut interp2, mxcalc=',mxcalc 5140 1774 k1=0 5141 1775 k2=0 -
Property
svn:keywords
set to
Note: See TracChangeset
for help on using the changeset viewer.
