{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_s1982-2016\n",
      "ACC10mem_fgco2_HIN1_s1982-2016_LY1_5x5 MSSS10mem_fgco2_HIN1_s1982-2016_LY1_5x5\n",
      "(35, 10, 36, 72)\n",
      "(35, 36, 72)\n",
      "p_fdr =  0.08225000000000006\n",
      "_s1982-2016\n",
      "ACC10mem_fgco2_HIN1-ANA1_s1982-2016_LY1_5x5 MSSS10mem_fgco2_HIN1-ANA1_s1982-2016_LY1_5x5\n",
      "(35, 10, 36, 72)\n",
      "(35, 10, 36, 72)\n",
      "(35, 36, 72)\n",
      "p_fdr =  0.03875000000000003\n",
      "_s1982-2016\n",
      "ACC10mem_fgco2_HIN2-HIN1_s1982-2016_LY1_5x5 MSSS10mem_fgco2_HIN2-HIN1_s1982-2016_LY1_5x5\n",
      "(35, 10, 36, 72)\n",
      "(35, 10, 36, 72)\n",
      "(35, 36, 72)\n",
      "p_fdr =  0.0\n",
      "_s1982-2016\n",
      "ACC10mem_fgco2_HIN1-HIST_s1982-2016_LY1_5x5 MSSS10mem_fgco2_HIN1-HIST_s1982-2016_LY1_5x5\n",
      "(35, 10, 36, 72)\n",
      "(35, 10, 36, 72)\n",
      "(35, 36, 72)\n",
      "p_fdr =  0.0\n",
      "_s1983-2016\n",
      "ACC10mem_fgco2_HIN1-PERS_s1983-2016_LY1_5x5 MSSS10mem_fgco2_HIN1-PERS_s1983-2016_LY1_5x5\n",
      "(34, 10, 36, 72)\n",
      "(34, 36, 72)\n",
      "(34, 36, 72)\n",
      "p_fdr =  0.17750000000000013\n",
      "_s1982-2012\n",
      "ACC10mem_fgco2_HIN1_s1982-2012_LY2-5_5x5 MSSS10mem_fgco2_HIN1_s1982-2012_LY2-5_5x5\n",
      "(31, 10, 36, 72)\n",
      "(31, 36, 72)\n",
      "p_fdr =  0.1\n",
      "_s1982-2012\n",
      "ACC10mem_fgco2_HIN1-ANA1_s1982-2012_LY2-5_5x5 MSSS10mem_fgco2_HIN1-ANA1_s1982-2012_LY2-5_5x5\n",
      "(31, 10, 36, 72)\n",
      "(31, 10, 36, 72)\n",
      "(31, 36, 72)\n",
      "p_fdr =  0.022500000000000017\n",
      "_s1982-2012\n",
      "ACC10mem_fgco2_HIN2-HIN1_s1982-2012_LY2-5_5x5 MSSS10mem_fgco2_HIN2-HIN1_s1982-2012_LY2-5_5x5\n",
      "(31, 10, 36, 72)\n",
      "(31, 10, 36, 72)\n",
      "(31, 36, 72)\n",
      "p_fdr =  0.0\n",
      "_s1982-2012\n",
      "ACC10mem_fgco2_HIN1-HIST_s1982-2012_LY2-5_5x5 MSSS10mem_fgco2_HIN1-HIST_s1982-2012_LY2-5_5x5\n",
      "(31, 10, 36, 72)\n",
      "(31, 10, 36, 72)\n",
      "(31, 36, 72)\n",
      "p_fdr =  0.0\n",
      "_s1986-2012\n",
      "ACC10mem_fgco2_HIN1-PERS_s1986-2012_LY2-5_5x5 MSSS10mem_fgco2_HIN1-PERS_s1986-2012_LY2-5_5x5\n",
      "(27, 10, 36, 72)\n",
      "(27, 36, 72)\n",
      "(27, 36, 72)\n",
      "p_fdr =  0.10875000000000008\n",
      "_s1982-2008\n",
      "ACC10mem_fgco2_HIN1_s1982-2008_LY6-9_5x5 MSSS10mem_fgco2_HIN1_s1982-2008_LY6-9_5x5\n",
      "(27, 10, 36, 72)\n",
      "(27, 36, 72)\n",
      "p_fdr =  0.1\n",
      "_s1982-2008\n",
      "ACC10mem_fgco2_HIN1-ANA1_s1982-2008_LY6-9_5x5 MSSS10mem_fgco2_HIN1-ANA1_s1982-2008_LY6-9_5x5\n",
      "(27, 10, 36, 72)\n",
      "(27, 10, 36, 72)\n",
      "(27, 36, 72)\n",
      "p_fdr =  0.023500000000000017\n",
      "_s1982-2008\n",
      "ACC10mem_fgco2_HIN2-HIN1_s1982-2008_LY6-9_5x5 MSSS10mem_fgco2_HIN2-HIN1_s1982-2008_LY6-9_5x5\n",
      "(27, 10, 36, 72)\n",
      "(27, 10, 36, 72)\n",
      "(27, 36, 72)\n",
      "p_fdr =  0.0\n",
      "_s1982-2008\n",
      "ACC10mem_fgco2_HIN1-HIST_s1982-2008_LY6-9_5x5 MSSS10mem_fgco2_HIN1-HIST_s1982-2008_LY6-9_5x5\n",
      "(27, 10, 36, 72)\n",
      "(27, 10, 36, 72)\n",
      "(27, 36, 72)\n",
      "p_fdr =  0.0\n",
      "_s1986-2008\n",
      "ACC10mem_fgco2_HIN1-PERS_s1986-2008_LY6-9_5x5 MSSS10mem_fgco2_HIN1-PERS_s1986-2008_LY6-9_5x5\n",
      "(23, 10, 36, 72)\n",
      "(23, 36, 72)\n",
      "(23, 36, 72)\n",
      "p_fdr =  0.08800000000000006\n"
     ]
    }
   ],
   "source": [
    "import sys ; sys.path.remove('/mnt/bcpu-ns9039k/ingo/jupyter/Modules'); sys.path.append('../../scripts')\n",
    "import norcpmTools as nt\n",
    "from netCDF4 import Dataset\n",
    "import numpy as np\n",
    "\n",
    "# prepare fgco2 data\n",
    "if False:\n",
    "    field='fgco2'\n",
    "    res=[5,5]\n",
    "    for experiment in ['historical']:\n",
    "        nt.writeEns(field,experiment,[1950,2029],[1,10])\n",
    "        nt.regrid2reg(field,experiment,[1950,2029],res=res,memRange=[1,10],ensave=False)\n",
    "    for experiment in ['dcppA-assim-i1','dcppA-assim-i2']:\n",
    "        nt.writeEns(field,experiment,[1950,2018],[1,10])\n",
    "        nt.regrid2reg(field,experiment,[1950,2018],res=res,memRange=[1,10],ensave=False)\n",
    "    for experiment in ['dcppA-hindcast-i1','dcppA-hindcast-i2']:\n",
    "        for syear in range(1960,2019):\n",
    "            year1 = syear + 1\n",
    "            yearn = syear + 10\n",
    "            nt.writeEns(field,experiment,[year1,yearn],[1,10],syear=syear)\n",
    "            nt.regrid2reg(field,experiment,[year1,yearn],res=res,memRange=[1,10],syear=syear,ensave=False,isMasked=False)\n",
    "\n",
    "# plot ACC \n",
    "fields = [['fgco2','fgco2']]\n",
    "resOptions=[[5,5]]\n",
    "leadRanges = [[0,0],[1,4],[5,8]]\n",
    "expOptions = [['HIN1',''],['HIN1','ANA1'],['HIN2','HIN1'],['HIN1','HIST'],['HIN1','PERS']]\n",
    "memRange = [1,10] \n",
    "doACC = True\n",
    "doMSSS = False\n",
    "product = 'hindcast'\n",
    "for field in fields:\n",
    "    fieldMod = field[0]\n",
    "    fieldObs = field[1]\n",
    "    if fieldMod == 'fgco2':\n",
    "        fieldObs = 'fgco2'\n",
    "        obsName = 'SOCCOM'\n",
    "        obsCoverage = [1982,2017]\n",
    "        obsFactor = -12/(1000*365*24*3600) # mod = kg C m-2 s-1 , obs = mol/m2/yr\n",
    "        obsOffset = 0.\n",
    "        levRange = [0,0]\n",
    "        landFill = True        \n",
    "    for leadRange in leadRanges:   \n",
    "        for res in resOptions:\n",
    "            lon = np.arange(res[0]/2,360,res[0])\n",
    "            lat = np.arange(-90+res[1]/2,90,res[1])\n",
    "            lon2, lat2 = np.meshgrid(lon,lat)\n",
    "            mskfdr = np.where(lat2 < 80, 1, 0)\n",
    "            tagRes = '_{:d}x{:d}'.format(res[0],res[1])\n",
    "            tagField = '_' + fieldMod if obsName[0:11] != 'dcppA-assim' else '_' + fieldMod + 'Perfect'\n",
    "            tagLead = '_LY{:d}'.format(leadRange[0]+1) if leadRange[0] == leadRange[1] else '_LY{:d}-{:d}'.format(leadRange[0]+1,leadRange[1]+1)\n",
    "            # extract data\n",
    "            for expOption in expOptions:\n",
    "                if expOption[0][0:3] == 'ANA' or expOption[1][0:3] == 'ANA':\n",
    "                    modCoverage = [1950,2018]\n",
    "                else:\n",
    "                    modCoverage = [1950,2029]\n",
    "                if product == 'analysis':    \n",
    "                    syear1 = np.max((modCoverage[0],obsCoverage[0]))\n",
    "                else:\n",
    "                    syear1 = np.max((1960,obsCoverage[0])) if not expOption[1] == 'PERS' else np.max((1960,obsCoverage[0]+1+leadRange[1]-leadRange[0]))\n",
    "                syearn = np.min((modCoverage[1],obsCoverage[1]))-leadRange[1]-1\n",
    "                syears = range(syear1,syearn+1)\n",
    "                tagYears = '_s{:d}-{:d}'.format(syears[0],syears[-1])\n",
    "                print(tagYears)\n",
    "                tagYears = '_s{:d}-{:d}'.format(syears[0],syears[-1])\n",
    "                if obsName == 'GlobColour':\n",
    "                    obs = np.flip(nt.readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,levRange=levRange,suffix=tagRes),axis=1)                \n",
    "                elif obsName[0:11] == 'dcppA-assim': \n",
    "                    obs = np.mean(nt.readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,memRange=[1,10],levRange=levRange,suffix=tagRes,ensave=False),axis=1)\n",
    "                else:\n",
    "                    obs = nt.readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,levRange=levRange,suffix=tagRes)     \n",
    "                if fieldObs == 'fgco2':\n",
    "                    obs = -obs\n",
    "                if expOption[0] == 'HIST':\n",
    "                    fld1 =  nt.readHindcastLY(fieldMod,'historical',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False)\n",
    "                elif expOption[0] == 'PERS':\n",
    "                    if obsName[0:11] == 'dcppA-assim':\n",
    "                        fld1 = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,memRange=[1,10],\n",
    "                                              levRange=levRange,suffix=tagRes,persistence='mean',ensave=False)  \n",
    "                    elif obsName == 'GlobColour':\n",
    "                        fld1 = np.flip(readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,\n",
    "                                                      levRange=levRange,suffix=tagRes,persistence='mean',ensave=False),axis=1)\n",
    "                    else:\n",
    "                        fld1 = readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,\n",
    "                                              levRange=levRange,suffix=tagRes,persistence='mean',ensave=False)\n",
    "                elif expOption[0] == 'HIN1':\n",
    "                    fld1 = nt.readHindcastLY(fieldMod,'dcppA-hindcast-i1',syears,leadRange,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False)\n",
    "                elif expOption[0] == 'HIN2':\n",
    "                    fld1 = nt.readHindcastLY(fieldMod,'dcppA-hindcast-i2',syears,leadRange,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False)\n",
    "                elif expOption[0] == 'ANA1':\n",
    "                    fld1 = np.squeeze(nt.readHindcastLY(fieldMod,'dcppA-assim-i1',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False))\n",
    "                elif expOption[0] == 'ANA2':\n",
    "                    fld1 = nt.readHindcastLY(fieldMod,'dcppA-assim-i2',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False)               \n",
    "                #\n",
    "                if expOption[1] == 'HIST':\n",
    "                    fld2 =  nt.readHindcastLY(fieldMod,'historical',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False)\n",
    "                elif expOption[1] == 'PERS':\n",
    "                    fld2 = nt.readHindcastLY(fieldObs,obsName,syears,leadRange,yearRange=obsCoverage,levRange=levRange,suffix=tagRes,persistence='mean',ensave=False)\n",
    "                elif expOption[1] == 'HIN1':\n",
    "                    fld2 = nt.readHindcastLY(fieldMod,'dcppA-hindcast-i1',syears,leadRange,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False)\n",
    "                elif expOption[1] == 'HIN2':\n",
    "                    fld2 = nt.readHindcastLY(fieldMod,'dcppA-hindcast-i2',syears,leadRange,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False)\n",
    "                elif expOption[1] == 'ANA1':\n",
    "                    fld2 = np.squeeze(nt.readHindcastLY(fieldMod,'dcppA-assim-i1',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False))\n",
    "                elif expOption[1] == 'ANA2':\n",
    "                    fld2 = nt.readHindcastLY(fieldMod,'dcppA-assim-i2',syears,leadRange,yearRange=modCoverage,memRange=memRange,levRange=levRange,suffix=tagRes,ensave=False)  \n",
    "                #\n",
    "                tagExp = '_' + expOption[0] if expOption[1] == '' else '_{:s}-{:s}'.format(expOption[0],expOption[1])\n",
    "                rfilePrefix = 'ACC10mem' + tagField + tagExp + tagYears + tagLead + tagRes \n",
    "                mfilePrefix = 'MSSS10mem' + tagField + tagExp + tagYears + tagLead + tagRes \n",
    "                print(rfilePrefix + ' ' + mfilePrefix)\n",
    "                #\n",
    "                print(fld1.shape)\n",
    "                if expOption[1] != '':\n",
    "                    print(fld2.shape)\n",
    "                print(obs.shape)\n",
    "                titleString = expOption[0] if expOption[1] == '' else expOption[0] + ' - ' + expOption[1]\n",
    "                title2 = tagLead[1:] if product == 'hindcast' else ''                    \n",
    "                if doACC: \n",
    "                    fld = nt.corrMultiArrayYeager(fld1,obs) if expOption[1] == '' else nt.corrMultiArrayDiffYeager(fld1,obs,fld2)\n",
    "                    nt.plotACC(lon=lon,lat=lat,fld=fld,filePrefix=rfilePrefix,lbLabelBarOn=False,\n",
    "                               title=' ',title2=' ',landFill=landFill,plottype='ACC')\n",
    "                if doMSSS: \n",
    "                    fld = nt.MSSSyeager(fld1,obs) if expOption[1] == '' else nt.MSSSyeager(fld1,obs,fld2)\n",
    "                    nt.plotACC(lon=lon,lat=lat,fld=fld,filePrefix=mfilePrefix,lbLabelBarOn=False,\n",
    "                               title=titleString,title2=title2,landFill=landFill,plottype='MSSS')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/mnt/bcpu-ns9039k/ingo/NorCPM1-paper/revision1/fig17\n"
     ]
    }
   ],
   "source": [
    "cd fig17\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
