;ZTS Aug 2023 -- Onset prediction for daily forecasts from C3S data store ;ZTS: 8 Aug 2023 -- CFSv2 added load "./calculateRunLenTimePosition.ncl" load "./CessationFun.ncl" begin ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;; CHANGEABLE PARAMETERS ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; modelNames = (/"eccc","ecmwf","metfr","cmcc","cfs"/) version = (/3,51,8,35,2/) ensmax = 52 ; ;fDatDir = "./precip/" nmodrunS = MODELSTRT nmodrunE = MODELLAST nmodels = nmodrunE - nmodrunS + 1 initYear = CURRENTYR startmon = OMONS ; start month for onset/cessation startday = ODAYS ; start date of onset seasrch endmon = OMONE ; last month for onset search endday = ODAYE ; last day for onset search SEASON = SEASONNAME hr = 00 ; daily forecasts are initialized from 00Z for C3S daily sMonOnset = ONSETSTART ; actual start month of search for onset ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; THESE ARE IMPORTANT ONSET/CESSATION CRITERIA. CHANGE FOR YOUR CLIMATE ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; heavyrain = RHEAVY ;mm/day ; CHANGE AS NEEDED rainzero = RZERO ;dry and wet rainy day threshold rtrsholdstart = RTHROLD ;total rainfall in numraindays for onset numraindays = NRD ;Num. of days in which total rain is computed ndayafteronset = NDAYAFTER ;total days DRY SPELL is calculated ndayrain = NDAYR ;lowest # of rainy days for onset limnorainydayafter = NORAINAFTER ;maximum len of drysp after onset in 21 days cessndrylen = NCESSN ;lenth of dry days for cessn after onset + 20 cessanrainamount = RCESSN ;total rainfall to the end of december 31 latS = LATS latN = LATN lonW = LONW lonE = LONE Res = RESL ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;; END OF CHANGEABLE PARAMETERS ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if (lonW .lt. 0.) then lonW = lonW + 360. end if if (lonE .lt. 0.) then lonE = lonE + 360. end if tmp := tochar(tostring(Res)) fDatDir = "./precip/regionalp" + tmp(2) + "/" analdir = "./analdat/" system("if ! test -d " + analdir +" ; then mkdir -p " + analdir + " ; fi") monthdef = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep", \ "Oct","Nov","Dec"/) initMonth = monthdef(startmon) writeoutany = True ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; this assignment use lat-lon defined earlier near the top mlat = toint((latN-latS)/Res + 1) mlon = toint((lonE-lonW)/Res + 1) REGlat = fspan(latS, latN, mlat) REGlon = fspan(lonW, lonE, mlon) REGlat@units = "degrees_north" REGlon@units = "degrees_east" REGlat!0 = "lat" REGlon!0 = "lon" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;; Define hindcast time variable to save computed values yrs = ispan(1990,2024,1) mm = yrs mm(:) = sMonOnset dd = mm dd(:) = startday hh = mm hh(:) = 00 mn = hh ss = hh units = "days since 1990-01-01 00:00" TIME = cd_inv_calendar(yrs,mm,dd,hh,mn,ss,units,0) TIME!0 = "time" nyears = dimsizes(TIME) time = cd_calendar(TIME,-3) time!0 = "time" time@units = "yyyymmddhh" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; prepare variables ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ensall = ispan(0,ensmax - 1, 1) ensall!0 = "ens" ensall&ens = ensall ensall@description = "nens=ensize + 1; 0 is for ens mean, 1:nens for members" jlatStrt = 0 jlatLast = mlat - 1 ilonStrt = 0 ilonLast = mlon - 1 ;;;; save onset for the ensmean and individual runs precip ; Note nyears is general: dimsizes(TIME); not model hindcast yrs onsetAll = new((/nyears, nmodels, ensmax, mlat,mlon/),float,-9999) onsetAll!0 = "time" onsetAll!1 = "models" onsetAll!2 = "ens" onsetAll!3 = "lat" onsetAll!4 = "lon" varx = "" do i= nmodrunS, nmodrunE varx = varx + " " + modelNames(i) end do models = ispan(0,nmodels - 1,1) models!0 = "models" models&models = models models@modelnames = varx onsetAll&time = time onsetAll&models = models onsetAll&ens = ensall ;ensmax ;ispan(0,ensmax-1,1);13/Jan/2022: 0 for ensemble mean, ; 1:enssize is for indivisual members. Since we plan to do ; probability based forecasts, we will have to use ; only ensemble mean onsetAll&lat = REGlat onsetAll&lon = REGlon cesationAll = onsetAll growlengAll = onsetAll ;drystartAll = onsetAll ; drystlenAll = onsetAll dryspLengFix = onsetAll ; max spells for startmon/day to endmon/day dryspLoc0Fix = onsetAll ; the first time location of max dryspell dryspLoc1Fix = onsetAll ; the second location if any of max dryspell dryspLengAct = onsetAll ; max spells subsequent to calculated onset dryspLoc0Act = onsetAll ; the first time location of max dryspell dryspLoc1Act = onsetAll ; the second time location of max dryspell, if any wetspLengFix = onsetAll ; max spells for startmon/day to endmon/day wetspLoc0Fix = onsetAll ; the first time location of max wetspell wetspLoc1Fix = onsetAll ; the second time location of max wetspell, if any wetspLengAct = onsetAll ; max spell subsequent to calculated onset wetspLoc0Act = onsetAll ; the first time location of max wetspell wetspLoc1Act = onsetAll ; the second time loc of max wet spell, if any ;;;;;;;; ;;;; assign definition and attributes onsetAll@long_name = "Onset date for " + SEASON cesationAll@long_name = "Cessation date for " + SEASON growlengAll@long_name = "Growing Length for " + SEASON dryspLengFix@long_name = "Max dryspell length between " + startday+"-"+monthdef(startmon)+ " to " + endday+"-"+endmon dryspLoc0Fix@long_name = "Time location of the 1st max drysp occuring during " + startday+"-"+monthdef(startmon)+ " to " + endday+"-"+endmon dryspLoc1Fix@long_name = "Time location of the 2nd max drysp occuring during " + startday+"-"+monthdef(startmon)+ " to " + endday+"-"+endmon dryspLengAct@long_name = "Max dryspell length within predicted onset and cessation" dryspLoc0Act@long_name = "Time location of the 1st max drysp within predicted onset and cessation" dryspLoc1Act@long_name = "Time location of the 2nd max drysp within predicted onset and cessation" wetspLengFix@long_name = "Max wetspell length between " + startday+"-"+monthdef(startmon)+ " to " + endday+"-"+endmon wetspLoc0Fix@long_name = "Time location of the 1st max wetsp occuring during " + startday+"-"+monthdef(startmon)+ " to " + endday+"-"+endmon wetspLoc1Fix@long_name = "Time location of the 2nd max wetsp occuring during " + startday+"-"+monthdef(startmon)+ " to " + endday+"-"+endmon wetspLengAct@long_name = "Max wetspell length within predicted onset and cessation" wetspLoc0Act@long_name = "Time location of the 1st max wetsp within predicted onset and cessation" wetspLoc1Act@long_name = "Time location of the 2nd max wetsp within predicted onset and cessation" ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;; monthly dry/wet spells ;;;;;;;;;;;;;;;;;;; mondryS = startmon mondryE = endmon nummon = mondryE - mondryS + 1 ; number of month dmonth = ispan(0,nummon - 1,1) dmonth!0 = "month" dmonth@start_month = mondryS dmonth@last_month = mondryE dryMonAll = new((/nyears,nummon, nmodels, ensmax, mlat,mlon/),float,-9999) dryMonAll!0 = "time" dryMonAll!1 = "month" dryMonAll!2 = "models" dryMonAll!3 = "ens" dryMonAll!4 = "lat" dryMonAll!5 = "lon" ;ensmax ;ispan(0,ensmax-1,1); 0 for ensemble mean, ; 1:enssize is for indivisual members. Since we plan to do ; probability based forecasts, we will have to use ; only ensemble mean dryMonAll&time = time dryMonAll&month = dmonth ;ispan(0,nummon - 1,1) dryMonAll&models = models dryMonAll&ens = ensall dryMonAll&lat = REGlat dryMonAll&lon = REGlon dryMonAll@long_name = "Maximum dry spell in a month " wetMonAll = dryMonAll wetMonAll@long_name = "Maximum wet spell in a month " ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Calculate onset for each model separately ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; do imodel = nmodrunS, nmodrunE MODEL = str_upper(modelNames(imodel))+ "v" +version(imodel) fname = fDatDir + "precDaily"+ initMonth+"Hindcast" + MODEL +".nc" a := addfile(fname, "r") tmp := a->prec ; (year,lead/day,ens,lat,lon) rain_tot := tmp(:,:,:,{latS:latN},{lonW:lonE}) ztemp := a->time ; yyyymmddhh; initial time in the Hindcast variable leadAll := a->lead ; Offset (in days) from each initial time ens := a->ens enssize = dimsizes(ens) nens = enssize + 1 ens := ispan(0, enssize,1) ens!0 = "ens" ens&ens = ens ens@description = "nens=ensize + 1; 0 is for ens mean, 1:nens for members" ;;;;;;;; intime := ztemp ;cd_calendar(ztemp,-3) ;yyyymmddhh initial time in the Hindcast variable hindYrs := intime/1000000 ;This is the actual years in a model mmddhh := intime - hindYrs*1000000 allmons := mmddhh/10000 ddhh := mmddhh - allmons*10000 allday := ddhh/100 allhrs := ddhh - allday*100 ;;;;;;;;;;;;;;; ;;;;;;; Loop over all years (initial time) ;;;;;; do iyr = 0, dimsizes(hindYrs) - 1 iyear = hindYrs(iyr) yrind = ind(yrs .eq. iyear ) ; this is yr index to save in variable ; yrs contains all possible model years imonth = allmons(iyr) iday = allday(iyr) ihr = allhrs(iyr) lead := leadAll ;a->lead--Offset (in days) from each initial time lead!0 = "time" lead&time = leadAll x1 = "days since "+iyear+ "-" + sprinti("%0.2i",imonth) x2 = "-" + sprinti("%0.2i",iday) + " " x3 = sprinti("%0.2i",ihr) + ":00:0" lead@units = x1 + x2 + x3 lead = (/ lead -1 /) ;rain is measured for last 24 hrs ending today utc_date := cd_calendar(lead, 0) ;gives 6 col for yr,MM,DD,HH,mm,ss fmonth := tointeger(utc_date(:,1)) ; all months of data for each day fday := tointeger(utc_date(:,2)) ; all days of data for a given year timelabelsDat := monthdef(fmonth)+"-"+fday ; full forecast time ;tmpS = monthdef(startmon)+ "-" + startday ; onset search start day tmpS = monthdef(sMonOnset)+ "-" + startday ; onset search start day tmpE = monthdef(endmon) + "-" + endday ; onset search last day print(MODEL+" with data dates from "+ timelabelsDat(0)+ " to " +timelabelsDat(dimsizes(timelabelsDat)-1)) print(MODEL + " Onset search is from " + tmpS +" to " + tmpE) indxOnsetStart = get1Dindex(timelabelsDat,tmpS) indxOnsetEnd = get1Dindex(timelabelsDat,tmpE) ; print(indxOnsetStart+ " " + indxOnsetEnd) timelabels := timelabelsDat(indxOnsetStart:indxOnsetEnd) ; print(dimsizes(timelabels)) delete([/fday,fmonth, tmpS, tmpE/]) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Create a daily rainfall variable; set it to 0 if less than 0.01mm ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; do iens= 0, enssize ; 0 save for mean, all others for members startindx = indxOnsetStart ; The onset computation in model starts ; from startmon/startday endindx = indxOnsetEnd ; Onset search ends at the specified date ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;; Now go through all ensemble to determine onset stuff ;;;;;;;;;;;;;;;;;;;;;;;;;; if (iens .gt. 0) then raindaily := rain_tot(iyr,startindx:endindx,iens-1,:,:) rainSpell := rain_tot(iyr,:,iens-1,:,:) ; print(dimsizes(raindaily)+" "+"should be nday x nlat x nlon" ) else raindaily := dim_avg_n_Wrap(rain_tot(iyr,startindx:endindx,:,:,:),1) rainSpell := dim_avg_n_Wrap(rain_tot(iyr,:,:,:,:),1) ;print(dimsizes(raindaily)+" "+"should be nday x nlat x nlon" ) end if ;;;;;;;;;;;;;;; ;;;; Compute monthly dry and wet spells ;;;;;;;;;;;;;;; mymonth = 0 do imond = startmon, endmon ; nummon - 1 iendday = days_in_month(iyear, imond) myS = monthdef(imond) + "-" + 1 ; day of imond myE = monthdef(imond) + "-" + iendday ; end day of imond istartindx = get1Dindex(timelabelsDat,myS) iendindx = get1Dindex(timelabelsDat,myE) ;print(istartindx + " " + iendindx) ;print(dimsizes(rainSpell)) tmp := rainSpell(istartindx:iendindx,:,:) ;printMinMax(tmp, 0) dryd = where(tmp .le. rainzero, 1,0) wetd = where(tmp .gt. heavyrain,1,0) drysp = dim_numrun_n(dryd,1, 0) ;unique dryday runs, dim 0 wetsp = dim_numrun_n(wetd,1, 0) ;unique wetday runs, dim 0 dryMonAll(yrind,mymonth, imodel, iens, :,:) = (/maxrunlen(drysp) /) wetMonAll(yrind,mymonth, imodel, iens, :,:) = (/maxrunlen(wetsp) /) delete([/istartindx, iendindx, iendday /]) delete([/dryd,wetd,drysp,wetsp/]) mymonth = mymonth + 1 end do ;;;;;;;;;;;;;;; ;;;; Continue with onset stuff ;;;;;;;;;;;;;;; startindx = 0 ;from here on the 0 index is the start search for onset endindx = indxOnsetEnd - indxOnsetStart ; raindaily = where(raindaily .lt. 0.1, 0.0, raindaily) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Define wet and dry days and maximum dry and wet spells ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; These are used for onset cessation wetdays := new(dimsizes(raindaily), integer, -999) drydays := wetdays wet5day := wetdays wetdays = where(raindaily .gt. rainzero, 1, 0) ; 1 for rainy day drydays = where(raindaily .le. rainzero, 1, 0) ; 1 for no rain wet5day = where(raindaily .gt. heavyrain, 1, 0) ; 1 for wet day ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;; ONSET/CESSATION for each ensemble member xonset = raindaily(0,:,:) xonset@_FillValue = -9999 xonset = xonset@_FillValue xcessn = xonset xgrowlen = xonset xdryspelAc = xonset xdryloc0Ac = xonset xdryloc1Ac = xonset xwetspelAc = xonset xwetloc0Ac = xonset xwetloc1Ac = xonset xdryspelFx = xonset xdryloc0Fx = xonset xdryloc1Fx = xonset xwetspelFx = xonset xwetloc0Fx = xonset xwetloc1Fx = xonset ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;; These are used for climatological max dry/wet spells ;;;;;; But will be overwritten if max dry/wet spell are computed ;;;;;; bewteen onset and cessation dates x := drydays(startindx:endindx,:,:) ; y := xRunTimeLen(x) xdryspelFx = (/ y[0] /) xdryloc0Fx = (/ y[1] /) xdryloc1Fx = (/ y[2] /) x := wet5day(startindx:endindx,:,:) ; y := xRunTimeLen(x) xwetspelFx = (/ y[0] /) xwetloc0Fx = (/ y[1] /) xwetloc1Fx = (/ y[2] /) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Main block for determining onset,cessation ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; izcount = 0 do j=jlatStrt, jlatLast do i=ilonStrt,ilonLast onsetfound = False do mwetstart=startindx, endindx ;start at 1st day of season mwetend = mwetstart + numraindays -1 if (mwetend .le. endindx ) then ; is the next 3 days contained rainydaysum = dim_sum_n(wetdays(mwetstart:mwetend,j,i),0) rainsumwet = dim_sum_n(raindaily(mwetstart:mwetend,j,i),0) ;print(rainsumwet + " rw nd " + rainydaysum) ndsum := rainydaysum rwsum := rainsumwet if (.not.ismissing(rwsum) .and. .not.ismissing(ndsum) .and. \ rainydaysum .ge. ndayrain .and. rainsumwet .ge. \ rtrsholdstart .and. .not. onsetfound) then itruestart = mwetstart ;- ndsum mwetendp1 = mwetend + 1 mdry = mwetendp1 + ndayafteronset - 1 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Now check if the rain does not stop after mwetend ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if (mdry .le. endindx) then x := drydays(mwetendp1:mdry,j,i) raindayafter = maxrunlen(dim_numrun_n(x,1, 0)) delete(x) if (.not. ismissing(raindayafter) .and. raindayafter .le. limnorainydayafter ) then onsetdate = 1.0*itruestart xonset(j, i) = onsetdate ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Cessation 1 starts here ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; icessnstart = toint(onsetdate) tmpraindaily = raindaily(:,j,i) tmpdrydays = drydays(:,j,i) itruecessn = CessationFun(icessnstart, \ cessndrylen, cessanrainamount, \ endindx,tmpraindaily, \ tmpdrydays) xcessn(j,i) = 1.0*itruecessn delete([/tmpraindaily,tmpdrydays/]) ;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;; 15 Aug 2019 -- Part I: What is the dry length immediately ;;;;;;;;;; following the onset but before cessation but ;;;;;;;;;; after the wet spell run breaks ;;;;;;;;;;ZTS: 29 Jan2022 -- Removed dryspells timing following onset ;;;;;;;;;; becasue usually a few days long. I now replaced ;;;;;;;;;; by calculating timing of maximum dry len afteronset ;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Now compute dryspell and maximum wet spell between the two dates ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; istart = toint(xonset(j,i)) iend = toint(xcessn(j,i)) ;we can assign iend = endindx if missing if (.not.ismissing(istart) .and. \ .not. ismissing(iend)) then x := drydays(istart:iend,j,i) y := xRunTimeLen(x) xdryspelAc(j,i) = (/ y[0] /) xdryloc0Ac(j,i) = (/ istart + y[1] /) xdryloc1Ac(j,i) = (/ istart + y[2] /) x := wet5day(istart:iend,j,i) y := xRunTimeLen(x) xwetspelAc(j,i) = (/ y[0] /) xwetloc0Ac(j,i) = (/ istart + y[1] /) xwetloc1Ac(j,i) = (/ istart + y[2] /) delete([/x,y/]) end if delete([/istart, iend/]) ;;;;;;;;;;;;;;;;;;;;;;;; ;; Compute growing length ;;;;;;;;;;;;;;;;;;;;;;;;; istart = toint(xonset(j,i)) iend = toint(xcessn(j,i)) if (.not.ismissing(istart) .and. \ .not. ismissing(iend)) then xgrowlen(j,i) = iend - istart +1 else xgrowlen(j,i) = xgrowlen@_FillValue end if ; end growing length delete([/istart, iend/]) onsetfound = True delete(onsetdate) break end if ;onset is satisfied (rainafter le limno..) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;This is line 434;;;;;;;;;;;;;;;; else ; mdry .le. endindx #366 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; mdry is not endindx, but this could be becasue of ; large ndayafteronset for regions of late onset ; that get their rain in December. Let us atleast try ; one more time by reducing ndayafteronset ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; mdry = itruestart + ndayafteronset - 10 - 1 if (mdry .le. endindx) then ; x := drydays(mwetendp1:mdry,j,i) raindayafter = maxrunlen(dim_numrun_n(x,1, 0)) delete(x) if (.not. ismissing(raindayafter) .and. raindayafter .le. limnorainydayafter) then izcount = izcount + 1 onsetdate = 1.0*itruestart xonset(j, i) = onsetdate onsetfound = True ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Cessation 2 starts here ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; icessnstart = toint(onsetdate) tmpraindaily = raindaily(:,j,i) tmpdrydays = drydays(:,j,i) itruecessn = CessationFun(icessnstart,\ cessndrylen, cessanrainamount, \ endindx, tmpraindaily, \ tmpdrydays) xcessn(j,i) = 1.0*itruecessn delete([/tmpraindaily,tmpdrydays/]) ;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;; 15 Aug 2019 -- Part II: What is the dry length immediately ;;;;;;;;;; following the onset but before cessation but ;;;;;;;;;; after the wet spell run breaks ;;;;;;;;;;ZTS: 29 Jan2022 -- Removed dryspells timing following onset ;;;;;;;;;; becasue usually a few days long. I now replaced ;;;;;;;;;; by calculating timing of maximum dry len afteronset ;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Now compute dryspell and maximum wet spell between the two dates ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; istart = toint(xonset(j,i)) iend = toint(xcessn(j,i)) ;we can assign iend = endindx if missing if (.not.ismissing(istart) .and. \ .not. ismissing(iend)) then x := drydays(istart:iend,j,i) y := xRunTimeLen(x) xdryspelAc(j,i) = (/ y[0] /) xdryloc0Ac(j,i) = (/istart + y[1]/) xdryloc1Ac(j,i) = (/istart + y[2]/) x := wet5day(istart:iend,j,i) y := xRunTimeLen(x) xwetspelAc(j,i) = (/ y[0] /) xwetloc0Ac(j,i) = (/istart + y[1]/) xwetloc1Ac(j,i) = (/istart + y[2]/) delete([/x,y/]) end if delete([/istart, iend/]) ;;;;;;;;;;;;;;;;;;;;;;;; ;; Compute growing length ;;;;;;;;;;;;;;;;;;;;;;;;; istart = toint(xonset(j,i)) iend = toint(xcessn(j,i)) if (.not.ismissing(istart) .and. \ .not. ismissing(iend)) then xgrowlen(j,i) = iend - istart +1 else xgrowlen(j,i) = xgrowlen@_FillValue end if ; end growing length delete([/istart, iend/]) break ; this is line 509 end if ;~447; onset is satisfied (raindayafter .le. limnorainydayafter ...) else xonset(j,i) = xonset@_FillValue end if ; ~443/369; reduced mdry one more time end if ; ~ 435/365 mdry .le. endindx; #514 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;print("j/i "+xlat(j,i)+" "+xonset(j,i)+" "+xlon(j,i) + " RR "+rainsumwet+" "+rainydaysum) ;> 3 days & > 25mm but no wet after") ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; end if ;~ 357; First onset condition completed; rainydaysum .ge. ndayrain .and. else ; ~353; mwetend .le. endindx; This is #522 xonset(j,i) = xonset@_FillValue end if ; 353 mend le endindx; this is #522 end do ;351 mwetstart=startindx, endindx end do ;~349 i=ilonStrt,ilonLast end do ;~348 j=jlatStrt, jlatLast ; this is # 527 print("Yr: " + iyear +", Ens: " + iens + " +21 days took it over "+timelabels(endindx)+" for "+izcount+" points ") ;;;;;;; Save data onsetAll(yrind, imodel,iens,:,:) = (/ tofloat(xonset) /) cesationAll(yrind, imodel,iens,:,:) = (/ tofloat(xcessn) /) growlengAll(yrind, imodel,iens,:,:) = (/ tofloat(xgrowlen) /) wetspLengAct(yrind,imodel,iens,:,:) = (/ tofloat(xwetspelAc) /) wetspLoc0Act(yrind,imodel,iens,:,:) = (/ tofloat(xwetloc0Ac) /) wetspLoc1Act(yrind,imodel,iens,:,:) = (/ tofloat(xwetloc1Ac) /) dryspLengAct(yrind,imodel,iens,:,:) = (/ tofloat(xdryspelAc) /) dryspLoc0Act(yrind,imodel,iens,:,:) = (/ tofloat(xdryloc0Ac) /) dryspLoc1Act(yrind,imodel,iens,:,:) = (/ tofloat(xdryloc1Ac) /) wetspLengFix(yrind,imodel,iens,:,:) = (/ tofloat(xwetspelFx) /) wetspLoc0Fix(yrind,imodel,iens,:,:) = (/ tofloat(xwetloc0Fx) /) wetspLoc1Fix(yrind,imodel,iens,:,:) = (/ tofloat(xwetloc1Fx) /) dryspLengFix(yrind,imodel,iens,:,:) = (/ tofloat(xdryspelFx) /) dryspLoc0Fix(yrind,imodel,iens,:,:) = (/ tofloat(xdryloc0Fx) /) dryspLoc1Fix(yrind,imodel,iens,:,:) = (/ tofloat(xdryloc1Fx) /) delete([/wetdays,drydays,wet5day,xonset,xcessn, raindaily/]) delete([/xgrowlen,startindx,endindx,xwetspelAc,xwetspelFx/]) delete([/xwetloc0Ac,xwetloc1Ac,xdryspelAc,xdryloc0Ac,xdryloc1Ac/]) delete([/xwetloc0Fx,xwetloc1Fx,xdryspelFx,xdryloc0Fx,xdryloc1Fx/]) end do ; ens ; ZTS 10May2020 end do ; iyr ; 0 to hindYrs ZTS: 19 Jan 2021 end do ; imodel ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if (writeoutany) then ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;; Monthly dry and wet spells ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; fname := analdir+"AllModels"+"_HindcastMonthlyDryWetSpells"+SEASON+".nc" system("/bin/rm -f " + fname) fout = addfile(fname,"c") print(fname) ;*************************************************** ; define file options [Version a033] ;*********************************************** setfileoption(fout, "prefill", False) setfileoption(fout, "suppressclose", True) setfileoption(fout, "definemode", True) ;=================================================================== ; create global attributes of the file (optional) ;=================================================================== fAtt = True ; assign file attributes fAtt@title = SEASON+" Monthly dry/wet spell" fAtt@source_file = " CMORPH " fAtt@Conventions = "None" fAtt@creation_date = systemfunc ("date") fAtt@start_month = startmon fAtt@Last_month = endmon fAtt@heavy_rain = heavyrain + " mm/day" fAtt@rain_zero = rainzero + " mm/day" fileattdef( fout, fAtt ) ; copy file attributes ;*********************************************** ; predefine coordinate information ;*********************************************** dimNames = (/"time", "month","models","ens", "lat", "lon" /) dimSizes = (/ nyears ,nummon, nmodels, ensmax, mlat, mlon /) dimUnlim = (/ True , False, False, False, False, False /) filedimdef(fout, dimNames , dimSizes, dimUnlim ) filevardef(fout, "time" , typeof(time), getvardims(time) ) filevarattdef(fout, "time", time) filevardef(fout, "month" , typeof(dmonth), getvardims(dmonth) ) filevarattdef(fout, "month", dmonth) filevardef(fout, "models" , typeof(models) , getvardims(models) ) filevarattdef(fout, "models", models) filevardef(fout, "ens" , typeof(ensall) , getvardims(ensall) ) filevarattdef(fout, "ens", ensall) filevardef(fout, "lat", typeof(REGlat), getvardims(REGlat) ) filevarattdef(fout, "lat", REGlat) filevardef(fout, "lon", typeof(REGlon), getvardims(REGlon) ) filevarattdef(fout, "lon", REGlon) ;*********************************************** ; predefine variable sizes ;*********************************************** filevardef(fout,"TIME",typeof(TIME), getvardims(TIME) ) filevarattdef(fout,"TIME",TIME) filevardef(fout,"dryspell",typeof(dryMonAll), getvardims(dryMonAll) ) filevarattdef(fout,"dryspell",dryMonAll) filevardef(fout,"wetspell",typeof(wetMonAll), getvardims(wetMonAll) ) filevarattdef(fout,"wetspell",wetMonAll) ;*********************************************** ; terminate define mode: not necessary but for clarity ;*********************************************** setfileoption(fout, "definemode", False) ;*********************************************** ; write data values to predefined locations ; (/ .../) operator transfer values only ;*********************************************** fout->time = (/ time /) fout->month = (/ dmonth /) fout->models = (/ models /) fout->ens = (/ ensall/) fout->lat = (/ REGlat /) fout->lon = (/ REGlon /) fout->TIME = (/ TIME /) fout->dryspell = (/ dryMonAll /) fout->wetspell = (/ wetMonAll /) delete(fout) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;; done writing monthly dry and wet spells ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; writeout onset/cessation/dryspell/wetspell in netcdf ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;tmpS = monthdef(startmon)+ "-" + startday ; onset search start day tmpS = monthdef(sMonOnset)+ "-" + startday ; onset search start day tmpE = monthdef(endmon) + "-" + endday xnames = "AllModels" fname = analdir+xnames+"_HindcastOnsetCessnDryWetGrowLen"+SEASON+".nc" system("/bin/rm -f " + fname) fout = addfile(fname,"c") ;*************************************************** ; define file options [Version a033] ;*********************************************** setfileoption(fout, "prefill", False) setfileoption(fout, "suppressclose", True) setfileoption(fout, "definemode", True) ;=================================================================== ; create global attributes of the file (optional) ;=================================================================== fAtt = True ; assign file attributes fAtt@title = SEASON+" Seasonal Characteristic " fAtt@source_file = "GCM daily forecasts " fAtt@Conventions = "None" fAtt@creation_date = systemfunc ("date") fAtt@start_month = sMonOnset fAtt@start_day = startday fAtt@heavy_rain = heavyrain + " mm/day" fAtt@rain_zero = rainzero + " mm/day" fAtt@numraindays = numraindays fAtt@ndayrainInNumr = ndayrain fAtt@rtrsholdstart = rtrsholdstart + " mm/day" fAtt@ndayafteronset = ndayafteronset fAtt@limnorainydayafter = limnorainydayafter fAtt@cessationRainThr = cessanrainamount + " mm/day" fAtt@icessationDrylen = cessndrylen fAtt@onsetStart = tmpS fAtt@onsetEnd = tmpE fAtt@onsetUnit = "days from " + tmpS fAtt@odrystart = tmpS fAtt@dayoffset = tmpS fAtt@nmodels = nmodels fAtt@modelNames = varx ; startday - 1 ; day_of_year(yyyy,startmon,startday) fileattdef( fout, fAtt ) ; copy file attributes ;*********************************************** ; predefine coordinate information ;*********************************************** dimNames := (/"time", "models","ens", "lat", "lon" /) dimSizes := (/ nyears ,nmodels, ensmax, mlat, mlon /) dimUnlim := (/ True , False,False, False, False /) filedimdef(fout, dimNames , dimSizes, dimUnlim ) filevardef(fout, "time" , typeof(time), getvardims(time) ) filevarattdef(fout, "time", time) filevardef(fout, "models" , typeof(models) , getvardims(models) ) filevarattdef(fout, "models", models) filevardef(fout, "ens" , typeof(ensall) , getvardims(ensall) ) filevarattdef(fout, "ens", ensall) filevardef(fout, "lat", typeof(REGlat), getvardims(REGlat) ) filevarattdef(fout, "lat", REGlat) filevardef(fout, "lon", typeof(REGlon), getvardims(REGlon) ) filevarattdef(fout, "lon", REGlon) ;*********************************************** ; predefine variable sizes ;*********************************************** filevardef(fout,"TIME",typeof(TIME), getvardims(TIME) ) filevarattdef(fout,"TIME",TIME) filevardef(fout,"onset",typeof(onsetAll), \ getvardims(onsetAll) ) filevarattdef(fout,"onset",onsetAll) filevardef(fout,"cessation",typeof(cesationAll), \ getvardims(cesationAll) ) filevarattdef(fout,"cessation",cesationAll) filevardef(fout,"growlength",typeof(growlengAll), \ getvardims(growlengAll)) filevarattdef(fout,"growlength",growlengAll) filevardef(fout,"wetspLengAct",typeof(wetspLengAct), \ getvardims(wetspLengAct)) filevarattdef(fout,"wetspLengAct",wetspLengAct) filevardef(fout,"dryspLengAct",typeof(dryspLengAct), \ getvardims(dryspLengAct)) filevarattdef(fout,"dryspLengAct",dryspLengAct) filevardef(fout,"wetspLoc0Act",typeof(wetspLoc0Act), \ getvardims(wetspLoc0Act)) filevarattdef(fout,"wetspLoc0Act",wetspLoc0Act) filevardef(fout,"dryspLoc0Act",typeof(dryspLoc0Act), \ getvardims(dryspLoc0Act)) filevarattdef(fout,"dryspLoc0Act",dryspLoc0Act) filevardef(fout,"wetspLoc1Act",typeof(wetspLoc1Act),\ getvardims(wetspLoc1Act)) filevarattdef(fout,"wetspLoc1Act",wetspLoc1Act) filevardef(fout,"dryspLoc1Act",typeof(dryspLoc1Act),\ getvardims(dryspLoc1Act)) filevarattdef(fout,"dryspLoc1Act",dryspLoc1Act) filevardef(fout,"wetspLengFix",typeof(wetspLengFix), \ getvardims(wetspLengFix)) filevarattdef(fout,"wetspLengFix",wetspLengFix) filevardef(fout,"dryspLengFix",typeof(dryspLengFix), \ getvardims(dryspLengFix)) filevarattdef(fout,"dryspLengFix",dryspLengFix) filevardef(fout,"wetspLoc0Fix",typeof(wetspLoc0Fix), \ getvardims(wetspLoc0Fix)) filevarattdef(fout,"wetspLoc0Fix",wetspLoc0Fix) filevardef(fout,"dryspLoc0Fix",typeof(dryspLoc0Fix), \ getvardims(dryspLoc0Fix)) filevarattdef(fout,"dryspLoc0Fix",dryspLoc0Fix) filevardef(fout,"wetspLoc1Fix",typeof(wetspLoc1Fix), \ getvardims(wetspLoc1Fix)) filevarattdef(fout,"wetspLoc1Fix",wetspLoc1Fix) filevardef(fout,"dryspLoc1Fix",typeof(dryspLoc1Fix), \ getvardims(dryspLoc1Fix)) filevarattdef(fout,"dryspLoc1Fix",dryspLoc1Fix) ;*********************************************** ; terminate define mode: not necessary but for clarity ;*********************************************** setfileoption(fout, "definemode", False) ;*********************************************** ; write data values to predefined locations ; (/ .../) operator transfer values only ;*********************************************** fout->time = (/ time /) fout->models = (/ models /) fout->ens = (/ ensall/) fout->lat = (/ REGlat /) fout->lon = (/ REGlon /) fout->TIME = (/ TIME /) fout->onset = (/ onsetAll /) fout->cessation = (/ cesationAll /) fout->growlength = (/ growlengAll /) fout->wetspLengAct = (/ wetspLengAct /) fout->dryspLengAct = (/ dryspLengAct /) fout->wetspLoc0Act = (/ wetspLoc0Act /) fout->dryspLoc0Act = (/ dryspLoc0Act /) fout->wetspLoc1Act = (/ wetspLoc1Act /) fout->dryspLoc1Act = (/ dryspLoc1Act /) fout->wetspLengFix = (/ wetspLengFix /) fout->dryspLengFix = (/ dryspLengFix /) fout->wetspLoc0Fix = (/ wetspLoc0Fix /) fout->dryspLoc0Fix = (/ dryspLoc0Fix /) fout->wetspLoc1Fix = (/ wetspLoc1Fix /) fout->dryspLoc1Fix = (/ dryspLoc1Fix /) delete(fout) ;;;;;;;; imodel = 0 MODEL = str_upper(modelNames(imodel))+ "v" +version(imodel) fname = fDatDir + "precDaily"+ initMonth+"Hindcast" + MODEL +".nc" a := addfile(fname, "r") ztemp := a->time ; yyyymmddhh; initial time in the Hindcast variable leadAll := a->lead ; Offset (in da intime := ztemp ; hindYrs := intime/1000000 mmddhh := intime - hindYrs*1000000 allmons := mmddhh/10000 ddhh := mmddhh - allmons*10000 allday := ddhh/100 allhrs := ddhh - allday*100 ;;;;;;;;;;;;;;; ;;;;;;; Loop over all years (initial time) ;;;;;; iyr = dimsizes(hindYrs)-1 iyear = hindYrs(iyr) imonth = allmons(iyr) iday = allday(iyr) ihr = allhrs(iyr) lead := leadAll ;a->lead--Offset (in days) from each initial time lead!0 = "time" lead&time = leadAll x1 = "days since "+iyear+ "-" + sprinti("%0.2i",imonth) x2 = "-" + sprinti("%0.2i",iday) + " " x3 = sprinti("%0.2i",ihr) + ":00:0" lead@units = x1 + x2 + x3 lead = (/ lead -1 /) utc_date := cd_calendar(lead, 0) ;gives 6 col for yr,MM,DD,HH,mm,ss fmonth := tointeger(utc_date(:,1)) ; all months of data for each day fday := tointeger(utc_date(:,2)) ; all days of data for a given year timelabelsDat := monthdef(fmonth)+"-"+fday ; full forecast time ;tmpS = monthdef(startmon)+ "-" + startday ; onset search start day tmpS = monthdef(sMonOnset)+ "-" + startday ; onset search start day tmpE = monthdef(endmon) + "-" + endday ; onset search last day indxOnsetStart = get1Dindex(timelabelsDat,tmpS) indxOnsetEnd = get1Dindex(timelabelsDat,tmpE) timelabels := timelabelsDat(indxOnsetStart:indxOnsetEnd) fout = analdir+ xnames +"OnsetCessnDatesClim"+SEASON+".txt" system("/bin/rm -f "+ fout) asciiwrite(fout, timelabels) delete(fout) end if ; if writeoutany end